Fracturing and Porosity Channeling in Fluid Overpressure Zones in the Shallow Earth’s Crust

GeoZentrum Nordbayern, University Erlangen-Nuremberg, Schlossgarten 5, 91054 Erlangen, Germany School of Earth and Environment, Institute of Geophysics and Tectonics, The University of Leeds, Leeds LS2 9JT, UK Eberhard Karls Universität Tübingen, Germany Institut de Physique du Globe de Strasbourg, UMR 7516, Université de Strasbourg/EOST, CNRS, 5 rue René Descartes, 67084 Strasbourg Cedex, France SFF PoreLab, The Njord Centre, Department of Physics, University of Oslo, P.O. Box 1048 Blindern, NO-0316 Oslo, Norway


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][2][3][4][5][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 σ ij ′ according to Terzaghi's law [9]: with σ ij being the total normal stress, P f the pore fluid pressure, and δ ij the Kronecker delta with the sign convention of positive compressive stress. 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][11][12][13][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][17][18][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 λ v is added that represents the loss of fluid from a pore into the surrounding rock [21,22] according to with the pore fluid factor λ v representing the ratio of pore fluid pressure (P f ) to vertical stress (σ v ), α P the poroelastic factor, and v 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 2 Geofluids 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 P is the local overpressure and K e 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.
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.
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 leftand 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.

Governing
Equations. The fluid phase is described by the fluid pressure P f 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 ρ s and ρ f are the solid and fluid densities, respectively; u s and u f 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 ϕ u f 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 K is calculated from the local porosity according to the Kozeny-Carman relation [36]: with r 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: 3 Geofluids with ρ 0 being the fluid density at some reference pressure. When ρ f and u f are substituted into equation (5), ∂ t ϕ 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 dV, with a solid of mass dm = ρ s dV, and interparticle force f e , fluid force f p , and gravitational loading f g , so that The normal (f n ) and shear force (f s ) acting on a particle from its neighbor is calculated from the relative displacements of the neighbor (i) normal (Δu i n ) and tangential (Δu i s ) to the connecting spring with reference to an equilibrium position. The total force over all connected springs is [38] f where k i n and k i s 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 Fluid Initial configuration (geometry, gravity, P f ) Boundary condition (add P f ) Fluid pressure evolution Time step Elastic relaxation 4 Geofluids 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 ∇P is the fluid pressure gradient, ρ n is the solid fraction, k is running over the four fluid grid nodes near the particle, and sðr i − r k Þ 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 g, the real volume of the particle V R , and the scale factor w (w = 0:74, [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 (U tot ) in the elastic network is where U t and U s are the strain energies for tension and shear, respectively. With the critical strain energy for failure being E ct and E cs for tension and shear, respectively, then where σ n is the normal stress, τ is the shear stress, σ 0 is the tensile strength, and τ 0 is the shear cohesion of the material. The equation typically describes an ellipse in the σ n − τ space [40].

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 ρ ðr 0 Þ or solid fraction in a fluid cell is calculated as and the overall solid velocity u ðr 0 Þ in a fluid cell is where subscript i stands for the particle number and runs through all particles n 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.

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 leftand 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.

Geofluids
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 v ! 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 l being the typical grain or fracture channel size considered and μ the dynamic water viscosity. With a fluid density ρ f of 1000 kg/m 3 , a channel size l 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 Re < 1:0.  (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 1000 × 1000 m (1.0 in model dimension), the overburden is 1000 m, the overburden density is 2400 kg/m 3 , 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/m 3 , the fluid viscosity is 0.001 Pa·s, the fluid compressibility is 4:5 × 10 −10 Pa −1 , 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 t 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).

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 highporosity 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.
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 righthand side having a normal hydrostatic pressure and the resulting effective solid stress of σ 3 . 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 x and y 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 6 Geofluids   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 leftand 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.   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.  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 righthand 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  Figure 6: Fluid pressure (in color) and fracturing (in black) for scenario 3 with fluid pressure increase in a horizontal layer that is faulted in (b) and (c). A horizontal layer shows a fracture evolution that is similar to scenario 2 with a semihorizontal branching fracture in the layer. Once the layer is faulted, the fluid pressure gradients curve and form vertical to conjugate shear fractures at the layer tips where the fault is positioned. Fracturing within the layer is also more complicated because the layer is dilating in two dimensions in the vertical as well as the horizontal direction. This leads to stronger brecciation of the layer.
11 Geofluids 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.

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.
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 highpressure 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 g, and depth z), the second term representing the fluid overpressure buildup below a seal as a function of influx of fluid (with q being the Darcy velocity, μ the fluid viscosity, and k 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 Q 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 highpressure 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 13 Geofluids 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.

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][47][48][49][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 14 Geofluids and porosity channels, fracture and porosity channel development in scenario 1 happens at the same time.

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.

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.

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.