2.1 Percolation theory and fracture activation under constant fluid pressure

In percolation theory, one starts with connected closed bonds and opens them with a certain probability. Percolation theory states that one infinite percolating cluster is formed when a threshold fraction, $F_{c}$ , of the bonds is opened (Stauffer Reference Stauffer1979). The value of $F_{c}$ depends on the topology of the formed network. For instance, $F_{c}=1/2$ for a two-dimensional square lattice of bonds and $F_{c}=1$ for a one-dimensional network. Moreover, several bulk properties of the percolating network scale with the difference between the fraction of opened bonds, $F$ , and $F_{c}$ when $F-F_{c}\ll F_{c}$ . That is $X\sim |F-F_{c}|^{a}$ where $X$ is a bulk property of interest and $a$ is a universal scaling exponent that only depends on the dimension of the network and not on its topological details.

Among several bulk properties that follow universal scalings are the percolating network’s correlation length, its porosity and its hydraulic conductivity. The correlation length of a network, $\unicode[STIX]{x1D709}$ , is defined as the average distance between two points in the network and it is a measure of how sparse the network is. Above the percolation threshold, the percolating cluster’s correlation length is the average radius of gyration of the clusters of closed bonds. Portions of the percolating network smaller than the correlation length look fractal, while the network looks homogeneous when the length scale is larger than the correlation length. As the fraction of opened bonds increases beyond $F_{c}$ , the correlation length decays as $\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D709}_{0}(F-F_{c})^{-\unicode[STIX]{x1D708}}$ where $\unicode[STIX]{x1D709}_{0}$ is the correlation length of the existing bonds in the system. For example, if the bonds form a square lattice, $\unicode[STIX]{x1D709}_{0}$ is equal to the lattice spacing. In two dimensions, $\unicode[STIX]{x1D708}=4/3$ while it is estimated to be about $0.88$ in three dimensions (Sahimi Reference Sahimi2011).

Similarly, the fraction of the bonds that are open and belong to the percolating cluster, $S$ , and the percolating network’s hydraulic conductivity, $K$ , scale as $S\sim \unicode[STIX]{x1D709}^{-\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D708}}$ and $K\sim k_{0}\unicode[STIX]{x1D709}^{-\unicode[STIX]{x1D716}/\unicode[STIX]{x1D708}}$ , respectively. $k_{0}$ is the permeability of the network when all the bonds are opened. For $D=2$ , $\unicode[STIX]{x1D6FD}=5/36$ and $\unicode[STIX]{x1D716}\approx 1.3$ , while, for $D=3$ , $\unicode[STIX]{x1D6FD}\approx 0.42$ and $\unicode[STIX]{x1D716}\approx 2.0$ (Stauffer & Aharony Reference Stauffer and Aharony1994; Kozlov & Lagues Reference Kozlov and Lagues2010; Wang et al. Reference Wang, Zhou, Zhang, Garoni and Deng2013). These scalings are valid when $\unicode[STIX]{x1D709}_{0}\ll \unicode[STIX]{x1D709}\leqslant B$ where $B$ is the system size. If $\unicode[STIX]{x1D709}\gg B\gg \unicode[STIX]{x1D709}_{0}$ , the bulk properties scale with $B$ , e.g.  $S\sim B^{-\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D708}}$ . For more details about percolation theory and its application to fluid transport, the reader is referred to Stauffer & Aharony (Reference Stauffer and Aharony1994), Sahimi (Reference Sahimi2011) and Hunt, Ewing & Ghanbarian (Reference Hunt, Ewing and Ghanbarian2014).

Returning to the fracturing process, following the slippage mechanism, a weak plane will undergo slippage when the stress tangent to the plane is larger than the product of a friction coefficient and an effective normal stress. If the stress field is constant, each discontinuity is associated with a critical fluid pressure, $P_{c}$ , the fluid has to reach in order for slippage to occur:

(2.1) $$\begin{eqnarray}P_{c}=\unicode[STIX]{x1D70E}_{n}-\frac{\unicode[STIX]{x1D70E}_{t}}{f_{r}},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{t}$ , $\unicode[STIX]{x1D70E}_{n}$ are the tangential and normal stresses acting on the fracture, respectively, and $f_{r}$ is the friction coefficient. The initial aperture of the pre-existing discontinuity is very small such that fluid flow is negligible. Upon slippage, its aperture becomes large enough to allow fluid flow. The creation of fluid filling volume is due to the mismatch between the asperities of the two surfaces induced by slippage. To simplify the model, the effects of stress perturbation due to surface displacement are ignored. Hence, the time variation in the aperture within the growing fracture is ignored and a characteristic hydraulic aperture is assumed. When the viscosity of the fluid saturating the weak planes is much smaller than that of the injected fluid, the fracturing fluid’s pressure controls the activation process as the fluid perturbs the pressure within the natural fracture.

If the viscous pressure drop is negligible and the fluid is injected at a constant pressure, the fraction of pre-existing weak planes whose critical pressure is below the fluid’s pressure can open. When the fractures’ critical pressures follow a homogeneous distribution, the fraction of fractures that have a critical pressure below the fluid pressure, denoted as $F$ , is given by:

(2.2) $$\begin{eqnarray}F(P)=\int _{P_{c_{min}}}^{P}f_{pc}(x)\,\text{d}x,\end{eqnarray}$$

where $P$ is the fluid pressure and $f_{pc}(x)$ is the probability distribution function of the fractures’ critical pressures. $f_{pc}$ becomes random when the critical pressures are uncorrelated with the fracture’s orientation. This can be the case when the stress field is heterogeneous over large length scales. $P_{c_{min}}$ is the minimum critical pressure in the distribution. As the injected fluid propagates through the network of activated fractures, some of the fractures that have critical pressures below the fluid pressure will remain unactivated. These fractures are connected to the activated fracture network via weak planes with critical pressures that are larger than the fluid pressure. Therefore, they will not be activated and the injected fluid will be constrained within the formed percolating network.

The conformity of the percolation problem and the mechanism of activating pre-existing weak planes under a constant pressure can be stated as follows. The fraction of opened bonds in percolation theory, $F$ , corresponds to the fraction of natural fractures whose critical pressure is lower than $P$ . The fraction of the opened bonds that belong to the percolating network, $S$ , represents the fraction of activated fractures. The threshold fluid pressure, $P_{cc}$ , needed to form a percolating network of activated fractures can be obtained by taking the inverse of (2.2) evaluated at $F=F_{c}$ . Since the value of $F_{c}$ depends on the topology of the network and its dimension, one expects the same dependence for the value of $P_{cc}$ . Additionally, the value of $P_{cc}$ depends on the functional form of $f_{pc}$ .

To relate the bulk properties of the activated network of fractures to the fluid pressure, equation (2.2) can be linearized so that:

(2.3) $$\begin{eqnarray}F-F_{c}=k_{f}(P-P_{cc}),\end{eqnarray}$$

where $k_{f}$ is a proportionality constant whose order of magnitude is equal to the reciprocal of the probability distribution’s standard deviation, $\unicode[STIX]{x1D6FF}_{pc}$ . For a normal distribution, for instance, $k_{f}=1/(\sqrt{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FF}_{pc})$ . This linearization applies when $F-F_{c}\ll F_{c}$ ; the regime where the universal scalings apply. Given (2.3), the dependence of the network properties on the fluid pressure can be written as:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{0}k_{\unicode[STIX]{x1D709}}[k_{f}(P-P_{cc})]^{-\unicode[STIX]{x1D708}}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle S=k_{s}[k_{f}(P-P_{cc})]^{\unicode[STIX]{x1D6FD}}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle K=k_{0}k_{k}[k_{f}(P-P_{cc})]^{\unicode[STIX]{x1D716}}, & \displaystyle\end{eqnarray}$$

where $k_{\unicode[STIX]{x1D709}}$ , $k_{s}$ and $k_{k}$ are dimensionless proportionality constants and depend on the geometry of the network and its dimension while the critical exponents are universal and only depend on the dimension of the network. $\unicode[STIX]{x1D709}_{0}$ is the correlation length of the existing natural fractures and its equivalent length scale in regular percolation is the lattice spacing.

2.2 Network formation under controlled injection rate

Previously, it has been shown that when the fluid pressure is constant, the dependence of the activation of natural fractures on the fluid pressure leads to the formation of a network whose bulk properties follow the same universality scalings found in regular percolation. Now, we consider the case where the injection rate is controlled and show that one can use the derived scalings, equations (2.4)–(2.6), locally to model fluid transport within the fractured rock.

There are four length scales that control the activation process: (i) the average length of the pre-existing weak planes, $L_{av}$ ; (ii) the correlation length of pre-existing natural fractures, $\unicode[STIX]{x1D709}_{0}$ , that is related to the second moment of the fractures’ length distribution (Balberg et al. Reference Balberg, Anderson, Alexander and Wagner1984); (iii) a characteristic correlation length of the cluster of activated natural fractures $\unicode[STIX]{x1D709}_{ch}$ ; and (iv) the cluster’s radius, $R$ . $\unicode[STIX]{x1D709}_{0}$ is the average distance between these fractures and is of the order of the fracture spacing. If the density of these fractures is large and their lengths are of the same order of magnitude, the average spacing and hence $\unicode[STIX]{x1D709}_{0}$ become of the same order as $L_{av}$ . In this case, the network of pre-existing natural fractures is homogeneous, that is its fractal dimension is equal to the Euclidean dimension of the network. However, if the density approaches a threshold value, $\unicode[STIX]{x1D709}_{0}\rightarrow \infty$ producing a fractal network of pre-existing natural fractures. In this paper, we are limiting ourselves to the case in which $\unicode[STIX]{x1D709}_{0}$ is finite and very close to $L_{av}$ . $\unicode[STIX]{x1D709}_{ch}$ is defined as the average radius of gyration of the clusters of unactivated fractures within the network of activated fractures. It is a function of the ratio of the pressure drop across the network of activated fractures to $\unicode[STIX]{x1D6FF}_{pc}$ , i.e.  $\unicode[STIX]{x1D709}_{ch}\sim (P_{ch}/\unicode[STIX]{x1D6FF}_{pc})^{-\unicode[STIX]{x1D708}}$ . $P_{ch}$  is the characteristic pressure drop across the cluster of activated fractures which will be derived later on. If the viscous pressure drop is negligible, $\unicode[STIX]{x1D709}_{ch}\rightarrow \infty$ and when an extremely viscous fluid is injected at a high rate, $\unicode[STIX]{x1D709}_{ch}$ is expected to approach $\unicode[STIX]{x1D709}_{0}$ .

Figure 1. Type of networks produced based on the importance of the viscous pressure drop. In all cases, the pre-existing natural fractures form a homogeneous network, i.e.  $\unicode[STIX]{x1D709}_{0}\sim O(L_{av})$ . (a) Shows a fractal network produced when the viscous pressure drop is not important and the fracturing process is controlled by the distribution of the critical pressures. (b) Shows a homogeneous network produced when the viscous pressure drop overcomes the variability of the critical pressures. (c) Shows a network that looks the same as a network near the percolation threshold at the length scale $\unicode[STIX]{x1D709}_{ch}$ . Over the length scale of the radius, the network looks heterogeneous.

Depending on the competition between the viscous pressure drop required to drive the flow and the pressure change induced by activating a pre-existing fracture, the growth dynamics can behave differently. When the viscous pressure drop is negligible when compared to the variability of the critical pressures, a scale-invariant fractal network is produced while a homogeneous network at all length scales is produced when the viscous pressure drop is dominant in the fracturing process. Between these two extreme cases, an intermediate regime can be identified where the viscous pressure drop is important over the cluster radius but can be neglected over the length scale, $\unicode[STIX]{x1D709}_{ch}$ . In this case, the produced network has a structure at small length scales similar to a percolating network formed near the percolation threshold and looks heterogeneous over the cluster radius. In all three cases, the injected fluid pressure must be at least equal to the threshold critical pressure, $P_{cc}$ .

In the first extreme case where the viscous pressure drop is negligible, i.e.  $\unicode[STIX]{x1D709}_{0}\sim O(L_{av})\ll R\ll \unicode[STIX]{x1D709}_{ch}$ , the network of activated fractures is fractal. Its structure is the same as a percolating network formed at the threshold value as shown in figure 1(a). In this case, the injected fluid activates the accessible fractures with the smallest critical pressure and its growth behaviour can be described by the invasion percolation without trapping algorithm (Wilkinson & Willemsen Reference Wilkinson and Willemsen1983). In this algorithm, each bond in a network is associated with a random resistance to opening. At each simulation step, a bond with the lowest resistance that is connected to the network of opened bonds is opened. Bonds with a low resistance connected to the network via a path of bonds with higher resistances remain closed. This algorithm produces a fractal network with the same fractal dimension, $1.9$ when $D=2$ , as a network generated at the percolation threshold following regular percolation rules.

A recent study by Tayeb et al. (Reference Tayeb, Sahimi, Aminzadeh and Sammis2013) has argued that a hydraulically fractured network of natural fractures, in the Geysers geothermal field in northeast California, has the same fractal dimension as a percolation cluster formed at the threshold value. The results were interpreted by a model developed by Sahimi, Robertson & Sammis (Reference Sahimi, Robertson and Sammis1993) where large scale heterogeneities in the resistances to activate a fracture lead to a random fracturing process similar to the process used to form a percolating network. As mentioned in § 2.2, heterogeneity in the stress field, at length scales larger than $L_{av}$ , can lead to a random distribution of critical pressures. In order for the randomness in the resistances to control the cluster growth and produce a fractal network, our model predicts that the variability of the critical pressures must be much larger than the viscous pressure drop required to drive the flow, i.e.  $\unicode[STIX]{x1D709}_{0}\approx L_{av}\ll R\ll \unicode[STIX]{x1D709}_{ch}$ .

The other extreme case to consider is when the viscous pressure drop within the fractures is dominant. If $\unicode[STIX]{x1D709}_{ch}\sim \unicode[STIX]{x1D709}_{0}\ll R$ , one expects that the injected fluid will easily overcome the largest resistance and hence will activate all natural fractures at the fluid front. The growth behaviour in this case can be described by a linear diffusion equation and the network produced via this behaviour looks homogeneous at all length scales as shown in figure 1(b). The fractal dimension of the produced network is the same as the Euclidean dimension of space embedding the network. This case is not relevant to fractured geological reservoirs as Sahimi (Reference Sahimi2011) has shown that all networks of activated fractures either form a scale-invariant fractal network or a network that is fractal at small length scales and homogeneous at large length scales.

Between such extreme cases, one can identify a regime where the viscous pressure drop across the cluster of activated fractures is important but it can be neglected over distances of order $\unicode[STIX]{x1D709}_{ch}$ , where

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D709}_{0}\ll \unicode[STIX]{x1D709}_{ch}\ll R.\end{eqnarray}$$

The produced network in this regime is expected to be near the percolation threshold throughout the cluster. Below $\unicode[STIX]{x1D709}_{ch}$ , it looks fractal while the network’s properties are heterogeneous at large length scales as shown in figure 1(c). In this paper, a model is developed to describe the growth of a network of activated fractures when (2.7) applies. If the viscous pressure drop is negligible over $\unicode[STIX]{x1D709}_{ch}$ and $\unicode[STIX]{x1D709}_{ch}\ll R$ , a local fluid pressure can be defined over that region. In § 2.1, it has been shown that when the many connected fractures feel a constant fluid pressure, their activation will lead to the formation of a percolating network. Furthermore, the properties of the activated network will scale with the fluid pressure in that region. Therefore, one can see that the local properties of the network can be described using (2.4)–(2.6) if the local correlation length is much larger than $\unicode[STIX]{x1D709}_{0}$ but is much smaller than the cluster radius. Since the viscous pressure drop required to drive the flow of the fluid is important over the cluster radius, the local properties of the network such as the fraction of activated fractures, $S$ , and the permeability of the network, $K$ , will evolve spatially and temporally as fluid is injected. Hence, one can couple fluid transport with the activation of natural fractures through the dependence of $S$ and $K$ on the fluid pressure. Since the fluid pressure has to be at least equal to $P_{cc}$ to form a percolating network, it is expected that the pressure at the edge of the cluster will be very close to $P_{cc}$ and the network will look fractal in that region.

Figure 2. A sketch of the cluster of activated fractures. The solid lines represent the network of activated fractures while the dashed bonds represent the weak planes with critical pressures below the local pressure that nonetheless remain unactivated. On the right-hand side, the fluid pressure is slightly larger than the threshold critical pressure to ensure the formation of a percolating network. As the fluid pressure increases, the density of the percolating network increases. The network on the left represents a denser percolating network than the one on the right because the local pressure is higher due to the viscous pressure drop required to drive the flow.

Figure 2 is an illustration of how the activated fracture network will look when (2.7) applies. The middle cartoon in the figure shows the entire cluster of activated fractures if the pre-existing natural fractures form a square lattice. The right-hand sketch shows how a network of activated fractures is expected to form when the local pressure is slightly larger than the threshold critical pressure, $P_{cc}$ . The solid lines represent the activated fracture network and the dashed lines are disconnected fractures which remain unactivated despite having critical pressures below the local fluid pressure. Since the pressure drop required to drive the flow is important over the cluster’s radius, the pressure in the interior region of the network is higher than $P_{cc}$ and, hence, the network there is more connected as shown in the left hand sketch. Based on this picture of how a network of activated fractures propagates, we shall introduce a continuum model that couples fluid flow with permeability and porosity evolution within a porous media.

2.3 Transport equation

The fractured rock is modelled as a dual porosity medium: a primary porosity due to interconnected activated fractures and a secondary porosity due to the pores of the rock matrix. Fluid flow in the secondary porosity is coupled with a pressure-driven flow within the primary porosity through a process called leak off. Each natural fracture within the primary porosity is assumed to be activated when the fluid pressure reaches a critical value. The activation process of a network of fractures with random resistances can be described by percolation theory. Hence, fluid flow within the primary porosity is controlled by the interconnectivity of the activated fractures but not by the detailed geometry of the formed network. Equation (2.7) applies and thus the primary porosity is allowed to change with time and position as weak planes are fractured. The secondary porosity is assumed to be constant and unaffected by the pressure of the leaked fluid. The time scale at which the primary porosity evolves is much larger than the time scale to completely activate a pre-existing weak plane.

To relate fluid flow with the growth of a cluster of fractures, a local mass balance that accounts for changes in the primary porosity upon the injection of the fracturing fluid can be used:

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}+V_{l}=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D700}$ is the primary porosity if all the local pre-existing fractures are activated, which is assumed to be constant when considering a homogeneous number density of the natural weak planes. $V_{l}$ accounts for the total leakage rate from the local network per unit medium volume. $\boldsymbol{q}$ is the superficial flux of the fluid through the network and it can be described using Darcy’s law where the permeability is given by (2.6). $\unicode[STIX]{x1D700}S$ is the primary porosity and its dependence on the fluid pressure is given by (2.5).

Two boundary conditions can be used to fully define the problem. One can specify the flux at the injection point and the fluid pressure at the edge of the cluster. Since a percolating network forms when the fluid pressure is at least equal to $P_{cc}$ , the fluid pressure at the edge of the cluster can be set to be equal to $P_{cc}$ . To find how the cluster radius grows with time, one can perform a volume integral on (2.8) and apply the flux at the injection point and zero flux past the moving front of the cluster, $R(t)$ , to get:

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\left(\int S\,\text{d}V\right)=q_{inj}-\int V_{l}\,\text{d}V,\end{eqnarray}$$

where $q_{inj}$ is the flux at the injection point.

The above system of equations governs the growth behaviour of a $D$ -dimensional network of activated fractures by coupling fluid transport and changes in the rock’s porosity. As can be seen from (2.8), the competition between the leakage term and the pressure-driven flow controls how the primary porosity evolves. To demonstrate how this competition can lead to two distinct self-similar growth behaviours, we will first analyse the growth of a one-dimensional network. For this simple case, an analytical solution for the length of the network can be obtained. Furthermore, the analysis for this network will prove useful in justifying some of the assumptions that will be used to solve for the two-dimensional case. The solution for a two-dimensional network will be presented in § 4.

3.1 Governing equations

Assuming that the one-dimensional network forms a rectangular fracture, equation (2.8) can be reduced to:

(3.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}x}q+2v_{l}=0.\end{eqnarray}$$

In this equation, $q$ is the fluid flux per width and the fluid flows in the $x$ direction whereas the leaking fluid flows perpendicular to the two fracture surfaces into the rock matrix. The leakage term in (2.8) becomes the local leakage velocity, $v_{l}$ , that is defined as the volumetric flow rate per surface area. The porosity term, $\unicode[STIX]{x2202}S/\unicode[STIX]{x2202}t$ , in (2.8) is set to be zero since the local porosity consists of the volume of the locally activated discontinuity. Therefore, it will not change with time when its aperture is assumed to be constant. This assumption holds when the pressure within the activated discontinuities does not overcome the far-field normal stress acting on asperities on the fracture surfaces that prop the fracture open and the effects of shear dilation following the initial activation of a discontinuity are ignored. Shear dilation results from stress perturbations due to slippage and is a function of the magnitude of displacement and the surface roughness (Willis-Richards, Watanabe & Takahashi Reference Willis-Richards, Watanabe and Takahashi1996; Olsson & Barton Reference Olsson and Barton2001). Since the network’s conductivity was described by an effective permeability when deriving (2.8), a hydraulic aperture, $b$ , will be used for the one-dimensional network. Since the competition between pressure-driven flow and leakage depends on a characteristic conductivity parameter, changes in the fracture’s aperture are not expected to change the qualitative growth behaviour of the fracture.

The pressure at the fluid front $L(t)$ must be equal to $P_{c}$ , in order for slippage to occur and maintain the continuous growth of the activated fracture. Hence, one can use the following boundary condition:

(3.2) $$\begin{eqnarray}P=P_{c},\quad \text{at }x=L(t).\end{eqnarray}$$

The other boundary condition required to solve (3.1) can be written as:

(3.3) $$\begin{eqnarray}q=\frac{Q}{w},\quad \text{at }x=0,\end{eqnarray}$$

through which the injection rate at the injection well $Q$ can be specified.

Finally, equation (2.9) to describe the growth of a rectangular fracture can be written as:

(3.4) $$\begin{eqnarray}wb\frac{\text{d}L}{\text{d}t}=Q-2w\int _{0}^{L}v_{l}\,\text{d}x,\end{eqnarray}$$

where $w$ is the fracture’s width and $b$ is the activated fracture’s aperture.

To completely define the problem, constitutive relations are needed to describe the flux within the fracture and the leakage velocity. Since the length of the activated fracture is of the order of metres while the aperture is of the order of micro- to millimetres, the flux within the fracture for a Newtonian fluid can be described by the cubic law (Zimmerman & Bodvarsson Reference Zimmerman and Bodvarsson1996):

(3.5) $$\begin{eqnarray}q=-\frac{b^{3}}{12\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}x},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the fracturing fluid viscosity. The leaked fracturing fluid is assumed to displace the fluid saturating the rock, whose viscosity is much smaller than that of the injected fluid. The pressure at the interface between the two fluids is assumed to be equal to the pressure of the saturating fluid, $P_{f}$ . If capillary effects are important, i.e. when the surface tension per pore width is comparable with $P_{c}$ , one can account for it in the definition of $P_{f}$ . Given a pressure gradient imposed by the difference between the pressure of the fracturing fluid at the fracture surface and $P_{f}$ , the leaking fluid is assumed to propagate following Darcy’s law. This process is denoted as fluid seepage and can be described by (Howard & Fast Reference Howard and Fast1957):

(3.6) $$\begin{eqnarray}v_{l}=\unicode[STIX]{x1D700}_{m}\frac{\text{d}l_{w}}{\text{d}t}=\frac{k_{m}}{\unicode[STIX]{x1D707}}\frac{P-P_{f}}{l_{w}},\end{eqnarray}$$

where $l_{w}$ is the distance between the fracture surface and the interface, and $\unicode[STIX]{x1D700}_{m}$ and $k_{m}$ are the rock’s porosity and permeability, respectively. If the viscous pressure drop is much smaller than $(P_{c}-P_{f})$ , one can assume a quasi-steady fluid pressure at the surface of the fracture and, hence, the solution of (3.6) yields:

(3.7) $$\begin{eqnarray}v_{l}=C\sqrt{\frac{P-P_{f}}{t-t^{\prime }}},\end{eqnarray}$$

where $C=\sqrt{(k_{m}\unicode[STIX]{x1D700}_{m})/2\unicode[STIX]{x1D707}}$ . The fracturing fluid starts to leak from a point $x$ in space at $t^{\prime }$ ; the time at which the fracture tip reaches that point. The range of validity of the quasi-steady assumption is discussed in appendix A.

By substituting (3.7) and (3.5) into (3.1) and (3.4), the system of equations governing the growth of the fracture becomes:

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle b\frac{\text{d}L}{\text{d}t}=Q/w-2C\int _{0}^{L}\sqrt{\frac{P-P_{f}}{t-t^{\prime }(x)}}\,\text{d}x, & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}P}{\unicode[STIX]{x2202}x^{2}}=\frac{24\unicode[STIX]{x1D707}C}{b^{3}}\sqrt{\frac{P-P_{f}}{t-t^{\prime }(x)}}, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle P=P_{c},\quad \text{at }x=L(t), & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}x}=\frac{-12\unicode[STIX]{x1D707}Q}{wb^{3}},\quad \text{at }x=0, & \displaystyle\end{eqnarray}$$

with the initial condition $L(0)=0$ .

The following characteristic parameters will be used to non-dimensionalize the above system of equation:

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle t_{ch}=\frac{\unicode[STIX]{x1D707}b^{2}}{2k_{m}\unicode[STIX]{x1D700}_{m}(P_{c}-P_{f})}, & \displaystyle\end{eqnarray}$$

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle L_{ch}=\frac{Q\unicode[STIX]{x1D707}b}{2wk_{m}\unicode[STIX]{x1D700}_{m}(P_{c}-P_{f})}, & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle P_{ch}=\frac{6\unicode[STIX]{x1D707}^{2}Q^{2}}{w^{2}b^{2}k_{m}\unicode[STIX]{x1D700}_{m}(P_{c}-P_{f})}, & \displaystyle\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle l_{w}^{ch}=\frac{b}{\unicode[STIX]{x1D700}_{m}}, & \displaystyle\end{eqnarray}$$

where $t_{ch}$ is the time at which the activated fracture extends for a large enough length, $L_{ch}$ , for the leakage rate to significantly slow down the growth rate. For typical rocks with permeabilities ranging from $10^{-19}$ to $10^{-16}~\text{m}^{2}$ and $P_{c}-P_{f}$ of the order of MPa, the characteristic time ranges from a few seconds for water injected in a high permeability rock to years for a viscous fluid such as a cross-linked gel whose viscosity can range from 100 to 1000 cP injected in an ultra-low permeable rock (Montgomery Reference Montgomery and Jeffery2013). The characteristic pressure, $P_{ch}$ , is the viscous pressure drop required to drive the flow of the fracturing fluid over $L_{ch}$ . $l_{w}^{ch}$ is the characteristic penetration length into the rock matrix. In order to avoid interference between fractures, $l_{w}^{ch}$ must be much smaller than the average spacing between two fractures. Since the aperture of typical fractures is of the order of a few millimetres, the penetration length will be of the order of a few centimetres which is much smaller than a typical $O$ (10) m fracture spacing.

Now, let $\bar{t}=t/t_{ch}$ , $\bar{P}=(P-P_{c})/P_{ch}$ , $\bar{L}=L/L_{ch}$ and substitute these dimensionless variables into (3.8) through (3.11) to get:

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\bar{L}}{\text{d}\bar{t}}=1-\int _{0}^{\bar{L}}\sqrt{\frac{\unicode[STIX]{x1D6E5}_{pc}\bar{P}+1}{\bar{t}-\bar{t^{\prime }}}}\,\text{d}\bar{x}, & \displaystyle\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}\bar{P}}{\unicode[STIX]{x2202}\bar{x^{2}}}=\sqrt{\frac{\unicode[STIX]{x1D6E5}_{pc}\bar{P}+1}{\bar{t}-\bar{t^{\prime }}}}, & \displaystyle\end{eqnarray}$$

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{P}=0,\quad \text{at }x=\bar{L}(\bar{t}), & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{P}}{\unicode[STIX]{x2202}\bar{x}}=-1,\quad \text{at }x=0, & \displaystyle\end{eqnarray}$$

with the initial condition of zero length at $t=0$ . $\unicode[STIX]{x1D6E5}_{pc}\equiv P_{ch}/(P_{c}-P_{f})$ measures the importance of pressure variation within the fracture on the leakage rate. For instance, the leakage rate becomes independent of the viscous pressure drop along the fracture when $\unicode[STIX]{x1D6E5}_{pc}\rightarrow 0$ .

Before introducing the numerical solution of the above system of equations for different values of $\unicode[STIX]{x1D6E5}_{pc}$ , several limiting cases will be discussed. We will show that when $\unicode[STIX]{x1D6E5}_{pc}=0$ , two similarity solutions for the fluid pressure can be obtained depending on the importance of the leakage rate in affecting the growth of the fracture. Thereafter, we show that the similarity exponents found in this case are retained even when one further simplifies the leakage rate by assuming that $t^{\prime }=0$ and $\unicode[STIX]{x1D6E5}_{pc}=0$ .

3.2 Similarity solutions when the pressure variation in the fracture is negligible

When $\unicode[STIX]{x1D6E5}_{pc}=0$ , the leakage velocity is not sensitive to fluid pressure variations in the fracture and therefore (3.16) and (3.17) decouple. In this case, one can obtain an analytical solution for the fracture length and identify two regimes where the growth behaves in a self-similar fashion. The complete solution of (3.16), when $\unicode[STIX]{x1D6E5}_{pc}=0$ , is given by:

(3.20) $$\begin{eqnarray}\bar{L}=\frac{1}{\unicode[STIX]{x03C0}}\left[\text{e}^{\unicode[STIX]{x03C0}\bar{t}}\text{erfc}\left(\sqrt{\unicode[STIX]{x03C0}\bar{t}}\right)+2\sqrt{\bar{t}}-1\right]\!.\end{eqnarray}$$

This analytical solution is valid when the ratio of the fluid velocity within the fracture to the characteristic leakage velocity is much larger than 1 but much smaller than the square root of the ratio of the flow resistance within the rock matrix to the resistance to flow within the fracture, i.e.

(3.21) $$\begin{eqnarray}1\ll \frac{\displaystyle \left(\frac{Q}{wb}\right)}{\displaystyle \left(\frac{(P_{c}-P_{f})k_{m}\unicode[STIX]{x1D700}_{m}}{\unicode[STIX]{x1D707}b}\right)}\ll \frac{b}{\sqrt{k_{m}\unicode[STIX]{x1D700}_{m}}}.\end{eqnarray}$$

Derivations of (3.20) and (3.21) are detailed in appendix A. Moreover, equation (3.20) has been identified by Carter to describe the growth of a tensile rectangular fracture where the fluid pressure is uniform inside the fracture (Howard & Fast Reference Howard and Fast1957). In his model, the fluid pressure is assumed to be equal to the pressure required to move the two surfaces apart creating a pressure-driven tensile fracture whereas this model assumes that the pressure is equal to the critical pressure for hydroshear slippage of a natural fracture.

The two regimes where the fracture length grows as a power law of time and (3.17) admits a similarity solution correspond to cases where most of the injected fluid goes toward growing the fracture and toward supplying the leaked fluid, respectively. When $\bar{t}\ll 1$ , the leakage terms can be neglected and (3.16)–(3.19) can be transformed to:

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\bar{L}}{\bar{t}}=\unicode[STIX]{x1D702}_{max}=1, & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}^{2}}=0, & \displaystyle\end{eqnarray}$$

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=0,\quad \text{at }\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{max}, & \displaystyle\end{eqnarray}$$

(3.25) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}}=-1,\quad \text{at }\unicode[STIX]{x1D702}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D702}=\bar{x}/\bar{t}$ , and $\unicode[STIX]{x1D719}=\bar{P}/\bar{t}$ . By solving the above system of equations, one can see that the fracture length grows as $\bar{L}=\bar{t}$ , and the pressure profile can be written as:

(3.26) $$\begin{eqnarray}\bar{P}=\bar{t}\left(1-\frac{\bar{x}}{\bar{t}}\right).\end{eqnarray}$$

The power law solution of the fracture’s length can also be obtained using the asymptotic behaviour of (3.20) when $\bar{t}\rightarrow 0$ .

As leakage becomes important, the growth rate slows down and when the fracture length is much larger than $L_{ch}$ , i.e.  $\bar{t}\gg 1$ , most of the injected fluid leaks off. In this regime, the left-hand term in (3.16) can be neglected and (3.16)–(3.19) can be transformed to:

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{max}\int _{0}^{1}\frac{1}{\sqrt{1-\unicode[STIX]{x1D712}^{2}}}\,\text{d}\unicode[STIX]{x1D712}=1, & \displaystyle\end{eqnarray}$$

(3.28) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D712}^{2}}=\frac{\unicode[STIX]{x1D702}_{max}^{2}}{\sqrt{1-\unicode[STIX]{x1D712}^{2}}}, & \displaystyle\end{eqnarray}$$

(3.29) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=0,\quad \text{at }\unicode[STIX]{x1D712}=1, & \displaystyle\end{eqnarray}$$

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D712}}=-\unicode[STIX]{x1D702}_{max},\quad \text{at }\unicode[STIX]{x1D712}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D702}_{max}=\bar{L}/\sqrt{\bar{t}}$ , $\unicode[STIX]{x1D702}=\bar{x}/\sqrt{\bar{t}}$ and $\unicode[STIX]{x1D719}=\bar{P}/\sqrt{\bar{t}}$ . Since $\bar{t^{\prime }}$ is the inverse of $\bar{L}(\bar{t})$ , $\bar{t^{\prime }}/\bar{t}=\unicode[STIX]{x1D712}^{2}$ where $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D702}/\unicode[STIX]{x1D702}_{max}$ . From (3.27), $\unicode[STIX]{x1D702}_{max}=2/\unicode[STIX]{x03C0}$ and, therefore, the length of the fracture grows as $\bar{L}=(2/\unicode[STIX]{x03C0})\sqrt{\bar{t}}$ . This can also be obtained using the asymptotic behaviour of (3.20) when $\bar{t}\rightarrow \infty$ . The solution of (3.28) shows that the self-similar pressure profile can be written as:

(3.31) $$\begin{eqnarray}\bar{P}=\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{2}\sqrt{\bar{t}}\left[\sqrt{1-\left(\frac{\unicode[STIX]{x03C0}\bar{x}}{2\sqrt{\bar{t}}}\right)^{2}}+\frac{\unicode[STIX]{x03C0}\bar{x}}{2\sqrt{\bar{t}}}\sin ^{-1}\left(\frac{\unicode[STIX]{x03C0}\bar{x}}{2\sqrt{\bar{t}}}\right)\right]-\bar{x}.\end{eqnarray}$$

3.3 Homogeneous leakage rate approximation

Next, we will analyse the case where $\bar{t^{\prime }}(x)=0$ and $\unicode[STIX]{x1D6E5}_{pc}=0$ . In this case, not only is the pressure variation within the fracture negligible but the leakage velocity is also homogeneous along the fracture and only depends on the age of the fracture, i.e. the time since it starts to grow. This approximation will prove useful when modelling the growth of the two-dimensional network of fractures when the time scale to grow a fracture is much smaller than the time scale to grow a whole network of activated fractures. We will show that simplifying the leakage rate in this way does not change the similarity scaling but only the numerical prefactor.

When the leakage rate is independent of position and pressure, the non-dimensional system of equations governing the fracture growth becomes:

(3.32) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\bar{L}}{\text{d}\bar{t}}=1-\frac{\bar{L}}{\sqrt{\bar{t}}}, & \displaystyle\end{eqnarray}$$

(3.33) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}\bar{P}}{\unicode[STIX]{x2202}\bar{x^{2}}}=\frac{1}{\sqrt{\bar{t}}}, & \displaystyle\end{eqnarray}$$

(3.34) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{P}=0,\quad \text{at }x=\bar{L}(\bar{t}), & \displaystyle\end{eqnarray}$$

(3.35) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{P}}{\unicode[STIX]{x2202}\bar{x}}=-1,\quad \text{at }x=0 & \displaystyle\end{eqnarray}$$

and $\bar{L}(0)=0$ .

By solving (3.32), one can find the solution for the fracture length that is given by:

(3.36) $$\begin{eqnarray}\bar{L}=\sqrt{\bar{t}}+{\textstyle \frac{1}{2}}\text{e}^{-2\sqrt{\bar{t}}}-{\textstyle \frac{1}{2}}.\end{eqnarray}$$

The self-similar solution in the regime where $\bar{t}\ll 1$ is the same as derived in § 3.2 since the difference introduced in this section is in the leakage rate which is negligible when $\bar{t}\ll 1$ . However, in the fluid loss dominated regime, the self-similar transform of (3.32) shows that $\bar{L}=\sqrt{\bar{t}}$ . Furthermore, the ordinary differential equation that describes the self-similar pressure profile, $\unicode[STIX]{x1D719}$ is given by:

(3.37) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}^{2}}=1, & \displaystyle\end{eqnarray}$$

(3.38) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=0,\quad \text{at }\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{max}, & \displaystyle\end{eqnarray}$$

(3.39) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}}=-1,\quad \text{at }\unicode[STIX]{x1D702}=0, & \displaystyle\end{eqnarray}$$

using the same similarity scaling used to derive (3.28). $\unicode[STIX]{x1D702}_{max}$ in this case is equal to 1 and the solution of the pressure profile in the leakage dominant regime is given by:

(3.40) $$\begin{eqnarray}\bar{P}=\sqrt{\bar{t}}\left(\frac{\bar{x}^{2}}{2\bar{t}}-\frac{\bar{x}}{\sqrt{\bar{t}}}+\frac{1}{2}\right).\end{eqnarray}$$

Figure 3 compares the growth of a fracture in the case of a uniform leakage velocity with the case where the leakage rate is a functional that depends on $t^{\prime }(x)$ . In both cases, $\unicode[STIX]{x1D6E5}_{pc}=0$ . As expected, neglecting $\bar{t}^{\prime }$ when calculating the leakage velocity results in over-predicting the fracture length, as shown in figure 3(a). Furthermore, figure 3(b) compares the self-similar pressure profile obtained from (3.31) and (3.40). Underestimating the leakage rate when assuming a uniform leakage velocity increases the required viscous pressure drop and, hence, the injection pressure is higher.

Figure 3. Comparison of the fracture growth with constant and space-dependent leakage velocity. In both cases, $\unicode[STIX]{x1D6E5}_{pc}=0$ . (a) The dimensionless fracture length versus the dimensionless time. The blue dashed line represents the solution for the fracture length in the case of a constant leakage, equation (3.36), while the red solid curve corresponds to the case where the leakage is space dependent, given by (3.20). (b) The similarity pressure profiles $\unicode[STIX]{x1D719}=\bar{P}/\sqrt{\bar{t}}$ for the single fracture model in the leakage dominated regime. The red solid line corresponds to the space-dependent leakage case, equation (3.31), and the blue dashed line represents the homogenous leakage case, equation (3.40).

To quantify the errors introduced when simplifying the leakage rate by neglecting its dependence on $t^{\prime }$ , consider a portion of the growing fracture with a length, $L_{joint}$ , that is much smaller than $L_{ch}$ where the self-similar solution in the short-time regime can be used to calculate $\hat{t}$ . $\hat{t}$ is the time at which the fluid front reaches $L_{joint}$ , $\hat{t}=t^{\prime }(L_{joint})$ . The homogeneous leakage rate predicts that the leakage rate of the whole joint will decay as $\hat{t}/\sqrt{t}$ while it decays as $2\sqrt{t}(1-\sqrt{1-\hat{t}/t})$ when one accounts for the dependence of the leakage rate on $t^{\prime }$ . Therefore, one can see that the homogeneous leakage rate represents the asymptotic solution of the integral leakage rate as $\hat{t}/t\rightarrow 0$ . Thus, this simplification applies when the time scale to grow a one-dimensional network of multiple fractures is much longer than the time scale to completely activate a single fracture of a finite length.

Figure 4. Comparison of the fracture growth with constant and space-dependent leakage velocity. (a) The solution for the length of a propagating single fracture for different values of $\unicode[STIX]{x1D6E5}_{pc}$ . As $\unicode[STIX]{x1D6E5}_{pc}$ increases the solution deviates from the analytical solution, which is obtained for the case of $\unicode[STIX]{x1D6E5}_{pc}=0$ , at large times when the leakage starts to become important. (b) The pressure profile inside the growing fracture for different values of $\unicode[STIX]{x1D6E5}_{pc}$ at the same fracture length. The blue thicker lines correspond to the case where $\unicode[STIX]{x1D6E5}_{pc}=1$ while the red thinner ones correspond to the pressure profile at different times for the case where $\unicode[STIX]{x1D6E5}_{pc}=0$ .

3.4 Numerical solution of the integro-differential equation

To capture the effects of the fluid pressure on the leakage rate and consequently on the fracture growth, equations (3.16)–(3.19) are solved numerically for different values of $\unicode[STIX]{x1D6E5}_{pc}$ . The solution is compared with the analytical solution (3.20) for the case where $\unicode[STIX]{x1D6E5}_{pc}=0$ . As shown in figure 4(a), as $\unicode[STIX]{x1D6E5}_{pc}$ is increased, the activated fracture grows at a slower rate. For smaller values of $\unicode[STIX]{x1D6E5}_{pc}$ , the deviation is less pronounced and it becomes important as $\unicode[STIX]{x1D6E5}_{pc}$ becomes larger than $O(1)$ . As the leakage starts to play a role, the effects of approximating the leakage rate as if it were independent of the fluid pressure inside the fractures become significant. At larger $\unicode[STIX]{x1D6E5}_{pc}$ , the leakage rate is higher and therefore the fracture propagates at a slower rate than would be predicted when the effects of pressure variation on the leakage velocity are neglected.

Additionally, the self-similar behaviour of the pressure profile is lost when the leakage rate depends on the local value of the fluid pressure. Figure 4(b) shows a comparison of the transient pressure profile for $\unicode[STIX]{x1D6E5}_{pc}=1$ and $\unicode[STIX]{x1D6E5}_{pc}=0$ at various fracture lengths. As expected, the initial pressure profile, when leakage is negligible, is not affected by the value of $\unicode[STIX]{x1D6E5}_{pc}$ but deviations occur as the fracture becomes long enough for leakage to become important. Leakage reduces the velocity of the injected fluid within the fracture and therefore decreases the pressure drop across the fracture. Since the pressure at the fracture tip is the same for both values of $\unicode[STIX]{x1D6E5}_{pc}$ , the injection pressure is lower when accounting for the effects of fluid pressure on the leakage velocity.

In summary, we have shown that an activated natural fracture’s length grows, when the fluid pressure’s effect on the leakage velocity is negligible, as a power law of time in two different regimes corresponding to negligible and dominant leakage rate. Furthermore, we have shown that homogenization of the leakage velocity along the fracture’s surface does not change the self-similarity behaviour but only the accuracy in predicting the length of the activated fractures. Finally, the homogenization of the leakage velocity can be used when the weak plane length is much smaller than  $L_{ch}$ . Although a deviation from the long-time similarity scaling occurs if $\unicode[STIX]{x1D6E5}_{pc}\geqslant O(1)$ , the effect of pressure variation in the fracture on the leakage rate is negligible for $\unicode[STIX]{x1D6E5}_{pc}\leqslant 0.2$ . These results will prove useful in the next section when we simplify the two dimensional model while maintaining the self-similar growth behaviour.

4.1 Derivation of the governing equations

Similar to the single fracture case, the system of equations given by (2.8) and (2.9) will be written in terms of the fluid pressure. The flux of the fluid is described using Darcy’s flow where the permeability depends on the fluid pressure. If the rock layer thickness is $H$ and $w/H=O(1)$ where $w$ is the average width of the fractures, the fracturing process can be described in two dimensions. However, if the rock layer has an unbounded thickness, the cluster grows in three dimensions. In cylindrical and spherical coordinates, the radially symmetric Darcy’s law is written as:

(4.1) $$\begin{eqnarray}\boldsymbol{q}=-\frac{K}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}r}\boldsymbol{e}_{r},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the fracturing fluid viscosity, and $K$ is the local permeability. The growth of the cluster in two and three dimensions is assumed to be radially symmetric since the distribution of critical pressures is considered to be statically homogeneous and isotropic.

In § 3.3, it has been shown that the rate of leakage from a growing rectangular fracture can be estimated by:

(4.2) $$\begin{eqnarray}\int _{0}^{L}v_{l}\,\text{d}x\approx \frac{2CLw\sqrt{P_{c}-P_{f}}}{\sqrt{t-t^{\prime }}},\end{eqnarray}$$

where $t^{\prime }$ is the time at which the fracture starts to grow. If a single fracture is considered, $t^{\prime }=0$ . This estimated leakage rate ignores the variation in the leakage velocity within the fracture and has been shown to be fairly accurate if the length of the fracture is much smaller than $L_{ch}$ and if $\unicode[STIX]{x1D6E5}_{pc}=0$ . If a representative number of rectangular fractures are continuously being activated within the cluster as it grows, one can replace the surface area in (4.2) by the surface area of the activated fractures. Thus, the local leakage rate per volume of the medium is given by:

(4.3) $$\begin{eqnarray}V_{l}=C_{l}\int _{0}^{S(r,t)}\frac{\text{d}s}{\sqrt{t-t^{\prime }(s,r)}},\end{eqnarray}$$

where $C_{l}=(2\unicode[STIX]{x1D700}/b)C\sqrt{P_{cc}-P_{f}}$ , and $t^{\prime }(s,r)$ is the time at which the local fraction of the activated fractures is equal to $s$ . The leakage rate within the activated fractures is assumed to be independent of the fluid pressure. This applies when the pressure drop across the network is much smaller than $(P_{cc}-P_{f})$ . This applies for fluids with moderate viscosity. However, when extremely viscous fluids are injected at a high flow rate, one needs to modify (4.3) by using (3.7) to account for the effects of fluid pressure on the leakage velocity. In this paper, such extreme cases are not considered since most fracturing processes use water derivatives that have modest viscosities. Typically, extremely viscous fluids are used to hold proppants which are not needed when activating pre-existing fractures (Montgomery Reference Montgomery and Jeffery2013).

Substituting (4.1) and (4.3) into (2.8) and (2.9), the governing equation describing the growth of the cluster in two dimensions becomes:

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}t}=\frac{1}{\unicode[STIX]{x1D707}r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(rK\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}r}\right)-C_{l}\int _{0}^{S}\frac{\text{d}s}{\sqrt{t-t^{\prime }(s,r)}}, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\left(\int _{0}^{R}Sr\,\text{d}r\right)=\frac{Q_{0}t^{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x1D714}}-C_{l}\int _{0}^{R}r\,\text{d}r\int _{0}^{S}\frac{\text{d}s}{\sqrt{t-t^{\prime }(s,r)}}. & \displaystyle\end{eqnarray}$$

For a time-dependent flow rate, $Q_{0}t^{\unicode[STIX]{x1D6FC}}$ , a boundary condition at the injection well can be written as:

(4.6) $$\begin{eqnarray}\lim _{r\rightarrow ~0}\frac{K}{\unicode[STIX]{x1D707}}r\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}r}=-\frac{Q_{0}}{\unicode[STIX]{x1D714}}t^{\unicode[STIX]{x1D6FC}},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}$ is a geometric factor that depends on the dimension of the cluster. For $D=2$ , $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}H$ and $\unicode[STIX]{x1D714}=4\unicode[STIX]{x03C0}$ when considering a three-dimensional cluster. Equation (4.6) applies when the cluster radius is much larger than the radius of the injection well. Since a percolating network forms when the fluid pressure is at least equal to $P_{cc}$ , a boundary condition at the edge of the cluster can be written as:

(4.7) $$\begin{eqnarray}P=P_{cc}\quad \text{on }r=R(t).\end{eqnarray}$$

This boundary condition ensures that a percolating network of activated fractures forms which in turn will allow the fluid to flow through the primary porosity.

Equation (4.4) to (4.5) is a complete set of equations coupling fluid flow with the growth of the network of activated fractures. The constitutive relations needed to relate $S$ and $K$ with the fluid pressure are given by (2.5) and (2.6), respectively. Finally, the required initial condition is that the radius is zero at $t=0$ .

To non-dimensionalize the governing equations, the following characteristic parameters will be used:

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle t_{ch}=\left(\frac{\unicode[STIX]{x1D700}}{C_{l}}\right)^{2}, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle R_{ch}=(C_{f}^{\unicode[STIX]{x1D716}+1}Q_{f}^{\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD}}t_{ch}^{-E_{2}})^{1/E_{1}}, & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle P_{ch}=\left(\frac{Q_{f}t_{ch}^{\unicode[STIX]{x1D6FC}}}{R_{ch}^{D-2}}\right)^{1/(\unicode[STIX]{x1D716}+1)}, & \displaystyle\end{eqnarray}$$

where $C_{f}=k_{0}k_{k}k_{f}^{\unicode[STIX]{x1D716}-\unicode[STIX]{x1D6FD}}/k_{s}\unicode[STIX]{x1D700}\unicode[STIX]{x1D707}$ and $Q_{f}=Q_{0}\unicode[STIX]{x1D707}/k_{0}k_{k}k_{f}^{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D714}$ . The exponents $E_{1}$ and $E_{2}$ are functions of the percolation exponents $D$ and $\unicode[STIX]{x1D6FC}$ and their values are given in table 1. Similar to the single fracture growth problem, $t_{ch}$ is the time required for the cluster to grow to a radius $R_{ch}$ at which the total leakage rate starts to play an important role in slowing down the cluster propagation rate. $P_{ch}$ is the characteristic flow induced pressure drop required to drive the flow over the radius $R_{ch}$ . One should note that the assumptions used to obtain the leakage rate in (4.3) and (2.3) are valid when $P_{ch}/(P_{cc}-P_{f})\ll 1$ .

Table 1. Summary of exponents used and their values in two and three dimensions. $\unicode[STIX]{x1D6FD}$  is equal to $5/36$ when $D=2$ and equal to $0.42$ in three dimensions. $\unicode[STIX]{x1D716}=1.3$ in two dimensions and 2.0 in three dimensions.

Now, let $\unicode[STIX]{x1D70F}=t/t_{ch},\unicode[STIX]{x1D70F}^{\prime }=t^{\prime }(s,r)/t_{ch},p=(P-P_{cc})/P_{ch},\unicode[STIX]{x1D70C}_{max}=R/R_{ch}$ and $\unicode[STIX]{x1D70C}=r/R_{ch}$ . When substituting the percolation relations given in (2.4)–(2.6) into (4.4)–(4.5), the non-dimensionalized form of the equations can be written as:

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p^{\unicode[STIX]{x1D6FD}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\left(\unicode[STIX]{x1D70C}p^{\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right)-\unicode[STIX]{x1D6FD}\int _{0}^{p}\frac{\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D6FD}-1}}{\sqrt{\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}^{\prime }}}\,\text{d}\unicode[STIX]{x1D701}, & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\left(\int _{0}^{\unicode[STIX]{x1D70C}_{max}}p^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D70C}\,\text{d}\unicode[STIX]{x1D70C}\right)=\unicode[STIX]{x1D70F}^{\unicode[STIX]{x1D6FC}}-\unicode[STIX]{x1D6FD}\int _{0}^{\unicode[STIX]{x1D70C}_{max}}\,\text{d}\unicode[STIX]{x1D701}\int _{0}^{p}\frac{\unicode[STIX]{x1D701}^{\unicode[STIX]{x1D6FD}-1}}{\sqrt{\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}^{\prime }}}\unicode[STIX]{x1D70C}\,\text{d}\unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle p=0\quad \text{on }\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{max}, & \displaystyle\end{eqnarray}$$

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{\unicode[STIX]{x1D70C}\rightarrow ~0}\unicode[STIX]{x1D70C}p^{\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}=-\unicode[STIX]{x1D70F}^{\unicode[STIX]{x1D6FC}}. & \displaystyle\end{eqnarray}$$

Before introducing the full solution of the differential equation, two regimes where similarity solutions arise will be introduced. Analogous to the single fracture model, the pressure-driven flow dominates over the total leakage rate when $\unicode[STIX]{x1D70F}\ll 1$ due to the small surface area the fluid can leak through. By neglecting the leakage term in (4.11) and (4.12), a similarity solution can be obtained. This regime is denoted as the short-time regime, although one should note that the time needs to be long enough for the radius of the cluster to grow much larger than $\unicode[STIX]{x1D709}_{0}$ so that the continuum description applies. The other case in which a similarity solution arises is when one can neglect the left-hand term in (4.11) and (4.12). This corresponds to a long-time regime where $\unicode[STIX]{x1D70F}\gg 1$ in which the total leakage rate dominates and most of the fluid injected leaks off. The time where $\unicode[STIX]{x1D70F}\sim 1$ corresponds to the crossover between the two similarity scaling and is reached when the leakage first starts to play a role in fluid transport.

In both similarity solutions, the radius of the cluster grows as a power law of time. For the short-time regime:

(4.15) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{max}=\unicode[STIX]{x1D702}_{max}^{s}\unicode[STIX]{x1D70F}^{\unicode[STIX]{x1D6FE}_{s}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{s}=((1+\unicode[STIX]{x1D6FC})/2)-\unicode[STIX]{x1D6FC}E_{3}$ and in the long-time regime:

(4.16) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{max}=\unicode[STIX]{x1D702}_{max}^{l}\unicode[STIX]{x1D70F}^{\unicode[STIX]{x1D6FE}_{l}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{l}=((1+2\unicode[STIX]{x1D6FC})/4)-\unicode[STIX]{x1D6FC}E_{3}$ . $\unicode[STIX]{x1D702}_{max}^{s}$ and $\unicode[STIX]{x1D702}_{max}^{l}$ are proportionality constants in the short and long-time regimes, respectively, and depend on the value of $\unicode[STIX]{x1D6FC}$ . The general expression and values of $E_{3}$ , $\unicode[STIX]{x1D6FE}_{s}$ and $\unicode[STIX]{x1D6FE}_{l}$ in both dimensions and for different $\unicode[STIX]{x1D6FC}$ values are given in table 1.

4.2 Short-time similarity solution

To find the similarity solution for the pressure profile in the short time regime, let $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70F}^{-\unicode[STIX]{x1D6FE}_{s}}$ and

(4.17) $$\begin{eqnarray}\unicode[STIX]{x1D719}=p\unicode[STIX]{x1D70F}^{-\unicode[STIX]{x1D6FF}_{s}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{s}=\unicode[STIX]{x1D6FC}/(\unicode[STIX]{x1D716}+1)$ ; a general expression for $\unicode[STIX]{x1D6FF}_{s}$ in $D$ dimensions is given in table 1. By substituting these scaled variables into (4.11)–(4.14) and neglecting the leakage terms, an ordinary differential equation is obtained:

(4.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}_{s}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D6FE}_{s}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FD}-1}\unicode[STIX]{x1D702}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}}=\frac{1}{\unicode[STIX]{x1D702}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D702}}\left(\unicode[STIX]{x1D702}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D716}}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}}\right), & \displaystyle\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\unicode[STIX]{x1D702}_{max}^{s}}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D702}\,\text{d}\unicode[STIX]{x1D702}=\frac{1}{\unicode[STIX]{x1D6FC}+1}, & \displaystyle\end{eqnarray}$$

(4.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=0\quad \text{on }\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{max}^{s}, & \displaystyle\end{eqnarray}$$

(4.21) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{\unicode[STIX]{x1D702}\rightarrow ~0}\unicode[STIX]{x1D702}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D716}}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}}=-1. & \displaystyle\end{eqnarray}$$

By inspecting the equations, one can find two asymptotic solutions in the limits where $\unicode[STIX]{x1D702}\rightarrow 0$ and when $\unicode[STIX]{x1D702}\rightarrow \unicode[STIX]{x1D702}_{max}^{s}$ . The asymptotic behaviour of $\unicode[STIX]{x1D719}$ near the injection point is given by:

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D719}=\left[(\unicode[STIX]{x1D716}+1)\ln \left(\frac{1}{\unicode[STIX]{x1D702}}\right)\right]^{1/(\unicode[STIX]{x1D716}+1)}\end{eqnarray}$$

and the behaviour of $\unicode[STIX]{x1D719}$ as $\unicode[STIX]{x1D702}\rightarrow \unicode[STIX]{x1D702}_{max}^{s}$ is given by:

(4.23) $$\begin{eqnarray}\unicode[STIX]{x1D719}=\left[\unicode[STIX]{x1D6FE}_{s}(\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D702}_{max}^{s}(\unicode[STIX]{x1D702}_{max}^{s}-\unicode[STIX]{x1D702})\right]^{1/(\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD})}.\end{eqnarray}$$

From (2.4), one can see that $\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D719}^{-\unicode[STIX]{x1D708}}$ . Therefore, $\unicode[STIX]{x1D709}$ is expected to approach $\unicode[STIX]{x1D709}_{0}$ near the injection well and $\unicode[STIX]{x1D709}\rightarrow \infty$ near the cluster edge. The criterion given by (2.7) is thus violated at both the centre and edge of the cluster. A detailed discussion about the thickness of the regions where the continuum approximation or use of percolation scalings break down is presented in appendix D. To find the solution for $\unicode[STIX]{x1D719}$ , a shooting method, described in appendix C, was used.

Figure 5(a) shows the short-time solution for $\unicode[STIX]{x1D719}$ for different values of $\unicode[STIX]{x1D6FC}$ . An approximate empirical function that fits the numerical solution for $\unicode[STIX]{x1D719}$ in this regime and has the correct asymptotic behaviours near the injection well and the cluster’s edge is given by:

(4.24) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{app}=C_{\unicode[STIX]{x1D6FC}}\frac{\displaystyle \left[(\unicode[STIX]{x1D716}+1)\ln \left(\frac{\unicode[STIX]{x1D702}_{max}^{s}}{\unicode[STIX]{x1D702}}\right)\right]^{1/(\unicode[STIX]{x1D716}+1)}+C_{\unicode[STIX]{x1D6FC}}^{-1}\unicode[STIX]{x1D702}^{3}}{\displaystyle C_{\unicode[STIX]{x1D6FC}}-C_{1}^{-1}+\left(C_{1}\left(1-\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{max}^{s}}\right)^{1/(\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD})}\right)^{-1}},\end{eqnarray}$$

where $C_{1}=[\unicode[STIX]{x1D6FE}_{s}(\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D702}_{max}^{s^{2}}]^{1/(\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD})}$ . $C_{\unicode[STIX]{x1D6FC}}$ is a fitting parameter that depends on $\unicode[STIX]{x1D6FC}$ and it is given by:

(4.25) $$\begin{eqnarray}C_{\unicode[STIX]{x1D6FC}}=a_{\unicode[STIX]{x1D6FC}}+b_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D6FC}+1)+c_{\unicode[STIX]{x1D6FC}}\ln (\unicode[STIX]{x1D6FC}+1)^{2},\end{eqnarray}$$

where $a_{\unicode[STIX]{x1D6FC}}=12.39,b_{\unicode[STIX]{x1D6FC}}=-5.42,c_{\unicode[STIX]{x1D6FC}}=6.93$ . The maximum relative absolute error when fitting the approximate function with the numerical solution for $\unicode[STIX]{x1D719}$ is less than 10 %. This expression along with relations for permeability and porosity in terms of pressure, as will be shown later in § 5, can be used as a constitutive relation to model gas and heat transport in hydraulically fractured rocks when fracturing fluid leakage was negligible during the stimulation process.

Figure 5. Comparison of the self-similar pressure behaviour for different values of $\unicode[STIX]{x1D6FC}$ in the short- and long-time regimes. (a) The similarity solution for different values of $\unicode[STIX]{x1D6FC}$ in the regime where $\unicode[STIX]{x1D70F}\ll 1$ . $\unicode[STIX]{x1D702}_{max}^{s}$ corresponds the value of $\unicode[STIX]{x1D702}$ when $\unicode[STIX]{x1D719}=0$ . As $\unicode[STIX]{x1D6FC}$ increases, $\unicode[STIX]{x1D702}_{max}^{s}$ increases. (b) The similarity solution for different values of $\unicode[STIX]{x1D6FC}$ in the regime where $\unicode[STIX]{x1D70F}\gg 1$ . $\unicode[STIX]{x1D702}_{max}^{l}$ is when $\unicode[STIX]{x1D719}=0$ and it decreases as $\unicode[STIX]{x1D6FC}$ increases.

Given the solution of $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D70C}_{max}$ and the characteristic parameters, one can calculate the porosity and permeability profiles of the network. Detailed discussion of the effects of different operating conditions on the properties of the network is presented in § 5.

4.3 Long-time similarity solution

For $\unicode[STIX]{x1D70F}\gg 1$ , a similarity solution can be obtained by balancing the pressure-driven flux with the leakage term in (4.11). This states that most of the injected liquid will be used to supply the leakage into the rock matrix and only a small amount of the injected fluid will be used to activate new fractures, so that we can neglect the left-hand term in (4.11) and (4.12). Since (4.11) is difficult to solve numerically, further assumptions will be used. In this section, we will derive the self-similar solution for a simplified problem and show that the required assumptions do not change the similarity scaling but only the quantitative accuracy in predicting the propagation rate. The derivation of the self-similar integro-differential equations without simplification is presented in appendix B.

In the single fracture model, it has been shown that when assuming a uniform leakage velocity throughout the fracture, the self-similar scalings were retained. Likewise, the dependence of the local leakage rate per surface area, in the multiple fracture model, on the time at which each fracture is activated will be ignored. We will assume that all fractures that are located at a particular position will have the same leakage velocity regardless of the time at which each fracture was activated. Their leakage velocity will depend on the time at which the fluid front reaches that position. This assumption will still preserve the scalings obtained from analysing (4.11) in the leakage dominant regime but is not expected to produce a quantitatively accurate solution. Since this assumption under-predicts the leakage rate, the predicted radius from the new similarity solution will be larger than the one obtained if one solves the actual system of equations for the full leakage rate.

Neglecting the dependence of the leakage velocity of each activated fracture on its activation time, equation (4.3) that describes the leakage rate can be approximated as:

(4.26) $$\begin{eqnarray}\int _{0}^{S(r,t)}\frac{\text{d}s}{\sqrt{t-t^{\prime }(s,r)}}\approx \frac{S(r,t)}{\sqrt{t-t^{\prime }(r)}}.\end{eqnarray}$$

In this equation, $t^{\prime }(r)$ is the time when the fracturing fluid first reaches a point in space $r$ . Using this approximated leakage term, equations (4.11) and (4.12) become:

(4.27) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p^{\unicode[STIX]{x1D6FD}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\left(\unicode[STIX]{x1D70C}p^{\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right)-\frac{p^{\unicode[STIX]{x1D6FD}}}{\sqrt{\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}^{\prime }}} & \displaystyle\end{eqnarray}$$

(4.28) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\left(\int _{0}^{\unicode[STIX]{x1D70C}_{max}}p^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D70C}\,\text{d}\unicode[STIX]{x1D70C}\right)=\unicode[STIX]{x1D70F}^{\unicode[STIX]{x1D6FC}}-\int _{0}^{\unicode[STIX]{x1D70C}_{max}}\frac{p^{\unicode[STIX]{x1D6FD}}}{\sqrt{\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}^{\prime }}}\unicode[STIX]{x1D70C}\,\text{d}\unicode[STIX]{x1D70C}. & \displaystyle\end{eqnarray}$$

To obtain the self-similar behaviour in the long-time regime, substitute $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70F}^{-\unicode[STIX]{x1D6FE}_{l}}$ , equation (4.16) and

(4.29) $$\begin{eqnarray}\unicode[STIX]{x1D719}=p\unicode[STIX]{x1D70F}^{-\unicode[STIX]{x1D6FF}_{l}},\end{eqnarray}$$

into (4.27)–(4.28) and neglect the left-hand terms. The resulting ordinary differential equation in two dimensions is written as:

(4.30) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FD}}}{\displaystyle \sqrt{1-\left(\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{max}^{l}}\right)^{1/\unicode[STIX]{x1D6FE}_{l}}}}=\frac{1}{\unicode[STIX]{x1D702}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D702}}\left(\unicode[STIX]{x1D702}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D716}}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D702}}\right), & \displaystyle\end{eqnarray}$$

(4.31) $$\begin{eqnarray}\displaystyle & \displaystyle 1=\int _{0}^{\unicode[STIX]{x1D702}_{max}^{l}}\frac{\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D702}\,\text{d}\unicode[STIX]{x1D702}}{\displaystyle \sqrt{1-\left(\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{max}^{l}}\right)^{1/\unicode[STIX]{x1D6FE}_{l}}}}, & \displaystyle\end{eqnarray}$$

where the boundary conditions are given by (4.20) and (4.21). $\unicode[STIX]{x1D6FF}_{l}=\unicode[STIX]{x1D6FC}/(\unicode[STIX]{x1D716}+1)$ and its general expression and value are given in table 1.

The asymptotic solution for $\unicode[STIX]{x1D719}$ near the injection point is given by (4.22) while $\unicode[STIX]{x1D719}$ as $\unicode[STIX]{x1D702}\rightarrow \unicode[STIX]{x1D702}_{max}^{l}$ is given by:

(4.32) $$\begin{eqnarray}\unicode[STIX]{x1D719}=\left(\frac{\unicode[STIX]{x1D6FA}[\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D716}+1)-1]}{(\unicode[STIX]{x1D6FE}_{l}\unicode[STIX]{x1D702}_{max}^{l})^{2}}\right)^{1/(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D716}-1)}\left[1-\left(\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{max}^{l}}\right)^{1/\unicode[STIX]{x1D6FE}_{l}}\right]^{\unicode[STIX]{x1D6FA}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}=3/(2(\unicode[STIX]{x1D716}+1-\unicode[STIX]{x1D6FD}))$ .

The pressure profile is obtained by solving (4.30) using finite difference on computational coordinates that are scaled with $\unicode[STIX]{x1D702}_{max}^{l}$ following a procedure similar to that discussed in appendix C. Figure 5(b) shows the solution for $\unicode[STIX]{x1D719}$ for different values of $\unicode[STIX]{x1D6FC}$ . Similar to the short-time regime, $\unicode[STIX]{x1D702}_{max}^{l}$ decreases as $\unicode[STIX]{x1D6FC}$ increases.

Figure 6. Comparing the self-similar pressure profile in the long-time regime for the case where the leakage is simplified, that is when (4.30) is solved, and for the full leakage term, when (B 6) is solved. In both cases, $\unicode[STIX]{x1D6FC}=0$ . The dotted-dashed blue curve corresponds to the case where the leakage is simplified while the solid black curve represents the solution of (B 6). As expected, ignoring the dependence of the leakage velocity on the time a fracture is activated over-predicts the cluster’s propagation rate; i.e. gives a larger value of $\unicode[STIX]{x1D702}_{max}^{l}$ .

To quantify the effects of simplifying the equations by homogenizing the local leakage velocity, the solution of (4.30) was compared to the solution one gets when solving the integro-differential equation presented in appendix B for the case where $\unicode[STIX]{x1D6FC}=0$ . Figure 6 shows the self-similar pressure profile for both cases. As expected, the value of $\unicode[STIX]{x1D702}_{max}^{l}$ for the simple case is larger leading to a higher propagation rate when under-estimating the leakage rate.

Figure 7. The growth of the dimensionless radius of the cluster for different values of $\unicode[STIX]{x1D6FC}$ . The transition from the short-time self-similar solution occurs near the cross-over transitional dimensionless time $\unicode[STIX]{x1D70F}=O(1)$ . The dotted-dashed lines were drawn using the theoretical scalings obtained from the similarity solutions.

4.4 Numerical solution of the partial differential equation

To validate the self-similar solutions and calculate the pressure profile in the transition regime, the simplified partial differential equations given by (4.27) and (4.28) were solved numerically for the two dimensional case. A robust numerical scheme similar to the method described by Zheng, Christov & Stone (Reference Zheng, Christov and Stone2014) was developed for this purpose. This method, which can be used for different values of $\unicode[STIX]{x1D6FC}$ in two and three dimensions, is described in appendix C. It was also used to find the self-similar solution for $\unicode[STIX]{x1D719}$ in the long-time regime by replacing $R$ with $\unicode[STIX]{x1D702}_{max}^{l}$ . Figure 7 shows the non-dimensional cluster radius as a function of time for different values of $\unicode[STIX]{x1D6FC}$ in a log–log plot. In all cases, the radius initially grows as $\unicode[STIX]{x1D70C}\sim \unicode[STIX]{x1D70F}^{\unicode[STIX]{x1D6FE}_{s}}$ where the value of $\unicode[STIX]{x1D6FE}_{s}$ depends on how the fluid is injected with time. After reaching the cross-over time, $\unicode[STIX]{x1D70F}=O(1)$ , the growth rate slows down gradually until the leakage dominates and a new power law behaviour with a different exponent, $\unicode[STIX]{x1D6FE}_{l}$ is established. Figure 8 shows the pressure profile at different times for different values of $\unicode[STIX]{x1D6FC}$ . Figure 8(a) shows the pressure profile for the constant injection rate case while figure 8(b) shows the pressure profile for the linearly ramped up injection rate. The increase in the fluid pressure as $\unicode[STIX]{x1D70C}\rightarrow 0$ is more noticeable when the fluid is injected with a time-dependent rate. This is due to the increase in the pressure drop required to drive the increased flow rate.

Figure 8. The pressure profile inside the cluster of open fractures at different times for the two-dimensional case for two different values of $\unicode[STIX]{x1D6FC}$ . (a) The pressure profile inside the cluster for the case of constant injection rate ( $\unicode[STIX]{x1D6FC}=0$ ). The pressure near the injection point does not change much as the cluster extends outwardly. (b) The pressure profile inside the network for the case of linearly ramping up the injection rate ( $\unicode[STIX]{x1D6FC}=1$ ). The pressure near the injection point increases as the cluster extends outwardly. In both self-similar regimes when $\unicode[STIX]{x1D6FC}=1$ , the pressure at a position $\unicode[STIX]{x1D70C}$ increases with time as $\unicode[STIX]{x1D70F}^{1/(\unicode[STIX]{x1D716}+1)}$ .