Mathematical Model of a Hyperbolic Hydraulic Fracture with Tortuosity

The aim of the research is to study the propagation of a hydraulic fracture with tortuosity due to contact areas between touching asperities on opposite crack walls. The tortuous fracture is replaced by a model symmetric partially open fracture with a hyperbolic crack law and a modified Reynolds flow law. The normal stress at the crack walls is assumed to be proportional to the half-width of the model fracture. The Lie point symmetry of the nonlinear diffusion equation for the fracture half-width is derived and the general form of the group invariant solution is obtained. It was found that the fluid flux at the fracture entry cannot be prescribed arbitrarily, because it is determined by the group invariant solution and that the exponent n in the modified Reynolds flow power law must lie in the range 2 5 n < < . The boundary value problem is solved numerically using a backward shooting method from the fracture tip, offset by 0 1 δ <  to avoid singularities, to the fracture entry. The numerical results showed that the tortuosity and the pressure due to the contact regions both have the effect of increasing the fracture length. The spatial gradient of the half-width was found to be singular at the fracture tip for 3 5 n < < , to be finite for the Reynolds flow law 3 n = and to be zero for 2 3 n < < . The thin fluid film approximation breaks down at the fracture tip for 3 5 n < < while it remains valid for increasingly tortuous fractures with 2 3 n < < . The effect of the touching asperities is to decrease the width averaged fluid velocity. An approximate analytical solution for the half-width, which was found to agree well with the numerical solution, is derived by making the approximation that the width averaged fluid velocity increases linearly with distance along the fracture.

the group invariant solution and that the exponent n in the modified Reynolds flow power law must lie in the range 2 5 n < < . The boundary value problem is solved numerically using a backward shooting method from the fracture tip, offset by 0 1 δ <  to avoid singularities, to the fracture entry. The numerical results showed that the tortuosity and the pressure due to the contact regions both have the effect of increasing the fracture length. The spatial gradient of the half-width was found to be singular at the fracture tip for 3 5 n < < , to be finite for the Reynolds flow law 3 n = and to be zero for 2 3 n < < . The thin fluid film approximation breaks down at the fracture tip for 3 5 n < < while it remains valid for increasingly tortuous fractures with 2 3 n < < . The effect of the touching asperities is to decrease the width averaged fluid velocity. An approximate analytical solution for the half-width, which was found to agree well with the numerical solution, is derived by making the approximation that the width averaged fluid velocity increases linearly with distance along the fracture.

Introduction
In hydraulic fracturing, fluid is pumped at high pressure into a crack in a rock mass in order to open the crack. Hydraulic fracturing has many applications. In mining, high pressure water is used to break rocks instead of explosives which leave small rock fragments in the air that could damage the lungs of miners [1].
In geothermal reservoirs, in their natural state, the cracks allow only a small flow of water. In order to increase the flow, high pressure water is pumped into the crack network at one borehole and extracted from another borehole [2]. Hydraulic fracturing is also used to enhance the extraction of oil and gas in large underground shale deposits [3]. Spence and Sharp [4] were the first to show that the equations of hydraulic fracture admit a similarity solution. These authors considered the enlargement of a lens shaped crack and a two-dimensional crack by a viscous fluid modelled by lubrication theory. The fluid pressure and the crack shape were connected by a singular integral equation from linear elasticity. The theory was applied to magma-driven propagation of cracks in geophysics [5] [6] [7].
The Cautchy principal value integral in the singular integral equation relating pressure to the crack shape is difficult to analyze. An important simplification was the Perkins-Kern-Nordgren (PKN) approximation in which the normal stress at the fluid-rock interface is proportional to the half-width of the fracture [8] [9]. The PKN approximation puts the differential equations of hydraulic fracture in a form which can be analyzed using the theory of Lie point symmentries and conservation laws.
The scientific literature on hydraulic fracturing is now very large. We will only comment briefly on the application of Lie point symmetries to hydraulic fracturing, because that approach will be used in this paper. Lie point symmetry analysis was first applied to hydraulic fracturing by Fitt, Mason and Moss who considered the propagation of a fracture in impermeable rock [1]. Fareo and Mason analyzed the propagation of a hydraulic fracture in permeable rock with fluid leak-off into the rock mass and investigated the effect of leak-off on the rate of propagation of a fracture [10]. Fareo and Mason also considered a fracture driven by a power law fluid in an impermeable rock mass [11]. Anthonyrajah, Mason and Fareo [12] compared laminar and turbulent fluid driven fractures using the wall shear stress model of Emerman, Turcotte and Spence [6] for the fluid and the PKN approximation instead of the Cautchy principal value integral model to relate the fluid pressure to the crack shape.
The effects of tortuosity on hydraulic fracturing due to asperities or surface roughness at the fluid rock interface and contact regions caused by touching asperities will be investigated in this paper. The hydraulic fracture with tortuosity will be replaced by a symmetric two-dimensional hydraulic fracture without asperities but with a modified Reynolds flow law, which accounts for the effect that the presence of asperities at the fluid rock interface has on the fluid flow, and Journal of Applied Mathematics and Physics with a modified crack law, which models the effect of the contact regions on the stress at the fluid-rock interface [2].
In this paper, we model the contact regions by the hyperbolic crack law which was introduced by Goodman [13]. This is motivated by the discussion by Fitt et al. [2] that the hyperbolic crack law is generally considered to be a more realistic model to describe the presence of contact regions (deformations formed by touching asperities) in a fracture than the linear crack law. Although the concept of tortuosity was investigated in [14] [15], the contact regions in these papers were modelled by the linear crack law which was proposed by Pine and Cundall [16] and further discussed by Fitt et al. [2].
In this paper, we aim to solve the problem of the evolution of a hydraulic fracture with tortuosity described by the hyperbolic crack law. The fluid flow in the fracture is described by a modified Reynolds flow law and the fluid pressure and crack shape are related by the PKN approximation. The aim of this work is to derive numerical and approximate analytical solutions for the evolution of the half-width and length of the model symmetric fracture which replaces the fracture with tortuosity. The problem is formulated mathematically and in dimensionless form in Section 2.6.
In Section 2, we present the derivation of the model describing fluid flow in a tortuous hydraulic fracture with contact regions modelled by the hyperbolic crack law. Section 3 outlines the derivation of the group invariant solution of the problem. In Section 4 we consider possible operating conditions at the fracture entry and investigate the existence of analytical solutions. Section 5 describes a numerical solution for the problem. In Section 6, the investigation of the width averaged fluid velocity leads to the derivation of an approximate analytical solution for the problem. In this Section, a comparison of the numerical and the approximate analytical solutions is made. Finally in Section 7 we summarize important findings and conclude the paper. Journal of Applied Mathematics and Physics

( )
, h t x . Since the fracture is long and thin, the lubrication approximation is imposed on the viscous flow in the fracture which gives that the fluid pressure in independent of the variable z. It follows that the fluid velocity components and the fluid pressure respectively have the form The fluid is incompressible with constant density ρ and dynamic viscosity µ . The body force due to gravity is neglected.

Review of the Flow Model
From lubrication theory for a two-dimensional hydraulic fracture the volume flux of fluid across the fracture per unit breadth, ( ) , Q t x , satisfies the Reynolds flow law [14]: and the width averaged fluid velocity is The half-width of the fracture ( ) which when expanded with the aid of (2.2) gives the nonlinear diffusion equation relating half-width of the fracture ( ) , h t x and the fluid pressure ( ) and

Hyperbolic Crack Law
In [14] the linear crack law was used to describe the contribution made by the contact areas to support the compressive normal stress, , at the fracture walls. In this paper we will use the hyperbolic crack law which was first introduced by Goodman [13] and further considered by Bandis et al. [21], Murphy et al. [22] and Fitt et al. [2].  [14] and [15], in this paper we will discuss only a partially open fracture. We consider the hyperbolic crack law [13], where k is a negative constant, min h is the minimum half-width of the fracture and max h is the maximum half-width of the fracture. We make the approximation that min 0 h = which simplifies (2.11) to which clearly shows that the compressive normal stress at the fluid-rock interface, , is supported by both the fluid pressure ( ) , p t x and the pressure due to contact regions ( ) 1 , p t x .

The PKN Approximation
Lastly, it is necessary to specify a relation between the compressive normal stress at the fluid-rock interface, where E is the Young's modulus and ν is the Poisson ratio of the rock encasing the fracture and B is the breadth of the fracture. The PKN approximation has the disadvantage that the stress intensity factor vanishes at the fracture tip, ( ) x L t = , because the half-width is zero at the tip but it can be applied to a single-sided fracture with non-zero initial length. This is unlike the relation,

Model Closure
From equations (2.13) and (2.15), the fluid pressure in the fracture is and therefore Since 0 k < the magnitude of the pressure gradient in the fluid along the fracture is increased by the contact regions. Equations (2.6) and (2.7) for the volume flux of fluid in the fracture and the width averaged fluid velocity respec- and (2.8) becomes the nonlinear diffusion equation for the half-width of the fracture, ( ) for an open fracture which was discussed in [14] and [15]. In this paper we consider only partially open fractures ( 0 k < ) for which the compressive normal stress at the fracture walls is supported by both the fluid pressure ( ) , p t x and the pressure ( ) 1 , p t x due to contact regions given by (2.14).

Equation (2.22) is a nonlinear diffusion equation of the form
where the diffusion coefficient but with diffusion coefficients different from (2.24) [25]. To the best of our knowledge there is no application to hydraulic fracturing with a diffusion coefficient of the form (2.24).

Boundary and Initial Conditions
Consider first the boundary conditions at the fracture tip, ( ) The second boundary condition at the fracture tip is that the volume flux of fluid vanishes: and therefore using (2.20), The boundary conditions (2.25) and (2.27) are moving boundary conditions [26].
Consider next the condition at the fracture entry, 0 x = . We consider a partially open fracture. Let The parameter β is prescribed and describes the extent to which the initial fracture is partially open. The boundary condition at the fracture entry is = is the flux of fluid per unit breadth into the fracture at the entry 0 x = . It is prescribed as far as the invariant solution will allow.

Dimensionless Governing Equations
In order to make all equations and boundary and initial conditions dimensionless we transform to dimensionless variables. Consider first the fluid pressure.

Equation (2.18) can be written as
, P t x is the difference between the pressure of the fluid in the fracture, ( ) , p t x , and the background pressure, 3) to (2.8) on which the theory is based, remain unchanged when expressed in terms of ( ) We choose characteristic quantities which are independent of n and n a so that the results for a range of values of n can be compared. We choose max h as the characteristic half-width. From (2.30) the characteristic value of ( ) The characteristic value of the x-coordinate is the initial length of the Journal of Applied Mathematics and Physics which is also the characteristic length of the fracture. Equation (2.5) was derived using lubrication theory and gives for the characteristic velocity along the fracture: We transform to the following dimensionless variables: is the total volume of the fracture per unit breadth and ( ) We suppress the stars on the dimensionless variables and parameters to keep the notation simple. Expressed in dimensionless form the half-width of the fracture, ( ) , h t x , satisfies the nonlinear diffusion equation subject to the boundary conditions and the initial/boundary condition In dimensionless variables, the volume flux of fluid in the fracture, the width averaged fluid velocity, the fluid pressure and the volume of the fracture are None of the characteristic quantities depend on n or n a . The solution depends on the two parameters n K and φ and on ( ) S t , the fluid source at the fracture entry. The parameter n K contains all the dependence on n and n a .
The physical significance of φ is pressure due to contact regions compressive normal stress at the fracture walls minus the background compressive stress φ = (2.46) This completes the formulation of the mathematical model of a tortuous fracture with the hyperbolic crack law.

Group Invariant Solution
We consider the propagation of a pre-existing hydraulic fracture with a non-zero initial length. Therefore methods used to derive a similarity solution for gravity currents [27] and for hydraulic fractures [4], both with a zero initial length, are not applicable in this work. Lie group analysis of partial differential equations, which can be used to analyze fractures with non-zero initial length, will be used to solve the problem [14] [10]. We derive the Lie point symmetries of the nonlinear diffusion Equation (2.37) which are then used to reduce the problem to a boundary value problem for an ordinary differential equation which describes the effects of surface roughness and contact regions on the evolution of a two-dimensional hydraulic fracture.

Lie Point Symmetries
The form of the Lie point symmetry generator of the nonlinear diffusion equa- which satisfies the determining equation where the subscripts t and x denote partial differentiation with respect to t and x and [ ] is the second prolongation of the generator X with [28] ( ) ( ) , 1, 2,

) Journal of Applied Mathematics and Physics
There is summation over the repeated index k from 1 to 2 and We find that the Lie point symmetry generator X for the non-linear diffusion c c c are constants. We found that the coefficient ( ) This is the first significant consequence of the from (2.24) for the diffusion coefficient.
The group invariant solution, h = Ψ , of the partial differential Equation which expands to the first order partial differential equation Consider the general case for which 2 0 c ≠ . By solving the differential equations of the characteristic curves of (3.10) the general solution is readily derived and since ( ) where ( ) F ξ is an arbitrary function of the similarity variable ξ given by .
Unlike the group invariant solution for the linear fracture law, ( ) , h t x depends only on the similarity variable ξ and does not depend on an arbitrary parameter. Substituting the solution (3.11) and (3.12) into (2.37) reduces the second order partial differential equation to a second order ordinary differential Since the constant 3 c does not appear in the ordinary differential Equation We consider the general case for which where o w is a constant. It follows that ( ) which can be written as ( ) where we have assumed that 1 0 c ≠ . Since the characteristic length of the fracture is the initial length of the fracture, ( ) 0 1 L = and therefore and the boundary condition at the fracture tip (3.15) becomes Expressed in terms of ( ) We see that the time dependence of fluid flux into the fracture at the entry,

( )
S t , cannot be prescribed arbitrarily but must be of the form ( We simplify the formulation by transforming from variables ξ and ( ) F ξ to variables u and ( ) We also define The boundary value problem transforms to the ordinary differential equation subject to the boundary conditions ( ) where the fluid flux source at the fracture entry is ( ) The remaining physical variables become ( )  The Lie point symmetry which generates the invariant solution is ( ) The physical significance of the source (3.42) is that the fluid pressure remains constant at the fracture entry for all time 0 t ≥ . From (3.47) and the boundary condition (3.38), From (3.44) and (3.38), it is clear that the half-width of the fracture also remains constant at the fracture entry for all time 0 t ≥ :

Asymptotic Solution
We now derive the asymptotic solution of the ordinary differential Equation The dominant terms as 1 u → are the terms with the smallest values for the exponent of ( ) 1 u − . To obtain an asymptotic solution these terms must balance and the remaining terms must vanish as 1 u → . The dominant terms in (4.2) are the second and the fourth terms with exponents ( ) The special case of 2 n = will be examined later. For 0 p > we require 2 n > . Equating coefficients of the second and fourth terms gives for 2 n > . The second term will vanish as 1 u → since 2 n > but the first term will vanish only if 2 5 n < < . It follows that the asymptotic solution for the ordinary differential Equation The asymptotic behaviour close to the fracture tip described by (4.6) is valid for all solutions of the hyperbolic hydraulic fracture described in this work. Since it follows that the boundary condition (3.41) at the fracture tip, 1 u = , is identically satisfied.
Consider now the special case of 2 n = . When The second term can no longer be balanced with the fourth term and it cannot be balanced with the first term for 0 p > . We balance the first term with the fourth term to give 1 2 p = . Equation (4.8) becomes Since we are considering 0 B ≠ , condition (4.9) requires Since the asymptotic solution for the hyperbolic crack law exists only for 2 5 n < < we conjecture that the solution of the boundary value problem, (3.37) to (3.41), exists only for 2 5 n < < . This conjecture will be investigated further when the numerical solution is derived in Section 5.
Consider now the behaviour of d d f u as 1 u → . From (4.6) for 2 5 n < < , and therefore The asymptotic results, (4.6) and (4.13), will be used in the numerical solution in Section 5 to offset the boundary condition at as ( ) x L t → , provided 2 5 n < < . Differentiating (4.14) with respect to x and evaluating the limit of the result as ( ) It is clear that the spatial gradient of the half-width of a two-dimensional hyperbolic model fracture is singular at ( ) x L t = for 3 5 n < < and the thin fluid film approximation breaks down at the fracture tip. However, we see that for the hyperbolic crack law, the spatial gradient of the half-width at the fracture tip is not singular for 2 3 n < ≤ . In the elasticity problem for an infinite body with a plane cut subject to symmetric loads Barenblatt describes similar behaviour near the cut tip [29].

Numerical Solution
In this Section, we present the numerical solution for the boundary value problem (3.37) to (3.42). The boundary value problem is solved with the backward shooting method from the neigbourhood of the fracture tip, 1 u δ = − where 1 δ  , to the fracture entry, 0 u = . Formulating a numerical scheme so that its step-size is adjusted according to the behaviour of the solution and the resulting truncating error in each step can greatly improve the accuracy of numerical solutions [30]. The boundary value problem presented in this work was solved using the Matlab ode45 solver, which applies to the problem the Runge Kutta method of order 4 and 5 with adaptive step-size [31]. A brief algorithm that shows how the problem is solved is provided below: In steps 3 and 4 of the algorithm the second order ODE is written as the system of two first order ordinary differential equations. Due to the PKN approximation the half-width of the model fracture is also constant for 0 t ≥ at the fracture entry and is given by (3.51). This is clearly shown in Figure Figure 3 and for a range of values of φ in Figure 4. The parameter β describes the extent to which the fracture si partially open while φ compares the contribution of the pressure due to the contact regions with the pressure due to the fluid in the fracture. From Figure 3 we observe that for a prescribed value of φ as the value of β increases, the fracture becomes longer. From Figure 3 and Figure 4 it is clear that for fixed values of β and φ , tortuosity increases, that is as n decreases from 5 to 2, the length of the model fracture increases. We also see from Figure 4 that there is less change in the length due to change in φ as the tortuosity increases.  Table 1. We see from Table 1 that as n decreases, that is as tortuosity increases, for fixed values of β and φ , the constant A increases. Thus a stronger fluid flux at the fracture    It is clear that the constant B determines the magnitude of the speed of propagation of the fracture. From Table 1, for fixed values of β and φ , the constant B decreases as n decreases. Thus the speed of propagation of the fracture increases as tortuosity of the fracture increases. This agrees with Figure 3 in which ( ) L t is plotted against time t for a range of values of n, β and φ . We also see from Table 1 Figure   3. For fixed values of β and n, A decreases as φ decreases and therefore the fluid flux at the entry is weaker as the pressure ratio φ decreases, consistent with Figure 5. Also for fixed values of β and n, B increases as φ decreases and therefore the speed of propagation of the fracture decreases as φ decreases, consistent with Figure 4 for ( ) L t .

Width Averaged Fluid Velocity
The ratio of the width averaged fluid velocity, x v , given by (3.46) to the speed of propagation of the fracture tip, ( ) Using the asymptotic solutions    Figure 6 that the velocity ratio decreases as φ increases, that is, as the effect of the pressure due to the contact regions increases. The effect of touching asperities is therefore to decrease the width averaged fluid velocity along the length of the tortuous fracture. We also see from Figure 6 that the velocity ratio linearly with u along the fracture for the values of the parameters considered.
The linear approximation is best for 4 n = and less applicable as n decreases and the fracture becomes more tortuous.
In order to derive approximate analytical solutions we will make the approximation that the variation of where C is an integration constant. Applying the boundary condition (6.8) gives and therefore Finally, imposing the boundary condition (6.6) on (6.11) gives a relation be-Journal of Applied Mathematics and Physics tween Â and B for prescribed values of n, φ , β and n K : The graph of Â against B is illustrated in Figure 7(a). We see that Â E < and that for ˆ0 A > , which from (3.42) describes fluid flux into the fracture at the The graph of B against Â is illustrated in Figure 7(b  In order to solve for the approximate analytical solution ( ) f u we use a non-orthogonal linesearch method described in [32]. We define the function (6.17), with the aid of (6.18), as

of Applied Mathematics and Physics
We then consider two coordinates, ( ) where k is a constant. Therefore in order to find the approximate analytical solution ( ) f u we proceed as follows: 1) Choose a value of m.
2) Let the first guess for 1 f be a vector given by the asymptotic solution (4.6).  Figure 8 shows that the graphs of the approximate analytical solution for the half-width approximately overlap the graphs of the numerical solution. We found that the approximate analytical solution is a better approximation for 3 n ≥ than for 2 3 n < ≤ . We also observe from Table 2 that with a specified value of B B = the parameter A obtained by the numerical method and the parameter Â obtained from Equation (6.15) are very closely matched with an error below 6% for 3 n ≥ and an increased error below 13% for 2 3 n < < . This is in agreement with Figure 6 in which the approximation that the curves are linear is clearly better for the less tortuous fracture with 4 n = and 3 n = than for the more tortuous fracture with 2.5 n =

3) Iterate with each solution for
. Results of the approximate analytical solution and the numerical solution are in good agreement as observed both in Figure 8 and in Table 2. In the absence of an exact analytical solution, it is clear that the approximate analytical solution is a reasonable alternative.

Conclusions
There are several novel features in this investigation. In the mathematical modelling of the tortuous hydraulic fracture, the hyperbolic crack law was used with the PKN approximation. Previous research [2] with the hyperbolic crack law used the elasticity model with the Cautchy principal value integral (2.5) for the Journal of Applied Mathematics and Physics stress in a two-sided fracture which has the advantage that the stress intensity  factor at the fracture tip is non-zero but is difficult to analyze mathematically.
With the PKN approximation, a nonlinear diffusion equation for a one-sided fracture was derived with a diffusion coefficient (2.24) not considered before in wards the solution of the constant pressure working condition [14]. A brief comparison of the results for the linear crack law and the hyperbolic crack law suggests that for very tortuous fractures, the constant pressure working condition is a dominant condition at the fracture entry.
The effect of the pressure due to contact regions is described by the pressure ratio φ . The length of the fracture after time t increases as φ increases for all values n which shows that the effect of touching asperities is to grow the length of the fracture. The effect of φ on the length of the fracture is more prominent in a very tortuous fracture with 2 3 n < < than in a less tortuous fracture with 3 n ≥ . The width averaged fluid velocity decreases as φ increases near the fracture entry but is almost independent of φ as the half-width decreases near the fracture tip. This shows that the effect of touching asperities is to decrease the width averaged fluid velocity in regions where the half-width is not small. As n decreases from 5 n = to 2 n = , the tortuosity of the fracture increases.
The value 3 n = describes the Reynolds flow law and even then there is some surface roughness. We saw that as n decreases from 4 n = to 2.5 n = the length of the model fracture after time t increases due to an increase in tortuosity. The tortuosity and the pressure due to contact regions both have the effect of increasing the fracture length.
For the hyperbolic crack law the spatial gradient of the half-width is singular at the fracture tip for 3 5 n < < and the thin fluid film approximation breaks down at the tip. However, for the Reynolds flow law 3 n = the spatial gradient of the half-width is finite while for 2 3 n < < it is zero and the thin fluid film approximation remains valid over the whole length of the fracture. As with the linear crack law, increasing tortuosity removes the singularity at the fracture tip.
The fluid flux at the fracture entry cannot be prescribed arbitrarily. It is determined by the invariant solution. The dependence on time must be of the form (3.42) and the constants A and B are determined numerically by the shooting method.
The approximate analytical solution is useful because there is no exact analytical solution for a hydraulic fracture with contact regions described by the hyperbolic crack law. There are no special cases which could lead to an analytical solution because constant pressure at the fracture entry is the only working condition. The approximate analytical solution is implicit in form and has to be solved iteratively and numerically in the final step. It is in good agreement with the numerical solution although the error increases slightly as the tortuosity increases. For the numerical solution neither of the constants, A or B, in the fluid Journal of Applied Mathematics and Physics flux at the fracture entry can be specified. For the approximate analytical solution either Â can be specified and B is obtained from (6.16) or B can be specified and Â is obtained from (6.15). This is because, built into the approximate analytical solution is the approximate result that the ratio of the width averaged fluid velocity to the speed of propagation of the fracture changes linearly along the fracture.