Abstract
At the time of energy transition, it is important to be able to predict the effects of fluid overpressures in different geological scenarios as these can lead to the development of hydrofractures and dilating high-porosity zones. In order to develop an understanding of the complexity of the resulting effective stress fields, fracture and failure patterns, and potential fluid drainage, we study the process with a dynamic hydromechanical numerical model. The model simulates the evolution of fluid pressure buildup, fracturing, and the dynamic interaction between solid and fluid. Three different scenarios are explored: fluid pressure buildup in a sedimentary basin, in a vertical zone, and in a horizontal layer that may be partly offset by a fault. Our results show that the geometry of the area where fluid pressure is successively increased has a first-order control on the developing pattern of porosity changes, on fracturing, and on the absolute fluid pressures that sustained without failure. If the fluid overpressure develops in the whole model, the effective differential and mean stress approach zero and the vertical and horizontal effective principal stresses flip in orientation. The resulting fractures develop under high lithostatic fluid overpressure and are aligned semihorizontally, and consequently, a hydraulic breccia forms. If the area of high fluid pressure buildup is confined in a vertical zone, the effective mean stress decreases while the differential stress remains almost constant and failure takes place in extensional and shear modes at a much lower fluid overpressure. A horizontal fluid pressurized layer that is offset shows a complex system of effective stress evolution with the layer fracturing initially at the location of the offset followed by hydraulic breccia development within the layer. All simulations show a phase transition in the porosity where an initially random porosity reduces its symmetry and forms a static porosity wave with an internal dilating zone and the presence of dynamic porosity channels within this zone. Our results show that patterns of fractures, hence fluid release, that form due to high fluid overpressures can only be successfully predicted if the geometry of the geological system is known, including the fluid overpressure source and the position of seals and faults that offset source layers and seals.
1. Introduction
Fracturing and the development of mineralized veins and breccia are an important process in the Earth’s crust linked to fluid flow and mineral deposits and have important applications related to the energy transition with reactive flow in geothermal systems and carbon capture and storage (CCS) in aquifers and decommissioned oil and gas fields and energy storage in sedimentary basins [1–6]. On the one hand, fluids can flow along open fractures and transport material that can precipitate and close fractures forming mineral veins (Figure 1). On the other hand, the fluid pressures are directly involved in the fracturing process itself through the mechanism of hydrofracturing. This process is highly dynamic with feedback between fluid overpressure (fluid pressure in excess of the hydrostatic pressure) and the solid material involving porosity and permeability changes due to the fracturing process and subsequent opening and closing of fractures [7, 8]. A rise in the local fluid pressure in a rock pore leads to an outward directed force that works against the solid and can thus reduce the solid stress to an effective stress according to Terzaghi’s law [9]: with being the total normal stress, the pore fluid pressure, and the Kronecker delta with the sign convention of positive compressive stress.
(a)
(b)
(c)
(d)
This effect is often visualized in the Mohr diagram of effective shear to normal stress where a rise in fluid pressure leads to a reduction of the mean stress and thus eventually to failure if the fluid pressure is high enough. This approach is based on assumptions and is a simplification of the effects of fluid pressures on solids and has a number of shortcomings that can lead to misinterpretation as pointed out by a number of authors ([10–14]; Cobbold and Rodriguez, 2007; [15]). The three main assumptions are as follows: (a) fluids are incompressible, (b) fluids are stationary, and (c) it is sufficient to consider fluid pressure only, disregarding fluid overpressure and fluid pressure gradient(s). However, in a natural geological scenario, these assumptions may not be appropriate; thus, conclusions based on these assumptions may be faulty and misleading. For example, fluid is generally about an order of magnitude more compressible than the surrounding solid. This compressibility is not necessarily important for the initial fracturing process but becomes important once the system evolves and the solid and fluid interact dynamically leading to opening and closing channels, for example, in the fluidization of sediments [16–19]. The assumption that the fluids are stationary is not realistic for a dynamic natural system. In a porous or fractured system, fluids will generally move leading to an evolution of the fluid pressure and fluid pressure gradients depending on the source, sink, and permeability of the system as was shown by Cobbold and Rodrigues [20]. This does not necessarily mean that the fluids move fast but that changes in pressure or a local fluid pressure increase from a source will eventually lead to evolving pressure gradients. Commonly, this shortcoming of the stationary fluid assumption is overcome by using Biot’s law where a linear pore fluid factor is added that represents the loss of fluid from a pore into the surrounding rock [21, 22] according to with the pore fluid factor representing the ratio of pore fluid pressure () to vertical stress (), the poroelastic factor, and the Poisson ratio.
Using this law in three dimensions assumes that all principal stresses are affected by the fluid pressure in the same way, which is a simplification that only works in restricted cases [10, 11, 13, 14]. The use of fluid pressure only in a dynamically evolving setting such as a sedimentary basin where through some dynamic processes, e.g., mineral reaction or oil maturation, fluid pressure increases successively and the fracture system may evolve is inappropriate. Here, only a lateral difference in fluid pressure or a fluid pressure gradient can lead to forces that act on the solid and eventually lead to fracturing associated with fluid pressure. Therefore, only an overpressure or a pressure gradient can be used in the effective stress law as was pointed out by Cobbold and Rodrigues [20]. An example is the evolution of effective stress in a fluid overpressure zone in a sedimentary basin below a seal, where Hillis [12] has shown that in this case, the lowest effective stress below the seal is linked to the vertical stress and the fluid pressure. That means that the effective stresses in the horizontal and vertical directions will approach zero, the differential and mean stress will decrease, and the system will fluidize. By fluidize, we mean the point where the system is only governed by a pressure without differential, shear, or deviatoric stresses. The reason for this effective stress pattern is that the fluid pressure increase is nearly constant for a given depth across the basin so that no horizontal fluid pressure gradients will develop. Gradients develop only in the vertical direction, and in order for failure to occur, the fluid overpressure will have to exceed the vertical stress. Cobbold and Rodrigues [20] show that the horizontal effective stress below a seal in an overpressurized sedimentary basin is where is the local overpressure and the dimensionless elastic proportionality factor (smaller than 1.0). The fluid overpressure affects the vertical and horizontal stresses differently, leading to a reduction of the differential stress and to horizontal failure so that “beef” veins develop (Figure 1). This scenario was illustrated by Cobbold and Rodrigues [20] in experiments. First, the system fluidizes, the effective stresses approach zero, and then a horizontal hydrofracture develops once the fluid pressure approaches the overburden stress (Figure 2). Ghani et al. [15, 23] use a numerical model to illustrate the difference between tectonic fractures and pure hydrofractures and derive similar results as Cobbold and Rodrigues [20] in their experiments.
(a)
(b)
(c)
In this contribution, we use an advanced version of the numerical model of Ghani et al. [15, 23] that can model the dynamic interplay between influx of compressive, nonstationary fluids, overpressure, fluid pressure gradients, associated fracturing, and dynamic feedback between these factors; hence, with this model, it is possible to simulate the dynamic evolution of a geological scenario where fluid pressure may increase locally or homogeneously. We study the effects of different layer geometries and elasticity on the evolution of the effective stress, fracturing, and dewatering of the overpressurized systems.
2. Numerical Setup
2.1. General Model
The numerical scheme used in this contribution is built on the work of Ghani et al. [15, 23] with a similar approach as that used for poroelasticity and fluid flow in granular media [24–26], during hydrofracture or aerofracture [27–31], for shearing of gouge layers in faults [32, 33], and for instabilities during sedimentation [16–19]. The model is implemented into the microstructural modeling environment “Elle” [34, 35].
In the 2-dimensional numerical scheme, the fluid and solid are treated on two different grids, with the solid being represented by an initially triangular elastic spring network and the fluid being represented by a square continuum grid (Figure 2(a)). The model represents a vertical square cross section of a crustal domain that is overlain by 1000 m of cover sediments. The upper model boundary is controlled by gravity (force boundary) whereas the two sides and the bottom are free-slip parallel to a wall but not allowed to move perpendicular to the wall (Figure 2(a)). The fluid pressure is horizontally periodic, wrapping on the left- and right-hand side of the model and has a constant value at the upper and lower boundary. Initial fluid pressure conditions are hydrostatic with fluid being inserted into the model over time in different geometrical configurations (Figure 2). We first present the governing equations of the fluid, followed by the solid and the connection of model assumptions and the scenarios modelled.
2.2. Governing Equations
The fluid phase is described by the fluid pressure of single nodes in a square grid. The inertia of the fluid is not considered assuming that the Reynolds number is low. Darcy’s law is used to describe the fluid movement through the solid. A diffusion equation is derived for the fluid pressure that contains mass and momentum conservation between the fluid and solid. The interstitial fluid flow is then expressed as a porosity-dependent pressure gradient. The continuity equations for solid and fluid at the scale of a grain diameter read [15] where and are the solid and fluid densities, respectively; and are the solid and fluid velocities, respectively, and is the local porosity of the solid. The Darcy equation can be used to calculate a local seepage velocity of the fluid for a pressure change according to the local permeability on a unit area: where is the fluid viscosity and P the fluid pressure. The permeability is calculated from the local porosity according to the Kozeny-Carman relation [36]: with being the grain radius. The fluid state equation is considered using the fluid compressibility , as a proportional approximation of the fluid density to pressure variation: with being the fluid density at some reference pressure. When and are substituted into equation (5), is eliminated and the following diffusion equation for the fluid overpressure is derived [15, 32, 33, 37]: with the left-hand side of the equation representing the Lagrangian derivative of the fluid pressure, the first term on the right-hand side the Darcy diffusion of the fluid pressure, and the third term the source term dealing with pressure change as a function of particle movement.
2.3. Solid
The solid is represented by an initially triangular elastic spring network with the translation of nodes as a function of the momentum exchange between the solid and fluid in a unit volume cell , with a solid of mass , and interparticle force , fluid force , and gravitational loading , so that
The normal () and shear force () acting on a particle from its neighbor is calculated from the relative displacements of the neighbor () normal () and tangential () to the connecting spring with reference to an equilibrium position. The total force over all connected springs is [38] where and are the spring constants for normal and shear displacement, respectively, and the sums are running over all six neighboring particles. Once springs break in the model, they are removed and the corresponding nodes will only experience a repulsive force. The repulsive force is calculated with the normal force from equation (11) but only for compressive interactions. The fluid forces that act on each node in the elastic grid are calculated from the difference of the neighboring fluid pressure cells and then applied on the area that the elastic node represents. where is the fluid pressure gradient, is the solid fraction, is running over the four fluid grid nodes near the particle, and is a smoothing function that satisfies the weighted distribution of particle mass relative to its position [15]. A gravitational vertical force is applied on the elastic nodes depending on the depth of the upper boundary (representing the overlying rocks) and gravitational forces from neighboring nodes. The gravity force on single particles in addition to the load of the overlying sediments is calculated from the particle density , the acceleration due to gravity , the real volume of the particle , and the scale factor (, [38]):
Failure in the elastic network will happen once a critical strain energy for shear and tensile failure is reached. In order to model combinations between the two failure modes, an elliptical energy model for mixed failure is used [39]. The strain energy () in the elastic network is where and are the strain energies for tension and shear, respectively. With the critical strain energy for failure being and for tension and shear, respectively, then where is the normal stress, is the shear stress, is the tensile strength, and is the shear cohesion of the material. The equation typically describes an ellipse in the space [40].
2.4. Interaction between Fluid and Solid
In the current configuration, the square fluid lattice cells are twice as large as the diameter of the elastic nodes of the triangular grid. This configuration is important in order to ensure that the solid fraction per fluid cell is calculated accurately to determine the larger-scale porosity and permeability. In order to communicate between the two lattices, a linear tent weight function is used to account for the differences in distance between the elastic and fluid nodes [15, 41]. The solid lattice passes the local porosity (or particle solid fraction) and particle velocity to the fluid, and the fluid passes the local fluid force back to the solid. The particle density or solid fraction in a fluid cell is calculated as and the overall solid velocity in a fluid cell is where subscript stands for the particle number and runs through all particles and the smoothing function satisfies the weighted distribution of particle mass relative to the fluid node position [15]. The simulation input and output data used to support the findings of this study are available from the corresponding author upon request. The basic software for the simulations can be found and downloaded at http://elle.ws, and the corresponding author will make the additional code available upon request.
2.4.1. Assumptions
The friction force between the fluid and the solid at the surface of the solid at a macroscopic scale is not considered, so that only the pressure gradient produces a drag force on the solid nodes (in the direction of the flow). This pressure gradient is nonetheless directly related to the fluid friction forces exerted at a small scale from the fluid over the solid along the pore boundaries, since the force associated with the pressure gradient exerted over the fluid is entirely transmitted to the pore boundaries in a Darcian flow with a negligible fluid inertia. The fluid is considered to be purely viscous; hence, a thermal evolution is not taken into account. We assume that the Kozeny-Carman relation that is derived from dense granular media can be applied to relate porosity and permeability in our setting. The local grain radius in the Kozeny-Carman equation can be used to adjust the relationship; however, in natural rocks, the relationship between porosity and permeability may be more complicated. In the model, we use a fixed grid geometry and a fixed node size for the solid. This configuration was chosen because it is the only configuration in 2d that produces linear elasticity on the large scale [24]. The fixed configuration has the potential disadvantage of grid effects. However, a strong variation of grid geometry and particle size in DEM models produces nonreliant material behaviors, which we aim to avoid. In order to minimize grid effects, a Gaussian variation of the spatial distribution of the breaking strengths is applied to the springs. We expect grid effects to become largest under shear fracturing especially when the grid weakness is oriented close to the actual fracture orientation. However, the model does still produce listric and curving faults [39, 41]. Fractures are not wrapping in the model, which leads to boundary effects on the left- and right-hand side and the upper and lower boundary. This results in high fracture densities close to the model boundaries effecting up to 10% of the simulation.
We also assume that the Darcy flow is valid in the simulations because the Reynolds number of the flow is small. In the simulations, the velocity of the fluid is the Darcy velocity divided by the porosity:
The average velocity of the fluid in the simulations is 10-5–10-7 m/s. The Reynolds number can be calculated as with being the typical grain or fracture channel size considered and the dynamic water viscosity. With a fluid density of 1000 kg/m3, a channel size of 10-5 m, and a dynamic water viscosity of 0.001 Pa s, Re is in the range of 10-3 to 10-5. Thus, the Darcy flow is still valid since .
2.5. Modelled Scenarios: Setup, Parameters, and Geological Relevance
The simulations start with an initial geometrical setup and the application of the gravitational force and hydrostatic fluid pressure. During each model time step, a change in fluid pressure in a number of cells is conducted. At each time step, the fluid pressure evolution is calculated as a function of the local solid fraction and solid node velocities, the pressure diffusion is applied, the elastic network is relaxed (finding a new equilibrium configuration), and failure is applied if the local critical strain energy exceeded in successive steps with intermediate relaxation until all bonds are stable (Figure 2(b)). The resolution of the simulations is grid cells for the elastic solid and grid cells for the fluid (Figure 2). The simulations run between 2 hours and 4 weeks due to the strong nonlinearity of the process. All simulations have the following settings: the size of the box is (1.0 in model dimension), the overburden is 1000 m, the overburden density is 2400 kg/m3, the Poisson ratio is 0.33, the internal angle of friction of the solid is set to 30 degrees [42], the mean breaking strength is 7 MPa with a Gaussian variation between a lower bound of 2 MPa and an upper bound of 20 MPa, the porosity is 0.1, the fluid density is 1000 kg/m3, the fluid viscosity is 0.001 Pa·s, the fluid compressibility is , and the Carman-Kozeny grain size is 10 μm. In all models, the injected fluid pressure input is given to 30 random nodes within the designated source zones (cf. Figure 2(c)) per time step, and the time step represents 1 hour. That means that the fluid pressure in the model is increased by 0.004 MPa/hour in the whole box. This fluid input is adjusted to a normalized fluid input per fluid cell in cases where the high-pressure zone is smaller than the whole box, for example, when fluid only enters the horizontal layer. The gravity boundary condition is stress driven and controlled by the weight of the particles. Stresses in the model are calculated using average stress tensors for model realms [24, 43, 44], making the stresses independent of grid sizes. Mean stresses for the simulations are calculated for the innermost 80% of the simulation box to avoid boundary effects.
We model three main scenarios. Scenario 1 represents a sedimentary basin where the fluid pressure is increased below a seal, for example, due to a diagenetic reaction or hydrocarbon maturation. Scenario 2 represents fluid pressure increase in a laterally confined zone. This could represent, for example, a sedimentary basin where the fluid enters from below into a confined zone or where fluid pressure is created due to a diagenetic reaction in a confined zone or cell. This scenario is also similar to a triaxial compression experiment where a sample under stress is injected by fluid. Scenario 3 represents a horizontal layer where fluid pressure increases, again for example, due to a diagenetic reaction. For scenario 3, we also present cases where the layer is offset (for example, by a fault) prior to fluid pressure increase (note that the offset zone or fault has no specific properties).
3. Results
3.1. Scenario 1 versus Scenario 2
Scenario 1 represents a simulation where the fluid pressure is raised in the whole simulation box mimicking an endless system in the horizontal layer, for example, a high fluid pressure region in a sedimentary basin below a seal (Figure 2). Figure 3(a) shows the stress patterns (average over the model) over time for scenario 1 (Figure 2(c), Table 1) and illustrates that in this geometry, the mean and the differential effective stress in addition to the vertical and horizontal effective stresses decrease until they all become zero. Hence, at this point in time, the system fluidizes (shear stresses vanish) before the horizontal and vertical stress switch so that the vertical stress is the smallest effective stress and the differential stress increases again. This then leads to a horizontal zone of high porosity as can be seen in Figures 3(b) and 3(c) and a horizontally aligned fracture zone in Figure 3(d). The high-porosity domain is in the lower region of the high fluid pressure zone (Figure 3(b)), and its geometry is relatively rough and wavy. This wavy character of the high-porosity domain is reflected by the fracture pattern (Figure 3(d)) that shows a breccia-like intergrowth of fractures. The dense fracture pattern at the left- and right-hand boundary represents an artifact and is an effect of the nonwrapping nature of the boundary.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Scenario 2 represents a simulation where the increase in fluid pressure is concentrated in a vertical zone in the center of the model with 10% of the model on the left- and right-hand side having a normal hydrostatic pressure and the resulting effective solid stress of . As an experiment, this scenario can be compared to a cylindrical rock deformation experiment with an increase in the fluid pressure under stress. Figure 3(e) shows the resulting effective stresses (average over the model) in the and direction of the scenario 2 model as well as the mean and differential stress. The differential stress in the simulation stays almost constant whereas the mean stress becomes almost zero, the effective stress in the vertical orientation decreases, and the stress in the horizontal orientation becomes tensile. In the Mohr circle diagram, this scenario would be illustrated by the circle moving towards the left-hand side without significantly changing its radius leading to tensile or hybrid tensile/shear failure depending on the shape of the failure curve. Figures 3(f), 3(g), and 3(h) show the fluid pressure, porosity, and mean stress with fractures towards the end of the simulation. The fluid pressure is high in the central part of the simulation and decreases towards the right- and left-hand side as well as towards the top. The porosity shows a localized increase in the central region of the model with fracture-like channels of very high porosity. The high-porosity zone in the center has a much smaller width than the high fluid pressure zone (25-30% of the model width compared to 80%) and is surrounded by a compacted low-porosity zone. Fractures develop in the whole zone where the fluid pressure is high. Two distinct fracture patterns can be distinguished: conjugate shear fractures that are dominating the left- and right-hand side of the high fluid pressure zone and tree-like horizontal to vertical fracture branches that form a breccia-like pattern and develop in the central part of the model. The dense fracture pattern at the lower and upper boundary of the model represents an artifact and is an effect of the nonwrapping nature of the boundary.
Figure 4 shows the evolution of the fracture pattern for both scenarios over time with the fluid pressure in the background. The fracture evolution in the first scenario starts roughly an order of magnitude later than that in scenario 2 (Figure 4(a)). Here, an initial fracture develops horizontally across the model in a very rough and wavy fashion in the middle of the high fluid pressure area (red zone in Figure 4(a)). Successive fractures develop progressively below the initial fracture and merge to form a breccia zone. We use the term breccia for any piece of host rock that has an angular shape and is completely surrounded by fractures. In scenario 2 (Figure 4(b)), the fracture patterns start an order of magnitude earlier than that in scenario 1 with the first fractures developing as mode I or mixed mode I and II fractures at the rim of the high fluid pressure zone on the left- and right-hand side of the model. They then form conjugate shear fracture sets and propagate from both sides towards the center of the high-pressure zone. Once the fractures have reached the center of the high-pressure zone, a second set of fractures develops within the center and forms mode I and tree-like fractures that merge to form a breccia. The porosity evolution of the two scenarios is illustrated in Figure 5.
(a)
(b)
(a)
(b)
The initial porosity for both scenarios is random and becomes bimodal during successive fracturing and channel opening. In scenario 1, the porosity is initially random and localizes once the horizontally aligned hydrofracture develops. The fracture opens and becomes a channel of high porosity whereas the surrounding host rock is compressed with relatively low porosity. The high-porosity channel then moves downwards in the simulation following successive fracture development towards the lower part of the simulation whereas the solid is pushed upwards by the fluid overpressure. In scenario 2, the porosity increases inwards towards the center of the high-pressure zone and decreases outwards towards the boundaries. At the same time, when the fracture pattern switches from the conjugate shear to the more tree- and breccia-like fractures, the high porosity localizes strongly in the center of the model and two compacted areas develop next to the high-porosity zone (Figure 5(b)). In terms of solid movement, the solid is pressed towards the right- and left-hand side of the model leading to an opposite directed increase in porosity or the creation of a stationary porosity zone or wave in the center.
3.2. Scenario 3: Horizontal Layer, Faulting, and Variation of Young’s Modulus
In the third scenario, the fluid pressure is increased in a horizontal layer. This scenario can mimic hydrocarbon generation or dewatering during diagenetic reactions in layers or pumping in of fluids into layers for energy storage or CCS. Figure 6(a) shows the fluid pressure and fracturing in the horizontal layer with the development of a horizontally aligned branching fracture in the upper part of the layer. The branching fractures connect to form a thin brecciated zone. The fracture is not exactly in the middle of the layer because of the fluid pressure difference below and above the layer. Figures 6(b) and 6(c) show the same horizontally aligned layer with an offset in the middle of the simulation where the layer in Figure 6(b) is offset by half of the thickness of the layer and the one in Figure 6(c) by the full thickness of the layer. The offset of the horizontal layers could, for example, be produced by a fault. Note however that in the current setup, this potential “fault” has no distinct properties other than offsetting the layer. The fracture pattern is significantly different from the nonoffset horizontal layer in the sense that two types of fractures develop: vertical or conjugate shear fractures and horizontal fractures. Vertical or conjugate shear fractures develop early at the offset end of the layer where a potential fault would be positioned and within the layer at localities where the local fluid pressure is very variable. Horizontal fractures within the layer develop later and become more irregular and less horizontally aligned than those in the unfaulted layer. The fractures within the layer shown in Figure 6(b) are still horizontally aligned but branch more than those of the unfaulted layer and connect across the fault. The fractures shown in Figure 6(c) are even more branched and brecciated and show only a semihorizontal alignment in the lower layer on the right-hand side. The fractures of the two faulted layer parts still connect across the fault.
Figure 7 shows the progressive development of fractures in the faulted layer and also differences that occur when Young’s modulus is varied with Figure 7(a) showing a soft material and Figure 7(b) a tough material (note that the layer and the host rock have the exact same Young’s modulus). The progressive fracture development shows that the fractures at the fault where the layers end develop first, followed by semihorizontal fracturing of the layers themselves and finally a brecciation of the layers. The fluid pressure in the offset layers also leads to the progressive development of a shear fracture that propagates into the surrounding matrix at the position where the potential fault would be. The differences in the variable Young’s modulus are mainly illustrated by the timing of the onset of fracturing, which is much earlier and thus at much lower fluid pressures in the soft material than in the hard material. Otherwise, the developing fracture patterns are very similar with a minor increase in fracture density in the hard material.
Figure 8 shows the porosity in 6 simulations with offset layers with different Young’s moduli. In all six cases, the porous channels in the horizontal layers connect across the fault and a compressed low-porosity zone develops in the hanging and the footwall of the fault next to the layers. The compressed zone is more pronounced in the soft than in the hard material. Another variation with a change in Young’s modulus is the timing with the mature porosity channels or waves developing at a much earlier stage in the soft materials than in the hard materials. The other main difference is that the soft materials develop a porosity channel that progresses into the fault outside of the layers, a feature that cannot be seen in the hard materials.
4. Discussion
4.1. Importance of Geometry of the Fluid Overpressure Zone
Our simulations show that the development of hydrofractures and effective stress fields depends on the geometry of the fluid overpressure zone and boundary conditions in simulations, by inference in experiments and nature. Scenario 2 mimics the classical effective stress model where the Mohr circle moves towards the left-hand side and reaches effective tensile stresses without changes in differential stress and thus Mohr circle radius. Several authors argued that this model is an oversimplification and cannot be applied in all cases in the Earth’s crust [10–15, 20, 23]. Hillis [12] and Cobbold and Rodrigues [20] describe our scenario 1 with the effective stresses becoming zero and a stress orientation flip so that the vertical stress becomes the smallest [23]. The stress flip happens because the vertical and horizontal stresses are linked if the system is vertically endless or confined by walls. That means that fluid overpressure reduces vertical and horizontal stresses at different rates bringing the overall effective stress to zero in the vertical and horizontal direction. Afterwards, the fluid can “push” upwards, the vertical effective stress becomes the smallest principal stress, and a horizontal hydrofracture or brecciated zone can develop. This leads to the so called beef veins (Figure 1(b), [20]), which are aligned horizontally. The horizontal stress is a function of the vertical stress and the fluid pressure (equation (3), [20]), a scenario that has been reported in sedimentary basins below seals [12]. In such a system, the fluid overpressure can become very high without the development of fractures (Figures 3 and 4). Once rock failure takes place, the fractured zones are aligned horizontally, or subhorizontally as in our simulations (Figure 3). In order to predict effective stresses and the fracture pattern accurately, fluid overpressure differences or fluid pressure gradients need to be taken into account. A simple rise in fluid pressure under a seal or in a layer will produce scenario 1. Figure 9 illustrates the difference in pattern formation with scenario 1 shown in Figure 9(a) (i). Here, fluid pressure increases either in the whole basin or in a layer. The gradients are aligned vertically, and the solid in the basin is pushed upwards and the layer is dilating. This will produce horizontally aligned hydrofractures, breccias, or beef veins. In this case, the symmetry of the system is only broken in one orientation due to the developing horizontally aligned fluid pressure gradients.
(a)
(b)
Scenario 2 will only develop if the zone where the fluid pressure rises has a vertical or subvertical boundary. The effect of pressure gradients is also illustrated in scenario 2 where the initial fractures develop as a function of the local fluid overpressure gradients at the rim of the high-pressure zone and then propagate progressively into the center of the model. The high-pressure zone is then dilating to form a stationary porosity wave with local mobile porosity channels, where different hydrofractures open and close depending on local fluid pressure variations. Because scenario 2 represents a zone of high pressure, the symmetry is broken in two orientations and the gradients develop around the high-pressure zone (Figure 9(a) (ii)). Now, the differential stress plays an important role and the zone will be dilating more in the direction of the lowest principal stress and thus form initially vertically aligned fractures. In addition, fractures in this zone will form due to two reasons: dilation that will produce a central fracture or breccia within the zone and the fluid pressure gradients from the zone to the surroundings (Figure 9(b)). If the gradients are high enough, the very early fractures will form at the rims of the high-pressure zone.
In scenario 1, the fracture only forms once the dilating high fluid pressure zone overcomes the lithostatic stress. This dilating zone moves downwards with an increase in overpressure so that more material will be pushed upwards. The initial position of the hydrofracture in the middle of the high-pressure zone was analytically predicted by Cobbold and Rodrigues [20] with the following equation: with the first term representing the hydrostatic fluid pressure (as a function of density , gravity , and depth ), the second term representing the fluid overpressure buildup below a seal as a function of influx of fluid (with being the Darcy velocity, the fluid viscosity, and the permeability), and the last term representing a source term, which could be fluid generation in a layer or zone as is the case in our simulations (with being the fluid production). The second term in equation (19) represents a linear fluid pressure gradient whereas the source term is quadratic and the ratio between the two terms determines the sharpness of the pressure peak [20]. Equation (19) is one-dimensional and thus is only partly applicable to a two- or three-dimensional fluid overpressure field, even though the second and third terms in the equation would be similar in the other dimensions.
The importance of two dimensions becomes obvious in scenario 3 which illustrates why an understanding of the interplay of local overpressure gradients and dilating high-pressure zones is important. Once the fluid-generating layer is faulted, the fluid pressure gradients and effective stresses become variable in two orientations (Figure 9). This leads to the initial development of fractures at the tips of the layers where they are faulted, because in this direction, the pressure gradients are acting in the direction of the lowest principal stress and are locally very high. This leads to early fracture development due to pressure gradients, which is later on followed by brecciation of the layers themselves which is happening due to the dilation (Figure 9). The developing fractures are initiated by the fluid overpressure in the layers but develop as a function of the effective solid stresses, the local fluid pressure gradients, and the dilating layers. One can also explain the fracture and porosity patterns by symmetry breaking. The initially random porosity has the highest symmetry. Once the symmetry is broken in one direction, the horizontally aligned hydrofractures and breccia zones develop (Figure 9). If the symmetry is broken in two directions (the layer is faulted or the high-pressure zone is localized), the pattern becomes more complicated and fractures arrange in a combination of vertical and horizontal patterns and combinations forming larger breccias.
4.2. Porosity Phase Transition
Our simulations show a transition in porosity from a random distribution to the development of stationary porosity waves with a dilating zone surrounded by compressed areas. This switch in geometry can be seen as a first-order phase transition where the ordered system (random porosity) reduces its symmetry and forms stationary waves. Such stationary waves are thought to form in sedimentary basins where they can form zebra dolomites [45]. The switch from random quasi-homogeneous porosity to more localized porosity is well known from a range of scenarios in geosciences with the development of compaction fronts or porosity waves in experiments, simulations, and observations in natural systems in porous media, sediments, and viscoelastic materials (for example, [46–50]). In our case within the dilating zones, the hydromechanical interaction leads to the creation of local porosity channels that are dynamic and move around. These channels represent the opening of connected hydrofractures, and they can close again if the local overpressure is reduced. The porosity channels are linked to hydrofractures, but not all hydrofractures develop into open channels nor is the development of fractures directly linked to the phase transition. The difference between fractures and porosity channels is illustrated in scenario 2 where the system fractures early on (Figure 4) but the porosity wave only appears later once the high-pressure zone is dilating (Figure 5). In contrast to this difference between fractures and porosity channels, fracture and porosity channel development in scenario 1 happens at the same time.
4.3. Faulting and Fluid Overpressure
Faults play an important role by changing layer geometries leading to the development of complex fluid pressure gradients and thus fracture patterns. If horizontal layers are intact (scenario 3a) or the high fluid pressure zone has the geometry of a simple sedimentary basin (scenario 1), then the layer or the basin fluidizes, the effective stresses approach zero, and horizontal fractures or breccia zones develop. However, if a fault disturbs this geometry (scenarios 3b and 3c), the effective stresses become complicated and a high-pressure zone dilates not only in the vertical direction but also in the horizontal direction. This leads to hydrofracture development within the layer-offsetting fault and to dilation of the fault at least right next to the high fluid pressure layer (Figure 8). The fault would most probably become a leak for high fluid pressure, and the fractures would dynamically open even a sealed fault.
4.4. Material Variations
The presented simulations show that a variation in Young’s modulus leads to an earlier fracturing in soft than in hard materials. The soft material shows stronger dilation than the hard material, so that the stationary porosity wave in the soft material develops under lower fluid overpressures than those in the hard material. This is in contrast to tectonic fractures (for example, during layer-parallel extension) that would develop earlier in the hard material. This difference in behavior between tectonically driven and fluid pressure-driven fracturing may be used as a proxy for fluid overpressure in layered systems.
5. Conclusion
In this contribution, a hydromechanical numerical model with a compressible fluid was used to simulate hydrofracturing and porosity evolution in high fluid pressure zones. Our simulations led to the following conclusions: (1)The geometry of the zone of high fluid pressure has a major control on rock stability, fracture patterns, fracture evolution, and porosity channeling(2)Hydrofracturing is driven by two main processes: dilation of high-pressure zones and high fluid pressure gradients(3)High fluid pressure generation in a sedimentary basin or a horizontal layer can lead to a reduction of both differential and mean stresses, followed by a switch of the lowest effective stress from a horizontal to vertical and a horizontal hydrofracture. Such a scenario produces high fluid overpressures and a horizontal porosity channel that progressively moves downwards(4)High fluid pressure generation in a confined zone surrounded by rocks with hydrostatic fluid pressure leads to a reduction of the mean stress and early fracturing due to pressure gradients followed by dilation in the direction of the lowest effective stress. During later stages, a central breccia zone with vertically aligned dynamic porosity channels follows(5)Horizontal layers that are offset by a fault and that develop internal high fluid pressures show a combination of scenarios with the fault developing into an early hydrofracture due to pressure gradients followed by dilation of the layer pieces leading to layer-parallel and perpendicular fractures and brecciation. Porosity channels that form in the layers connect across the faults(6)All simulations show a phase transition from random porosity to stationary and nonstationary porosity waves with dilating and compressed zones. In addition, dynamic porosity channels develop within the dilating zones when hydrofractures open and close(7)Fluid can best escape in vertical high-pressure zones with vertical porosity channels and in the faulted layer case where the fault itself develops into a porosity channel(8)In order to understand fluid pressure evolution and resulting fractures and drainage in the Earth’s crust, fluid pressure gradients and exact geometries of seals, fluid sources, and faults have to be taken into account
Data Availability
The simulation input and output data used to support the findings of this study are available from the corresponding author upon request. The basic software for the simulations can be found and downloaded at http://elle.ws, and the corresponding author will make the additional code available upon request.
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
RT acknowledges the support of the INSU ALEAS program of the LIA France-Norway D-FFRACT, of the Research Council of Norway through its Centres of Excellence funding scheme (project number 262644), and DK and RT acknowledge the support of the ITN FlowTrans, a funding from the European Union’s Seventh Framework Programme FP7 for research under grant agreement no. 316889.