Notes on Propagation of 3D Buoyant Fluid-Driven Cracks

It has been observed in laboratory experiments in gelatin (Takada, 1990; Heimpel and Olson, 1994; Taisne and Tait, 2009) that the fracture breadth is nearly stationary when the fracturing process near the fracture head ($z=\ell(t)$), that “generates” the final breadth of the crack $b_{*}$, is dominated by the solid toughness, as opposed to the dissipation in the viscous fluid flow in the crack. This is symptomatic of buoyancy-driven cracks which tend to have maximum opening in the head region, which, in turn, minimizes the viscous losses there.

We turn to study such buoyancy-driven cracks with a small aspect (breadth-to-length) ratio,

$$b_{*}\ll\ell(t)$$

(10)

The determination of the actual value of the stationary breadth $b_{*}(z)$ will require explicit consideration of the structure of the “head region”, which is delayed until the next section. We only note at this point that the breadth is expected to be invariant with depth (fracture height) if the rock properties are depth-independent. Conversely, if the fracture toughness and/or the rock modulus and/or rock density (and thus density mismatch between the fluid and the rock) are depth dependent (which is conceivable if kilometers high dikes/fracture are considered), then the breadth is expected to vary with depth. Thus, in what follows, we allow for the rock modulus and the rock toughness are to be some predetermined functions of depth.

For a finger-like crack, $b_{*}\ll\ell$, and “slowly” varying opening along the crack length, $\partial w/\partial z\sim w/\ell$, the convolution kernel in (2) can be approximated as

$$\frac{1}{[(x-x^{\prime})^{2}+(z-z^{\prime})^{2}]^{3/2}}\approx\frac{2}{(x-x^{\prime})^{2}}\delta_{Dirac}(z-z^{\prime})$$

which then allows to reduce (2) to its plane-strain equivalent (e.g. Bilby and Eshelby, 1968)

$$p(x,z)=-\frac{\bar{E}(z)}{4}\int_{-b_{*}(z)}^{b_{*}(z)}\frac{w(x^{\prime},z)}{(x-x^{\prime})^{2}}dx^{\prime}$$

(11)

Once again, this approximation breaks down in the vicinity of the crack “tips” (i.e., near $z=0$ and $z=\ell$).

Since the pressure is approximately equilibrated along the breadth, (5), $p(x,z)=p(z)$, eq. (11) yields classical elliptical crack opening that can be expressed as

$$w(x,z)=\frac{4\overline{w}(z)}{\pi}\sqrt{1-\frac{x^{2}}{b_{*}^{2}(z)}},$$

(12)

where $\overline{w}(z)$ is the breadth-averaged opening

$$\overline{w}(z)=b_{*}(z)\frac{p(z)}{\bar{E}(z)},$$

(13)

and time $t$ was omitted from the list of arguments in $w$, $\overline{w}$, and $p$.

Integrating the fluid flow equations along the crack breadth (and using the zero lag condition (6)) yields

$$\frac{\partial b_{*}\overline{w}}{\partial t}+\frac{\partial b_{*}\overline{q}}{\partial z}=0,$$

(14)

where $\overline{q}=\overline{wv}$ is the cross-section averaged flow rate

$$\overline{q}=-\frac{\overline{w}^{3}}{\bar{\mu}}\left(\frac{\partial p}{\partial z}-\Delta\rho g\right)$$

(15)

where $\bar{\mu}=\pi^{2}\mu$ is a dynamic viscosity parameter, as defined earlier in (1).

The zero fluid lag condition applied at the “tip” yields

$$\lim_{z\rightarrow\ell}(\overline{q}/\bar{w})=d\ell/dt$$

(16)

while the global fluid volume balance reduces to:

$$V=\int_{0}^{\ell}2b_{*}\overline{w}dz$$

(17)

The closure of the approximate model (12-17), which is in so far equivalent to the classical PKN model of a hydraulic fracture with restricted breadth (Perkins and Kern, 1961; Nordgren, 1972; Kemp, 1990), requires specifying an additional constraint or boundary condition. Naturally, the latter has to do with the crack propagation condition. Unfortunately, since the solution in the head of the fracture is not resolved by the current approximation, the propagation condition in the form (8) can not be applied directly. Instead, the “net” propagation criterion for the fracture head can be established from estimating the global energy release rate in extending the finger crack and equating it to the fracture energy $K_{Ic}^{2}/E^{\prime}$ (Sarvaramini and Garagash, 2015; Garagash, 2022). This “net” propagation constraint is:

$$p=\frac{\bar{K}}{\sqrt{b_{*}}}\quad\text{at}\quad z=\ell$$

(18)

We emphasize again that the approximate model of a finger crack includes the crack breadth $b(z)$. Therefore, when $b(z)$ is known a priori, such as in the case of a hydraulic fracture propagation in a layer sandwiched between tougher/stiffer layers and/or layers subjected to higher compressive stress, this model is complete (Sarvaramini and Garagash, 2015; Garagash, 2022). Otherwise, as is true for a buoyant crack in a homogeneous rock, $b_{*}(z)$ is not known a priori, and is a part of the solution in the fracture head (Section 4).

Let us nondimensionalize the field variables with respect to the scales

	$\displaystyle x_{*}$	$\displaystyle=z_{*}=L_{*}\quad t_{*}=\frac{L_{*}}{v_{*}}\quad w_{*}=\frac{\bar{K}\sqrt{L_{*}}}{\bar{E}}\quad p_{*}=\frac{\bar{K}}{\sqrt{L_{*}}}$
		$\displaystyle v_{*}=\left(\frac{\bar{K}}{\bar{E}}\right)^{2}\frac{p_{*}}{\bar{\mu}}\quad V_{*}=w_{*}L_{*}^{2}\quad Q_{*}=\frac{V_{*}}{t_{*}}$		(19)

which are predicated on a single lengthscale $L_{*}$, chosen here as the buoyancy length (Lister and Kerr, 1991)

$$L_{*}=\left(\frac{\bar{K}}{\Delta\rho g}\right)^{2/3}$$

(20)

Expanded form of the scales (19) with (20) is:

	$\displaystyle x_{*}$	$\displaystyle=z_{*}=\left(\frac{\bar{K}}{\Delta\rho g}\right)^{2/3}\quad t_{*}=\frac{\bar{\mu}\bar{E}^{2}}{\Delta\rho g\bar{K}^{2}}\quad w_{*}=\frac{\bar{K}^{4/3}}{\bar{E}\,(\Delta\rho g)^{1/3}}\quad p_{*}=\left(\Delta\rho g\bar{K}^{2}\right)^{1/3}$
		$\displaystyle v_{*}=\frac{\bar{K}^{8/3}(\Delta\rho g)^{1/3}}{\bar{\mu}\bar{E}^{2}}\quad V_{*}=\frac{\bar{K}^{8/3}}{\bar{E}\,(\Delta\rho g)^{5/3}}\quad Q_{*}=\frac{\bar{K}^{14/3}}{\bar{\mu}\bar{E}^{3}(\Delta\rho g)^{2/3}}$		(21)

Elasticity, fluid continuity and the Poiseuille’s equations can be rewritten in units of (19), or, conversely, in the normalized form as

$$\overline{w}=b_{*}\,p\qquad\frac{\partial\overline{w}}{\partial t}=-\frac{\partial\overline{q}}{\partial z}\qquad\overline{q}=-\overline{w}^{3}\left(\frac{\partial p}{\partial z}-1\right)$$

(22)

Plugging the first and the third into the second in (22) we have

$$\frac{\partial p}{\partial t}=b_{*}^{2}\frac{\partial}{\partial z}\left[p^{3}\left(\frac{\partial p}{\partial z}-1\right)\right]$$

(23)

This non-linear PDE is subjected to the following global and boundary conditions:

$$V=2b_{*}^{2}\int_{0}^{\ell}pdz\quad\text{or}\quad 2b_{*}\bar{q}_{|z=0}=\frac{dV}{dt}$$

(24)

$$(\overline{q}/\bar{w})_{|z=\ell}=\frac{d\ell}{dt}\quad\text{and}\quad p_{|z=\ell}=\frac{1}{\sqrt{b_{*}}}$$

(25)

The normalized solution $p(z,t)$ and $\ell(t)$ of (23-25) depends on the normalized injected fluid volume $V(t)$, and normalized fracture breadth $b_{*}$.

The formulation for the case when the rock properties and the fracture breadth vary with depth is given in Appendix A.

Numerical method of lines is used to solve the non-linear PDE (23) with boundary conditions (24) and (25), as detailed in Appendix B. The numerical solution reveals at large enough times an asymptotic structure comprised of a hydrostatic head region, where viscous stresses are negligible compared to the elastic ones and solution is stationary in the frame moving with the crack tip, and a viscous tail region, where the elastic stress is negligible. These two asymptotic head and tail solutions are discussed below.

The asymptotic solution with the full hydrostatic head attached to the viscous tail is realized when the neck pressure (opening), i.e., as previously defined maximum pressure in the tail occurring at the juncture of the tail and the head, is negligible compared to the prevailing pressure in the head, i.e. $p_{\mathrm{neck}}\ll 1$. This asymptotic framework can be approximately extended to the case when $p_{\mathrm{neck}}$ is not vanishingly small by “attaching” the partial head (“circumcised” at the value of pressure equal to $p_{\mathrm{neck}}$) to the viscous tail, such that the partial head size and volume are

$$\ell_{\mathrm{head}}=\frac{1}{\sqrt{b_{*}}}-p_{\mathrm{neck}}\qquad V_{\mathrm{head}}=b_{*}-b_{*}^{2}p_{\mathrm{neck}}^{2}$$

where $p_{\mathrm{neck}}$ is stationary for a constant injection rate, (29), and decreasing with time for a constant volume release, (30). In the latter case, since $V_{\mathrm{head}}$ evolves with time, so does $V_{\mathrm{tail}}$, which now satisfies an implicit equation

$$V-V_{\mathrm{tail}}=b_{*}-b_{*}^{2}\left(\frac{V_{\mathrm{tail}}}{4b_{*}^{4}t}\right)^{2/3}\qquad(V=\mathrm{const})$$

Corresponding solution with non-stationary volume in the tail is no longer strictly self-similar, yet, we hypothesize that the constant $V_{tail}$ solution formally used with a slowly varying $V_{tail}$ may still yield an adequate approximation…

Here we look for the full solution in the fracture head (which approximate solution we studied in the previous section), including the yet unknown value of the stationary breadth $b_{*}$ away from the tip and the unknown “build-up” $b(z)$ in the near-tip region of a priori unknown vertical extent $\lambda$, where the fracture breadth gradually expands towards the terminal value $b_{*}$ to accommodate the vertical propagation (Figure 4). We further make use of local coordinate $\bar{z}=z-(\ell-\lambda)$ with the origin situated at the boundary $\bar{z}=0$ between the laterally expanding ($\bar{z}>0$, $b(\bar{z})<b_{*}$) and laterally stationary ($\bar{z}<0$, $b(\bar{z})=b_{*}$) parts of the fracture edge.

Using the scales (19), the normalized form of the elasticity equation (2) is

$$p(x,z)=-\frac{1}{8}\int_{\ell_{\mathrm{head}-\lambda}}^{\lambda}d\bar{z}^{\prime}\int_{-b(\bar{z}^{\prime})}^{b(\bar{z}^{\prime})}\frac{w(x^{\prime},\bar{z}^{\prime})}{[(x-x^{\prime})^{2}+(\bar{z}-\bar{z}^{\prime})^{2}]^{3/2}}dx^{\prime}$$

(33)

while the net pressure is hydrostatic in the head

$$p(x,z)=\bar{z}+p_{0}\qquad(x,\bar{z})\in\text{crack footprint}$$

(34)

where $p_{0}$ is yet unknown constant. The normalized conditions along the expanding, (9), and the stationary parts of the fracture edge are, respectively,

$$w_{\perp}(r\rightarrow 0)=\frac{4}{\pi}\sqrt{r}\qquad(0<\bar{z}<\lambda)$$

(35)

$$b(\bar{z})=b_{*}\quad\text{and}\quad w(x\rightarrow b_{*},\bar{z})<\frac{4}{\pi}\sqrt{b_{*}-x}\qquad(\ell_{\mathrm{head}}-\lambda<\bar{z}<0)$$

(36)

where $w_{\perp}(r)$ is, as previously, the opening distribution in the cross-section normal to the fracture edge. Finally, these are complemented by the condition of smooth fracture edge, specifically at the juncture $\bar{z}=0$,

$$db/d\bar{z}=0\quad\text{at}\quad\bar{z}=0$$

(37)

We use the piecewise-constant displacement discontinuity method, as further detailed in Appendix C, to numerically solve elasticity integral equation (33) together with constraints (35-37) for the terminal fracture breadth and the reference pressure value, respectively,

$$b_{*}\approx 0.396,\qquad p_{0}\approx 1.344,$$

(38)

the crack opening $w(x,\bar{z})$ show on Figure 5, and the fracture breadth distribution

$$b(\bar{z}\geq 0)\approx b_{*}\sqrt{1-A(\bar{z}/\lambda)^{2}-(1-A)(\bar{z}/\lambda)^{4}}\quad\text{with}\quad A\approx 0.6967$$

in the laterally-expanding part (of the head region) of extent

$$\lambda\approx 0.552$$

The fracture breadth is stationary, $b(\bar{z}<0)=b_{*}$, in the remainder of the head region of extent

$$\ell_{\mathrm{head}}-\lambda\approx 1.504$$

The true extent of the head region is therefore

$$\ell_{\mathrm{head}}\approx 2.055$$

(39)

We note that the net pressure at the tail-end of the head, i.e. at $\bar{z}=-(\ell_{\mathrm{head}}-\lambda)$, has a negative value $\approx-0.1$, which owes to the non-local 3D crack effects there. (Departure from the PKN assumption, which would have required zero net-pressure there).

The volume of the head is

$$V_{\text{head}}\approx 0.408$$

(40)

The PKN-approximation of the solution in the head, as shown on on Figure 5b by dashed line, corresponds to $\overline{w}^{\mathrm{PKN}}=b_{*}p$, where $p$ in the head is taken from the full head solution (34) with (38). The PKN approximation of the advancing front location $\bar{z}_{\mathrm{front}}^{\mathrm{PKN}}\approx 0.245$ follows from the propagation condition $p(\bar{z}_{\mathrm{front}}^{\mathrm{PKN}})=1/\sqrt{b_{*}}$, while the back (reseeding) front $\bar{z}_{\mathrm{back}}^{\mathrm{PKN}}\approx-1.344$ corresponds to the closure condition $p(\bar{z}_{\mathrm{back}}^{\mathrm{PKN}})=0$. Although the validity of the PKN assumption is at best questionable in the head region, which extent $\ell_{\mathrm{head}}$ is only about 2.5 times larger than its terminal breadth $2b_{*}$, (39), the PKN solution for the breadth-averaged opening starts to track the full numerical solution some distance behind the tip (Figure 5b). What is however even more remarkable is that the equivalent PKN volume of the head region $V_{\mathrm{head}}^{\mathrm{PKN}}=b_{*}\approx 0.396$ (see (28) and (38)), is in very good agreement ($\sim 3$% different) with the full numerical solution (40). This lands support to the PKN solution for the finger crack (including the hydrostatic head, viscous tail, and the transition between the two) developed in Section 3, as long as the head region is at least partially developed, see constraints (31) for the constant injection rate and (32) for the constant volume release cases.

Draw comparisons to Weertman’s 2D pulse solution (Weertman, 1971)….

$$w=\frac{2\Delta\rho g}{E^{\prime}}L_{c}^{2}\left(1+\frac{Z}{L_{c}}\right)^{3/2}\left(1-\frac{Z}{L_{c}}\right)^{1/2},\qquad p=\Delta\rho g\left(Z+L_{c}/2\right)$$

where

$$L_{c}=\frac{L_{pulse}}{2}=\left(\frac{K_{Ic}}{\sqrt{\pi}\Delta\rho g}\right)^{2/3}$$

is the half-length of the pulse ($Z=\pm L_{c}$ correspond to advancing/reseeding tips in the coordinate system used in the Weertman’s solution). In our notation, $L_{c}=L_{*}/2^{1/3}$ where $L_{*}=(\bar{K}/\Delta\rho g)^{2/3}$ is a buoyancy length, $\bar{K}=\sqrt{2/\pi}K_{Ic}$ and $\bar{E}=E^{\prime}/\pi$.

Using scaling (19), i.e. scales $L_{*}$, $w_{*}=\bar{K}\sqrt{L_{*}}/\bar{E}$, and $p_{*}=\bar{K}/\sqrt{L_{*}}$, we have

$$w=\frac{2^{1/3}}{\pi}\left(1+2^{1/3}Z\right)^{3/2}\left(1-2^{1/3}Z\right)^{1/2},\qquad p=Z+2^{-4/3}$$

and non-dimensional pulse length is $2\times 2^{-1/3}=2^{2/3}$.

In order compare the 2D and 3D pulse solution, we need to relate the corresponding coordinate frames, i.e., $Z$ and $\bar{z}$, respectively. We choose to do so by requiring that the net pressure distributions in the two solutions are identical, i.e., $Z+2^{-4/3}=\bar{z}+p_{0}$, which leads to $Z=\bar{z}+\left(p_{0}-2^{-4/3}\right)\approx\bar{z}+0.947$.

The normalized volume of the Weertman’s 2-D pulse is equal to $0.5$, which approximates well the equivalent 2-D-volume of the 3-D pulse, $V_{\text{head }2D}=V_{\text{head }}/2b_{*}\approx 0.515$. This result, taken together with the approximate correspondence of the 2-D and 3-D viscous tail solutions obtained earlier, suggests that the 2-D framework provides a good approximation (within 10% in terms of the dike length, width, and head vs. tail volume partition) of the 3-D dike solution as long as the appropriate (obtained from the 3-D hydrostatic head solution) value of the stationary breadth $b_{*}$ is used to convert the dike’s volume to its 2-D proxy.