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 . 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 . Hydraulic fracturing is also used to enhance the extraction of oil and gas in large underground shale deposits .
Spence and Sharp  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   .
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  . 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 . 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 . Fareo and Mason also considered a fracture driven by a power law fluid in an impermeable rock mass . Anthonyrajah, Mason and Fareo  compared laminar and turbulent fluid driven fractures using the wall shear stress model of Emerman, Turcotte and Spence  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 with a modified crack law, which models the effect of the contact regions on the stress at the fluid-rock interface .
In this paper, we model the contact regions by the hyperbolic crack law which was introduced by Goodman . This is motivated by the discussion by Fitt et al.  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  , the contact regions in these papers were modelled by the linear crack law which was proposed by Pine and Cundall  and further discussed by Fitt et al. .
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.
2. Model formulation
We will first review the fluid flow results that were derived in . We will then close the fluid flow model by using both the hyperbolic crack law, which describes the presence of contact regions in the fracture, and the Perkins-Kern  and Nordgren  approximation, which relates the normal stress at the fracture walls with the half-width of the fracture. Lastly, we will express the governing equations in dimensionless form.
A tortuous hydraulic fracture with asperities at the fluid-rock interface is modelled. The tortuous fracture is replaced by a symmetric two-dimensional hydraulic fracture without asperities, as shown in Figure 1, but with a modified Reynolds flow law and a modified crack law. The x-axis is along the length of the model symmetric fracture, the z-axis is along the width and the y-axis is along the breadth. All quantities are independent of y. The half-width of the fracture is
Figure 1. Hydraulic fracture without asperities on the crack walls or contact regions.
. 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.
2.1. 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, , satisfies the Reynolds flow law :
and the width averaged fluid velocity is
The half-width of the fracture satisfies
which when expanded with the aid of (2.2) gives the nonlinear diffusion equation relating half-width of the fracture and the fluid pressure :
For hydraulic fractures with high tortuosity the Reynolds flow law has been shown to give theoretical results which diverge from experimental results    . It follows that the Reynolds flow law is not satisfactory in describing very tortuous fractures. Several authors have considered various flow models with the aim of accounting for high tortuosity in the fracture. We will consider a law that governs one of these models, the modified Reynolds flow law, in which in given in (2.2) is replaced by where n is a dimensionless constant and has dimensions . Both n and are obtained through experiments . Equations (2.2), (2.3) and (2.5) become:
2.2. Hyperbolic Crack Law
In  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  and further considered by Bandis et al. , Murphy et al.  and Fitt et al. .
A tortuous fracture can either be open without contact regions or partially open with contact regions. In an open fracture, the fluid pressure is sufficient to support the compressive normal stress at the fracture walls, , and therefore
For a partially open fracture, the fluid pressure is insufficient to support the compressive normal stress at the fracture walls and therefore
From (2.9) and (2.10), it is clear that the effective stress, , vanishes for an open fracture and is negative for a partially open fracture. While both an open and a partially open fracture were discussed in  and , in this paper we will discuss only a partially open fracture. We consider the hyperbolic crack law ,
where k is a negative constant, is the minimum half-width of the fracture and is the maximum half-width of the fracture. We make the approximation that which simplifies (2.11) to
Although in practice can never be zero even when the fracture is under the largest possible compressive stresses, is a reasonable approximation since  . It is clear that Equation (2.12) satisfies the negative effective stress condition, . Equation (2.12) can be rearranged to give
which clearly shows that the compressive normal stress at the fluid-rock interface, , is supported by both the fluid pressure and the pressure due to contact regions .
2.3. The PKN Approximation
Lastly, it is necessary to specify a relation between the compressive normal stress at the fluid-rock interface, , and the half-width of the fracture, . We will use the Perkins-Kern-Nordgren (PKN) approximation, which is widely used in the oil and gas industry to relate the normal stress at the crack walls to the half-width of the fracture,  :
where is the compressive stress at infinity (far-field compressive stress) within the rock and
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, , 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,
which describes a two sided fracture expanding from a point source, where the bar on the integral sign denotes the Cauchy principal value and G and are the shear modulus and Poisson ratio of the rock mass respectively. Equation (2.17) was first introduced by Spence and Sharp  and used by Fitt et al.  when considering a geothermal reservoir. While it is possible to derive a similarity solution for the integro-differential equation for the half-width of the fracture using the relation (2.17), the resulting boundary value problem is difficult to solve    . By expanding (2.17), Adachi and Peirce have shown that the PKN approximation is a good approximation away from an ε-neighbourhood of the fracture tip where .
2.4. Model Closure
From equations (2.13) and (2.15), the fluid pressure in the fracture is
Since 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 respectively become
and (2.8) becomes the nonlinear diffusion equation for the half-width of the fracture, :
When , the effective stress vanishes and the condition (2.9) for an open fracture is satisfied. Equations (2.20) to (2.22) reduce to the governing equations for an open fracture which was discussed in  and . In this paper we consider only partially open fractures ( ) for which the compressive normal stress at the fracture walls is supported by both the fluid pressure and the pressure due to contact regions given by (2.14).
Equation (2.22) is a nonlinear diffusion equation of the form
where the diffusion coefficient satisfies
and and are positive constants. Many problems in heat propagation, diffusion of fluids in porous media and ground water flow are of the form (2.23) but with diffusion coefficients different from (2.24) . To the best of our knowledge there is no application to hydraulic fracturing with a diffusion coefficient of the form (2.24).
2.5. Boundary and Initial Conditions
Consider first the boundary conditions at the fracture tip, , where is the length of the fracture. The half-width of the fracture vanishes at the 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 .
Consider next the condition at the fracture entry, . 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
where is the flux of fluid per unit breadth into the fracture at the entry . It is prescribed as far as the invariant solution will allow.
2.6. 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
We will work with instead of . The pressure is the difference between the pressure of the fluid in the fracture, , and the background pressure, , in the rock. Since is a constant, equations (2.3) to (2.8) on which the theory is based, remain unchanged when expressed in terms of instead of .
We choose characteristic quantities which are independent of n and so that the results for a range of values of n can be compared. We choose as the characteristic half-width. From (2.30) the characteristic value of is . The characteristic value of the x-coordinate is the initial length of the fracture 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, , 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 . The solution depends on the two parameters and and on , the fluid source at the fracture entry. The parameter contains all the dependence on n and . The physical significance of is
This completes the formulation of the mathematical model of a tortuous fracture with the hyperbolic crack law.
3. 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  and for hydraulic fractures , 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  . 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 equation (2.37) is
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 
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 Equation (2.37) with and is
and 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, , of the partial differential Equation (2.37) satisfies the condition
which expands to the first order partial differential equation
Consider the general case for which . By solving the differential equations of the characteristic curves of (3.10) the general solution is readily derived and since we obtain
where is an arbitrary function of the similarity variable given by
Unlike the group invariant solution for the linear fracture law, 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 equation for :
Since the constant does not appear in the ordinary differential Equation (3.13) we choose so that when . Thus
Consider next the boundary conditions. At the fracture tip the half-width of the fracture vanishes and (2.38) is satisfied. Expressed in terms of the group invariant solution (3.11), (2.38) becomes
Differentiating (3.15) with respect to “t” gives
We consider the general case for which is not a constant function. It follows that and therefore
which implies that
where is a constant. It follows that
which can be written as
where we have assumed that . Since the characteristic length of the fracture is the initial length of the fracture, and therefore
and the boundary condition at the fracture tip (3.15) becomes
Expressed in terms of the boundary condition (2.39) becomes
and the condition (2.41) becomes
Writing the source boundary condition (2.40) in terms of the group invariant solution (3.11) gives
We see that the time dependence of fluid flux into the fracture at the entry, , cannot be prescribed arbitrarily but must be of the form (3.28) determined by the invariant solution. The constants A and are prescribed if no restrictions are placed on them by the invariant numerical solution. It follows from (3.26) and (3.29) that
Since and are prescribed, (3.30) is a boundary condition on at . Expressed in terms of , the dimensionless fluid variables (2.42) to (2.45) become
We simplify the formulation by transforming from variables and to variables u and where
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 . 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 :
4. Asymptotic Solution
We now derive the asymptotic solution of the ordinary differential Equation (3.37) in an ε-neighbourhood of the fracture tip where . This asymptotic solution will be used in the numerical solution to offset the boundary condition at the fracture tip from to , in order to avoid singularities at the fracture tip ( )  .
We consider an asymptotic solution of the form
where p and C are constants to be determined with to satisfy the boundary condition and since for . While the resulting asymptotic solution is expected to satisfy the boundary condition at the fracture tip, (3.40), it is not expected to satisfy the boundary conditions at the fracture entry.
Substituting (4.1) into the ordinary differential Equation (3.37) and dividing the result by Cp gives:
The dominant terms as are the terms with the smallest values for the exponent of . To obtain an asymptotic solution these terms must balance and the remaining terms must vanish as . The dominant terms in (4.2) are the second and the fourth terms with exponents and , respectively. Equating these exponents gives
The special case of will be examined later. For we require . Equating coefficients of the second and fourth terms gives
and Equation (4.2) reduces to
for . The second term will vanish as since but the first term will vanish only if . It follows that the asymptotic solution for the ordinary differential Equation (3.37) is
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, , is identically satisfied.
Consider now the special case of . When , (4.2) becomes
The second term can no longer be balanced with the fourth term and it cannot be balanced with the first term for . We balance the first term with the fourth term to give . Equation (4.8) becomes
Since we are considering , condition (4.9) requires
Thus we have the result
which is the asymptotic solution for the linear crack law with .
Since the asymptotic solution for the hyperbolic crack law exists only for we conjecture that the solution of the boundary value problem, (3.37) to (3.41), exists only for . This conjecture will be investigated further when the numerical solution is derived in Section 5.
Consider now the behaviour of as . From (4.6) for ,
The asymptotic results, (4.6) and (4.13), will be used in the numerical solution in Section 5 to offset the boundary condition at to where .
Since the half-width of the fracture (3.44), with given by (4.6), is
as , provided . Differentiating (4.14) with respect to x and evaluating the limit of the result as gives
It is clear that the spatial gradient of the half-width of a two-dimensional hyperbolic model fracture is singular at for 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 . In the elasticity problem for an infinite body with a plane cut subject to symmetric loads Barenblatt describes similar behaviour near the cut tip .
5. 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, where , to the fracture entry, . 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 . 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 . A brief algorithm that shows how the problem is solved is provided below:
1) Specify parameters ,
2) Choose an initial guess for the shooting parameter B,
3) Solve the ODE (3.37) subject to the initial conditions and ,
Adjust B using the bisection method
Solve the ODE (37) subject to the initial conditions and ,
where is the asymptotic solution (4.6), is its derivative (4.12), , is the numerical solution obtained, is the moving target (3.38) due to the changing parameter B and the tolerance within which a solution is accepted is .
From steps 1 and 2 of the algorithm, it is clear that all parameters in the problem are specified. We fix for all values of n although depends on n. We consider analysis of partially open fractures, therefore we choose the parameter in the range . We expect the pressure difference at the fracture entry, , given by Equation (3.50) to be positive, therefore
where . A random initial guess of the shooting parameter is made.
In steps 3 and 4 of the algorithm the second order ODE is written as the system of two first order ordinary differential equations.
where . The boundary value problem is solved as an initial value problem by shooting a solution from the fracture tip to the fracture entry. The differential Equations (5.3) and (5.4) however have singularities at the fracture tip. We use (4.6) and (4.12) for and as . From (4.13) we see that (4.12) is singular at for . Thus (5.3) is singular at for . Also in (5.4) we find that
which is singular at for ,
which is singular at for and
which is singular at for . It follows that the system of equations (5.3) and (5.4) is singular for as . While there are no singularities in and for , the algorithm is unable to solve the problem when the initial condition is used in the entire range . We therefore offset the boundary conditions at the fracture tip from to for the entire range . Therefore the system of Equations (5.3) and (5.4) is solved by shooting a solution from , where , to the fracture entry. The initial conditions on and are evaluated at using (4.6) and (4.12):
The shooting method is applied iteratively, with the parameter B updated at each iteration by the bisection method, until the boundary condition at the fracture entry (3.38),
is satisfied. The boundary condition (3.41) is identically satisfied while the value of the constant A in the source flux of fluid at the entry (3.42) can now be determined from the boundary condition (3.39),
The numerical solution for , obtained from solving the boundary value problem (3.37) to (3.42), is substituted into equations (3.43) to (3.48) which give the properties of the fracture as well as the behaviour of the fluid inside the fracture.
We now analyze the half-width of the model fracture . Figure 2 shows the half-width of a partially open model fracture plotted against x for a range of dimensionless times and for , and . From (3.50) a tortuous fracture with the hyperbolic crack law admits only one working condition of constant pressure difference at the fracture entry. Since the far-field compressive stress is also constant this implies for the hyperbolic fracture the pressure of fluid injection at the fracture entry, , is constant. Due to the PKN approximation the half-width of the model fracture is also constant for at the fracture entry and is given by (3.51). This is clearly shown in Figure 2 where . Both the half-width for and the length of the model fracture grow as time increases resulting in the volume of the model fracture per unit breadth, , which from (2.45) is the area under the half-width curve, growing with increasing time t.
The behaviour of at the fracture tip given by (4.15) is clearly illustrated in Figure 2. The chosen values, , and , lie in the ranges , and and we see that is negative infinity, finite and zero at the tip.
The length, , of a partially open fracture is plotted against time for a range of values of in 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.
The numerical solution of the boundary value problem (3.37) to (3.42) was obtained for a range of values of , and n with the diffusion constant prescribed as . The values obtained by the shooting method for A and B in the source flux of fluid at the entry, (3.42), are presented in 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 entry, , is required to propagate a very tortuous fracture. From (3.43), the speed of propagation of the fracture is
Figure 2. Numerical solution for the half-width of a partially open model fracture, with , plotted against x with the diffusion constant and for increasing values of time t with (a) , (b) , (c) .
Figure 3. Numerical solution of the length of a partially open model fracture plotted against time t for , , , with , diffusion constant and (a) , (b) , (c) .
Figure 4. Numerical solution of the length of the fracture of a partially open model fracture plotted against time t for , , , with , diffusion constant and (a) , (b) , (c) .
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 is plotted against time t for a range of values of n, and .
We also see from Table 1 that for fixed values of n and , A decreases as decreases and therefore the flux of fluid at the entry is weaker when the width of the partially open fracture is smaller. For fixed values of n and , the constant B increases as decreases and therefore the speed of propagation of the fracture decreases for partially open fractures of smaller width, in agreement with 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 .
6. Width Averaged Fluid Velocity
The ratio of the width averaged fluid velocity, , given by (3.46) to the speed of propagation of the fracture tip,
Using the asymptotic solutions (4.6) for and (4.12) for it can be verified that
The width averaged fluid velocity therefore tends to the speed of propagation of the fracture tip as . Using the boundary conditions (3.38) and (3.39) for and at we find that
The curves of the width averaged fluid velocity ratio (6.2) plotted against u for therefore have end points at and at .
In Figure 6, the velocity ratio curves (6.2) for a partially open fracture with and , and respectively are plotted against u for , with the diffusion constant and for , and . The curves all pass through the point . In Figure 6(c) for example, the end points at for , and are ,
Table 1. Numerical values of the parameters A and B in the source fluid flux at the fracture entry for .
Figure 5. Numerical solution of the source fluid flux at the entry of a partially open model fracture plotted against time t for , diffusion constant and (a) , (b) , (c) .
Figure 6. Velocity ratio curves of a partially open model fracture, , plotted against u with the diffusion constant and for parameter values , and respectively with (a) , (b) , (c) .
and respectively and are in good agreement with the exact result which, with the aid of Table 1, gives , and . We see from 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 increases approximately
linearly with u along the fracture for the values of the parameters considered. The linear approximation is best for 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 along the fracture is linear in u for
  . In the approximate solution, we will replace , A and B by , and . By using the equation of the straight line between the points and and using Equation (6.2) we obtain
for . Equation (6.5) is a first order ordinary differential equation for subject to the boundary conditions
By using the differential Equation (6.5) and boundary condition (6.6) we see that the boundary condition (6.7) is identically satisfied. Integrating the ordinary differential Equation (6.5) gives
where C is an integration constant. Applying the boundary condition (6.8) gives
Finally, imposing the boundary condition (6.6) on (6.11) gives a relation between and for prescribed values of n, , and :
Equation (6.12) can be re-expressed as
If is prescribed then
The graph of against is illustrated in Figure 7(a). We see that and that for , which from (3.42) describes fluid flux into the fracture at the entry ( ), then it is necessary that . If is prescribed then
The graph of against is illustrated in Figure 7(b). The inequalities and are again apparent.
where and . Equation (6.17) is implicit because appears on both sides. Therefore in order to find the approximate analytical solution , a numerical method is required in which (6.17) is iterated until convergence is reached. In (6.17) the numerator has the factor and the denominator has the factor . However, the iterative method used to solve for requires that the iterated function (6.17) is neither factorised nor a power law. The parameters n, , and are prescribed. Unlike the numerical solution, either or can be prescribed. If is prescribed then is obtained from (6.16) while if is prescribed then is obtained from (6.15). In this work we prescribe and therefore the expression appearing in (6.17) is
For , , and we found that if we let
Figure 7. Graphs of parameters and of the approximate analytical solution with , , and : (a) Graph of against , (b) Graph of against .
then the approximate analytical value is in reasonable agreement with the numerical value since the fracture is very tortuous with .
In order to solve for the approximate analytical solution we use a non-orthogonal linesearch method described in . We define the function (6.17), with the aid of (6.18), as
We then consider two coordinates, and , with which to calculate a gradient m:
We re-express (6.20) to give :
where k is a constant. Therefore in order to find the approximate analytical solution we proceed as follows:
1) Choose a value of m.
2) Let the first guess for be a vector given by the asymptotic solution (4.6).
3) Iterate with each solution for used as in the next iteration,
4) The iteration is continued until convergence is reached at the kth step where the solution .
We choose which in a few iterations gives an approximate analytical solution . The gradient m is required to be negative because the spatial gradient of the upper half-width of the fracture is negative.
For a very tortuous partially open fracture with we found that more iterations are required to obtain the approximate analytical solution for the half-width than for a less tortuous fracture with . The computational time required to find the approximate analytical solution for the half-width of a partially open fracture with , and is approximately 0.7 seconds. 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 than for . We also observe from Table 2 that with a specified value of 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 and an increased error below 13% for . 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 and than for the more tortuous fracture with . 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.
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  with the hyperbolic crack law used the elasticity model with the Cautchy principal value integral (2.5) for the stress in a two-sided fracture which has the advantage that the stress intensity
Figure 8. Numerical solution () and approximate analytical solution () for the half-width of a partially open model fracture, , plotted against x with
diffusion constant and for , with (a) , (b) , (c) .
Table 2. Numerical values A compared to approximate analytical values for with , and various values of and .
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 hydraulic fracturing. Although an exact analytical solution could not be derived, an existing model was developed based on the width-averaged fluid velocity which led to a new approximate analytical solution. In the mathematical analysis of the problem, the Lie point symmetry of the nonlinear diffusion equation had the special property that the coefficient vanished which had a significant effect on the solution. A shooting method was developed to solve the problem numerically and a recently established non-orthorgonal line search method was applied to obtain an iterative solution for the implicit approximate analytical solution. The following discussion gives a summary of the results obtained.
For the hyperbolic crack, law the only working condition at the fracture entry is constant pressure. From the PKN approximation, the half-width of the fracture remains constant at the fracture entry. Mathematically these results follow because in the Lie point symmetry, which generates the group invariant solution, the coefficient vanishes. The exponent n in the power law in the modified Reynolds flow law lies in the finite range . While the linear crack law yielded a full range of working conditions of both fluid injection and extraction at the fracture entry and , it was observed that for a very tortuous fracture, all solutions for a range of working conditions converged towards the solution of the constant pressure working condition . 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 than in a less tortuous fracture with . 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 to , the tortuosity of the fracture increases. The value describes the Reynolds flow law and even then there is some surface roughness. We saw that as n decreases from to 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 and the thin fluid film approximation breaks down at the tip. However, for the Reynolds flow law the spatial gradient of the half-width is finite while for 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 flux at the fracture entry can be specified. For the approximate analytical solution either can be specified and is obtained from (6.16) or 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.
D.P.M. acknowledges financial support from the National Research Foundation, South Africa, under Grant No 132189.
M.R.R.K. acknowledges financial support from the National Research Foundation, Pretoria, South Africa.
 Fitt, A.D., Mason, D.P. and Moss, E.A. (2007) Group Invariant Solution for a Pre-Existing Fluid-Driven Fracture in Impermeable Rock. Zeitschrift fur Angewandte Mathematik und Physik, 58, 1049-1067.
 Fitt, A.D., Kelly, A.D. and Please, C.P. (1995) Crack Propagation Models for Rock Fracture in a Geothermal Energy Reservoir. SIAM Journal on Applied Mathematics, 55, 1592-1608.
 Emerman, S., Turcotte, D.L. and Spence, D.A. (1986) Transport of Magma and Hydrothermal Solutions by Laminar and Turbulent Fluid Fracture. Physics of the Earth and Planetary Interiors, 36, 276-284.
 Spence, D.A., Sharp, P. and Turcotte, D.L. (1987) Buoyancy-Driven Crack Propagation: A Mechanism for Magma Migration. Journal of Fluid Mechanics, 174, 135-153.
 Fareo, A.G. and Mason, D.P. (2011) Group Invariant Solutions for a Pre-Existing Fluid-Driven Fracture in Permeable Rock. Nonlinear Analysis: Real World Applications, 12, 767-779.
 Fareo, A.G. and Mason, D.P. (2013) Group Invariant Solution for a Pre-Existing Fracture Driven by a Power Law Fluid in Impermeable Rock. Communications in Nonlinear Science and Numerical Simulation, 18, 3298-3316.
 Anthonyrajah, M., Mason, D.P. and Fareo, A.G. (2013) Propagation of a Turbulent Fluid Fracture. International Journal of Non-Linear Mechanics, 54, 105-114.
 Kgatle, M.R.R. and Mason, D.P. (2014) Propagation of a Linear Two-Dimensional Hydraulic Fracture with Tortuosity. International Journal of Non-Linear Mechanics, 63, 39-53.
 Kgatle, M.R.R. and Mason, D.P. (2016) Linear Hydraulic Fracture with Tortuosity: Conservation Laws and Fluid Extraction. International Journal of Non-Linear Mechanics, 79, 11-25.
 Pine, R.J. and Cundall, P.A. (1985) Applications of the Fluid-Rock Interaction Programme (FRIP) to the Modelling of Hot Dry Geothermal Energy Systems. Proceedings of the International Symposium on Fundamentals of Rock Joints, Björkliden, 15-20 September 1985, 293-302.
 Qian, J., Chen, Z., Zhan, H. and Guan, H. (2011) Experimental Study of the Effect of Roughness and Reynolds Number on Fluid Flow in Rough-Walled Single Fractures: A Check of Local Cubic Law. Hydrological Processes, 25, 614-622.
 Witherspoon, P.A., Wang, J.S.Y., Iwai, K. and Gale, J.E. (1980) Validity of Cubic Law for Fluid Flow in a Deformable Rock Fracture. Water Resources Research, 16, 1016-1024.
 Bandis, S., Lumsden, A.C. and Barton, N.R. (1983) Fundamentals of Rock Joint Deformation. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 20, 249-268.
 Murphy, H., Dash, Z., Zyvoloski, G., White, A. and Wing, M. (1983) Fundamentals of Rock Joint Deformation. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 20, 249-268.
 Cook, N.G.W. (1992) Natural Joints in Rock: Mechanical, Hydraulic and Seismic Behaviour and Properties under Normal Stress. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 29, 198-223.
 Adachi, J.I. and Peirce, A.P. (2008) Asymptotic Analysis of an Elasticity Equation for a Finger-Like Hydraulic Fracture. Journal of Elasticity, 90, 43-69.
 King, J.R. and Bowen, M. (2001) Moving Boundary Problems and Non-Uniqueness for Thin Film Equations. European Journal of Applied Mathematics, 12, 321-356.
 Huppert, H.E. (1982) The Propagation of Two-Dimensional and Axisymmetric Viscous Gravity Currents over a Rigid Horizontal Surface. Journal of Fluid Mechanics, 121, 43-58.
 Ibragimov, N.H. and Anderson, R.L. (1994) One-Parameter Transformation Groups. In: Ibragimov, N.H., Ed., CRC Handbook of Lie Group Analysis of Differential Equations, Volume 1: Symmetries, Exact Solutions and Conservation Laws, CRC Press, Boca Raton, 7-14.
 Pallares, M.M.R. and Rodriguez, C.W. (2014) Methodology for Solving the Divergence of Fixed Point Method for the Solution of a Nonlinear Equation. ARPN Journal of Engineering and Applied Sciences, 9, 2177-2182.