Dynamic Development of Hydrofracture

Abstract

Many natural examples of complex joint and vein networks in layered sedimentary rocks are hydrofractures that form by a combination of pore fluid overpressure and tectonic stresses. In this paper, a two-dimensional hybrid hydro-mechanical formulation is proposed to model the dynamic development of natural hydrofractures. The numerical scheme combines a discrete element model (DEM) framework that represents a porous solid medium with a supplementary Darcy based pore-pressure diffusion as continuum description for the fluid. This combination yields a porosity controlled coupling between an evolving fracture network and the associated hydraulic field. The model is tested on some basic cases of hydro-driven fracturing commonly found in nature, e.g., fracturing due to local fluid overpressure in rocks subjected to hydrostatic and nonhydrostatic tectonic loadings. In our models we find that seepage forces created by hydraulic pressure gradients together with poroelastic feedback upon discrete fracturing play a significant role in subsurface rock deformation. These forces manipulate the growth and geometry of hydrofractures in addition to tectonic stresses and the mechanical properties of the porous rocks. Our results show characteristic failure patterns that reflect different tectonic and lithological conditions and are qualitatively consistent with existing analogue and numerical studies as well as field observations. The applied scheme is numerically efficient, can be applied at various scales and is computational cost effective with the least involvement of sophisticated mathematical computation of hydrodynamic flow between the solid grains.

Appendices

Appendices

1.1 Appendix-A: ADI—2D Pressure Diffusion

The ADI method is time implicit. With symmetric discretization in time i.e., between a forward and backward step, this methods is unconditionally stable and the precision is better than with a purely forward in time implicit method [Press, 1992]. The two-dimensional pressure diffusion Eq. (9) can be rewritten as

$$ \frac{{\partial P(\overrightarrow {r} ,t)}}{\partial t} = \left( {1 + \beta P} \right)\frac{{K(\overrightarrow {r} ,t)}}{{\beta \mu \phi (\overrightarrow {r} ,t)}}\left[ {\frac{{\partial^{2} P(\overrightarrow {r} ,t)}}{{\partial x^{2} }} + \frac{{\partial^{2} P(\overrightarrow {r} ,t)}}{{\partial y^{2} }}} \right] - \frac{1}{{\beta \phi (\overrightarrow {r} ,t)}}g(\overrightarrow {r} ,t) $$

(A1)

where $g(\overrightarrow {r} ,t)$ is the source term and $\overrightarrow {r}$ stands for position in space. This is a second-order parabolic partial differential equation. Corresponding to the time and space discretization of the 2D pressure continuum using forward difference with time on the left hand side and central difference with space on the right hand side of equation (A1).

$$ \frac{{P_{i,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} - P_{i,j}^{n} }}{\Updelta t} = (1 + \beta P)\frac{{k_{i,j} }}{{\beta \mu \phi_{i,j} }}\left[ {\frac{{P_{i + 1,j} - 2P_{i,j} + P_{i - 1,j} }}{{(\Updelta x)^{2} }} + \frac{{P_{i,j + 1} - 2P_{i,j} + P_{i,j - 1} }}{{(\Updelta y)^{2} }}} \right] - \frac{1}{\phi \beta }g_{i,j} $$

(A2)

where suffixes $i,\,j$ and $n$ are the indices in the x, y, and t directions respectively.

The main idea of the ADI method is to reduce the 2-D problem into a succession of two one-dimensional problems by proceeding one time step from $n$ to $n + 1$ in two sub-time steps (Fig. 11). The first half-step ($n$ to $n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}$) is taken implicitly in the x-direction and explicitly in the y-direction followed by the second half-step ($n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}$ to $n + 1$) that is taken implicitly in the y-direction and explicitly in the x-direction.

Fig. 11

Schematic diagram of ADI solution of finite-difference pressure continuum, after Wang and Chen (2001)

Detailed differential equations in stage-1 for each $j$ at marched time $n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}$ and the corresponding tridiagonal system of equations for the respective one-dimensional problem can be derived in form of matrix equation of dimension I:

$$ - \alpha_{i,j} P_{i + 1,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} + (1 + 2\alpha_{i,j} )P_{i,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} - \alpha_{i,j} P_{i - 1,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} = \gamma_{i,j} P_{i,j + 1}^{n} + (1 - 2\gamma_{i,j} )P_{i,j}^{n} + \gamma_{i,j} P_{i,j - 1}^{n} - \frac{\Updelta t}{2\phi \beta }g_{i,j} $$

(A3)

$$ \begin{gathered} \left[ {\begin{array}{*{20}c} 1 & 0 & 0 & \ldots & \ldots & 0 \\ { - \alpha_{i,j} } & {1 + 2\alpha_{i,j} } & { - \alpha_{i,j} } & 0 & \ldots & \ldots \\ 0 & \ldots & \ldots & \ldots & \ldots & 0 \\ \ldots & \ldots & 0 & { - \alpha_{i,j} } & {1 + 2\alpha_{i,j} } & { - \alpha_{i,j} } \\ 0 & \ldots & \ldots & 0 & 0 & 1 \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {P_{0,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} } \\ {P_{1,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} } \\ \ldots \\ \ldots \\ {P_{I,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} } \\ \end{array} } \right] = \hfill \\ \left[ {\begin{array}{*{20}c} 1 & 0 & 0 & \ldots & \ldots & 0 \\ {\gamma_{i,j} } & {1 - 2\gamma_{i,j} } & {\gamma_{i,j} } & 0 & \ldots & \ldots \\ 0 & \ldots & \ldots & \ldots & \ldots & 0 \\ \ldots & \ldots & 0 & {\gamma_{i,j} } & {1 - 2\gamma_{1,j} } & {\gamma_{i,j} } \\ 0 & \ldots & \ldots & 0 & 0 & 1 \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {P_{i,0}^{n} } \\ {P_{i,1}^{n} } \\ \ldots \\ \ldots \\ {P_{i,J}^{n} } \\ \end{array} } \right] - \frac{\Updelta t}{2\phi \beta }g_{i,j} \hfill \\ \end{gathered} $$

(A4)

$$ i = 0,1, \ldots ,I; \, \quad j = 0,1, \ldots ,J $$

where

$$ \alpha_{i,j} = (1 + \beta P)\frac{{K_{i,j} \Updelta t}}{{2\mu \beta \phi_{i,j} (\Updelta x)^{2} }}\;\;\;\;{\text{and }}\gamma_{i,j} = (1 + \beta P)\frac{{K_{i,j} \Updelta t}}{{2\mu \beta \phi_{i,j} (\Updelta y)^{2} }} $$

By analogy, stage-II of the ADI method for each $i$ at time $n + 1$, is expressed in tridiagonal system of dimension J:

$$ - \gamma_{i,j} P_{i,j + 1}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} + (1 + 2\gamma_{i,j} )P_{i,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} - \gamma_{i,j} P_{i,j - 1}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} = \alpha_{i,j} P_{i + 1,j}^{n} + (1 - 2\alpha_{i,j} )P_{i,j}^{n} + \alpha_{i,j} P_{i - 1,j}^{n} - \frac{\Updelta t}{2\phi \beta }g_{i,j} $$

(A5)

$$ \begin{gathered} \left[ {\begin{array}{*{20}c} 1 & 0 & 0 & \ldots & \ldots & 0 \\ { - \gamma_{i,j} } & {1 + 2\gamma_{i,j} } & { - \gamma_{i,j} } & 0 & \ldots & \ldots \\ 0 & \ldots & \ldots & \ldots & \ldots & 0 \\ \ldots & \ldots & 0 & { - \gamma_{i,j} } & {1 + 2\gamma_{i,j} } & { - \gamma_{i,j} } \\ 0 & \ldots & \ldots & 0 & 0 & 1 \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {P_{i,0}^{n + 1} } \\ {P_{i,1}^{n + 1} } \\ \ldots \\ \ldots \\ {P_{i,J}^{n + 1} } \\ \end{array} } \right] = \hfill \\ \left[ {\begin{array}{*{20}c} 1 & 0 & 0 & \ldots & \ldots & 0 \\ {\alpha_{i,j} } & {1 - 2\alpha_{i,j} } & {\alpha_{i,j} } & 0 & \ldots & \ldots \\ 0 & \ldots & \ldots & \ldots & \ldots & 0 \\ \ldots & \ldots & 0 & {\alpha_{i,j} } & {1 - 2\alpha_{i,j} } & {\alpha_{i,j} } \\ 0 & \ldots & \ldots & 0 & 0 & 1 \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {P_{0,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} } \\ {P_{1,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} } \\ \ldots \\ \ldots \\ {P_{I,j}^{{n + {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} } \\ \end{array} } \right] - \frac{\Updelta t}{2\phi \beta }g_{i,j} \hfill \\ \end{gathered} $$

(A6)

$$ i = 0,1, \ldots ,I; \, j = 0,1, \ldots ,J $$

Implementing the Gauss-algorithm with a Dirichlet boundary condition, the derived tridiagonal system of Eq. (A4) is solved $J$ times and Eq. (A6) by $I$ times.

Appendix B

See Figs. 11 and 12

Fig. 12

Different patterns of hydrofractures in a situation where a constant point source is injected in a homogeneous media under different remote stresses: a with relatively smaller $\sigma_{y}$, the results show the initial development of circular fracturing at the source location (centre) with elongated fractures oriented parallel to the axis of the applied stress. In contrast to the pattern shown in a, in b a larger $\sigma_{y}$ dominates the overall pattern and results only in sub-vertical oriented fractures parallel to the main stress axis and through the source location (central). c The figure shows the state of stress field at the onset of fracturing, where the red color code represents high differential stress and blue low differential stress

Appendix C

See Fig. 13

Fig. 13

This figure shows a time series illustrating the influence of the local pore overpressure on brittle failure in a pure shear stress regime. In this case the extension fractures develop at the source location (centre) and link up with shear fractures towards the edges of the simulation box (time steps T increase from left to right)