Steady-State Flow Through a Subsurface Reservoir with a Displaced Fault and its Poro-elastic Effects on Fault Stresses

We consider steady-state single-phase confined flow through a subsurface porous layer containing a displaced, fully conductive fault causing a sudden jump in the flow path, and we employ (semi-)analytical techniques to compute the corresponding pressures and fault stresses. In particular, we obtain a new solution for the pressure field with the aid of conformal mapping and a Schwarz–Christoffel transformation. Moreover, we use an existing technique to compute the poro-elastic stress field with the aid of inclusion theory. The additional resistance to fluid flow provided by a displaced fault, relative to the resistance in a layer without a fault, is a function of dip angle, fault throw divided by reservoir height, and reservoir width divided by reservoir height. Fluid flow has a larger effect on fault stresses in case of injection than in case of depletion, where injection with up-dip flow results in increased zones of fault slip near the bottom of the reservoir. Opposedly, injection with down-dip flow results in increased slip near the top of the reservoir. An order-of-magnitude estimate of the effect of steady-state flow across displaced faults in the Groningen natural gas reservoir shows that the effect on fault stresses is probably negligible. A similar estimate of the effect in low-enthalpy geothermal doublets indicates that steady-state flow may possibly play a small role, in particular close to the injector, but site-specific assessments will be necessary to quantify the effect. We present a semi-analytical method to compute the pressure and flow fields in reservoirs with a displaced fault. In case of reservoir depletion, the effect of steady-state flow on the induced fault stresses is expected to be negligible. In case of fluid injection, steady-state flow may significantly affect the induced fault stresses and the potential slip patches. We present a semi-analytical method to compute the pressure and flow fields in reservoirs with a displaced fault. In case of reservoir depletion, the effect of steady-state flow on the induced fault stresses is expected to be negligible. In case of fluid injection, steady-state flow may significantly affect the induced fault stresses and the potential slip patches.


Introduction
Numerous deep-subsurface activities have been associated with induced seismicity, including geothermal energy production, hydrocarbon production, underground storage of natural gas, CO 2 or hydrogen, and wastewater injection (Segall 1989;Suckale 2009;Elsworth et al. 2016;Foulger et al. 2018;Muntendam-Bos et al. 2022).Induced seismicity may lead to serious societal concerns and damage to housing and infrastructure.This has caused the early termination of several subsurface projects, such as the Groningen gas field in the Netherlands (Muntendam-Bos et al. 2022) and the enhanced geothermal project in Basel, Switzerland (Häring et al. 2008;Mignan et al. 2015).Therefore, understanding the factors leading to induced seismicity is crucial.
Changes in pore pressure caused by fluid injection or production alter the stresses in the rock through poro-elastic effects and may cause slip along existing faults.Such faults have typically been generated through deformations driven by plate tectonics over geological timescales, and are nearly always displaced, i.e., they have a nonzero offset (a.k.a. as fault throw); see Fig. 1 which displays a reservoir intersected by a so-called normal fault where tectonic extension in combination with gravity has caused the rock to the left of the inclined fault to move down relative to the rock to the right of the fault.Figure 2 displays the same reservoir but now intersected by a so-called reverse fault, a result from tectonic compression which caused the rock to the right of the fault to be pushed up relative to the rock to the left.
Fig. 1 Schematic representation of a porous and permeable reservoir with height h and width w, embedded in a non-permeable domain and intersected by a displaced normal fault.The fault is characterized by dip and offset t f and is located at the center of the reservoir.The combined pore pressure P(x, y, t) = p 0 (y) + p(x, y, t) is assumed to be quasi steady state with prescribed values P L (y, t) = p 0 (y) + p L (t) and P R (y, t) = p 0 (y) + p R (t) at the left-and right-hand sides, respectively For a nonzero fault throw, pore pressure changes induce large shear stresses at the reservoir corners along the fault.The importance of the fault throw for the onset of seismicity was first recognized by Roest and Kuilman (1994).This effect was recently investigated in more detail by van den Bogert (2015), Buijze et al. (2017), van Wees et al. (2017), van den Bogert (2018), and Buijze et al. (2019a) through numerical studies.Such numerical models can account for complexities such as spatial heterogeneities, various friction laws, and multi-phase flow.However, these models are computationally expensive, partially because adequately resolving the large stress peaks at the reservoir corners requires an extremely fine mesh.
Analytical or semi-analytical solutions can provide a fast estimate of the stresses induced by pore pressure changes.These solutions cannot handle the complexities that numerical models can, but the advantage is that the stress peaks can be resolved in great detail (Wu et al. 2021).Commonly used analytical methods are inclusion theory (Eshelby 1957;Rudnicki 2002) and the closely related nucleus of strain concept (Goodier 1937;Geertsma 1966Geertsma , 1973)).Inclusion theory has been applied to investigate the poro-elastic effects of incremental reservoir pressures caused by injection or depletion (Segall 1985(Segall , 1989;;Soltanzadeh and Hawkes 2008), or of thermal stresses (Segall and Fitzgerald 1998).More recently, the effect of the fault throw on Coulomb stresses and fault slip was investigated (semi-)analytically (Jansen et al. 2019;van Wees et al. 2019;Wu et al. 2021;Jansen and Meulenbroek 2022).
However, all these studies consider a spatially uniform near-steady-state increase or decrease in reservoir pressure and disregard the effect of pressure gradients and the corresponding flow across the fault.Postma and Jansen (2018) considered the effect of pressure transients on fault stresses but did not take into account the particular effect of fault throw.Anderson (2006) performed an in-depth analysis of the flow in a reservoir with a displaced fault but focused on flow in the horizontal plane (along and across the fault) while using a Dupuit assumption for the vertical flow.Zbinden et al. (2017) addressed the effects of flow in a displaced fault with such a large throw that the two reservoir blocks no longer share a common boundary.
Here, we will address the effect of steady single-phase flow across a fully conductive displaced fault between reservoir blocks in direct communication, while focusing on the Fig. 2 Schematic representation of a porous and permeable reservoir intersected by a displaced reverse fault.All other details are as stated in Fig. 1 vertical and horizontal flow as caused by the "jump" in the flow conduit.In particular, we aim to quantify the pressure field in the reservoir and its effects on the stresses in the fault.We pursue a semi-analytical approach such that the stress peaks at the reservoir corners can be computed fast and accurately.We consider a configuration with two wells, placed, respectively, on the left-hand and right-hand sides of the reservoir, indicated by P L and P R in Figs. 1 and 2. Such a setting is relevant for, but not restricted to, geothermal doublets.

Geometry
We consider a two-dimensional (2D) faulted reservoir with height h and width w and a geometry as depicted in Fig. 1.The fault is characterized by dip and fault throw t f .In this paper, we specify dip angles up to 180°, which differs from the convention in geology to use dip angles of up to 90° in combination with a dip direction.Thus, our definition of a fault dip of 110° is equivalent to a dip of 70° in the opposite direction.The impermeable domain outside the reservoir extends to infinity.Locations inside and outside the reservoir are indicated with an x − y coordinate system with its origin at the reservoir center.Two wells penetrate the reservoir, one on the left-hand side and one on the right-hand side of the domain.We assume that the pore pressure can be controlled uniformly across the entire depth of the reservoir.The parameter values used in this study are mostly identical to those used in Jansen and Meulenbroek (2022) and are roughly based on the properties of the Groningen natural gas reservoir in The Netherlands; see Table 1.

Incremental Pore Pressure
Several assumptions are required to obtain a semi-analytical solution for the incremental pressure field p(x, y, t).We consider steady-state single-phase flow and assume that the permeability of the surrounding rock is negligible compared to the permeability of the reservoir, such that fluid flow outside the reservoir can be disregarded.Additionally, the fault is assumed to be infinitely thin and completely conductive.We further assume that Darcy's law is applicable and that the permeability is homogeneous inside the reservoir.Under these assumptions, the governing equations for fluid flow reduce to Laplace's equation.Quasi-steady-state incremental pressures are prescribed at the vertical boundaries on the left-and right-hand sides of the reservoir, denoted by p L (t) and p R (t) , respectively (Fig. 1).All other boundaries are no-flow boundaries.A fully analytical solution seems to be out of reach for this problem, but a semi-analytical solution can be found as follows.Since the governing equation is the Laplace equation, the geometry of the faulted reservoir can be transformed to a simpler domain through conformal mapping.The flow problem in the faulted reservoir somewhat resembles the problem of confined aquifer flow over a stepped base as considered by Bakker and Post, (2022, p.181).In the stepped base problem, the pressure field can indeed also be solved through conformal mapping.However, already in that simpler problem the inverse transformation (from the complex potential plane to the physical plane) cannot be obtained analytically.We face a similar difficulty while, moreover, the geometry of the flow domain is more complex such that we need to revert to the Schwarz-Christoffel method to transform the problem domain.A detailed description of the solution strategy is given in Appendix A.
The result can be expressed in terms of the average incremental pressure p(t) and the pressure drop Δp = p R (t) − p L (t) , which yields where the exact definition of J(x, y) is given in Eq. (A18) in Appendix A. Note that this definition tacitly assumes that the pressure drop Δp is time-independent, whereas the aver- age incremental pressure p(t) is a (slow) function of time, with positive and negative val- ues of p corresponding to injection and depletion, respectively.We will drop the explicit dependence on time from the notation in the remainder of this paper.
The semi-analytical solution can also be applied to anisotropic reservoirs.Often, the vertical permeability is smaller than the horizontal permeability due to laminations within the reservoir.This difference in flow resistance can be represented by an equivalent isotropic domain with scaled height and width (Muskat 1937;Bakker and Post 2022).This transforms the anisotropic flow equation to Laplace's equation, which allows the application of the semi-analytical solution presented here to the scaled domain.
In practice, the Schwarz-Christoffel algorithm shows convergence issues when some reservoir edges are much larger than other edges.This issue arises for our reservoir geometry since w ≫ h .To minimize such convergence issues, we only apply the Schwarz-Christ- offel transformation to the region −2 ≤ x ≤ 2 and assume that the streamlines are paral- lel at this distance from the fault.The justification of this assumption and discretization of the subdomains is shown in Appendix A.4.The threshold of ±2h serves to make the Schwarz-Christoffel solution applicable to a large range of dip angles and fault throws, but this condition can be relaxed for some combinations of dip angle and fault throw.

Incremental Stresses
The stresses induced by pore pressure changes are computed using inclusion theory, similar to the work of Jansen et al. (2019) and Wu et al. (2021).The reservoir is treated as an inclusion in an infinite space with homogeneous elastic properties.Changes in pore pressure cause compaction or expansion of the reservoir.Since the reservoir is constrained by the surrounding rock, stresses develop in-and outside of the inclusion.Here, we only give a brief summary of the governing equations.A more rigorous derivation is given in Appendix B.
(1) p(x, y, t) = p(t) + J(x, y)Δp The displacements are computed by integrating Green's functions over the inclusion domain.These Green's functions describe the solution of the mechanical equilibrium equation for the displacements at a point (x, y) due to a change in pore pressure at a point ( , ) .The main difference with the analytical work of Jansen et al. (2019) and Wu et al. (2021) is that the pressure is now non-uniform and should be included inside the integral, an approach that was earlier taken by Segall (1985).The expression for the incremental displacements then becomes where D = (1−2 ) 2 (1− )G , is Biot's coefficient, is Poisson's ratio, and G the shear modulus.Expressions for the Green's functions g i are given in Appendix B. The strains are obtained from taking the spatial derivatives of the displacements The incremental stresses are computed from the strains using Hooke's law.Outside the inclusion, the total strain ij can be used directly, but inside the inclusion the contribution of the eigenstress needs to be subtracted.For isotropic elastic properties, this yields where is Lamé's first parameter, ij is the Kronecker delta, and Ω is a modified Kronecker delta which equals 1 inside the reservoir and 0 outside.The subscript kk implies summation.The incremental stresses given by Eq. ( 4) are transformed to obtain the incremental normal and shear stresses at the fault according to whereafter the incremental effective normal stress on the fault is given by Fault slip is governed by combined (i.e., initial plus incremental) shear and effective normal stresses, which are defined as where superscripts 0 indicate initial values.Fault slip occurs when where Σ sl is the slip stress, defined as with c indicating fault cohesion and st the static friction coefficient inside the fault, and where it should be kept in mind that a negative normal stress corresponds to compression.Alternatively, the criterion for fault slip can be expressed in terms of the combined pre-slip Coulomb stress in which case slip corresponds to positive values of Σ C .
Here, we only consider incremental stresses while we assume cohesion to be absent.In that case we can define the incremental Coulomb stress as Regions with positive incremental Coulomb stresses (  C > 0 ) indicate zones where the fault moves to failure.This does not necessarily mean that slip will occur here, as this also depends on the initial Coulomb stress distribution along the fault.
The integral for the displacements in Eq. ( 2) can be split into a component governed by the average incremental pressure p and a component affected by the pressure drop Δp .The result can be written as where the dimensionless pressure drop p is defined as Equation ( 14) combined with Eq. ( 3) can be used to assess the significance of fluid flow on the strains and stresses in the entire domain while its subsequent application in Eqs. ( 4) to (7) allows for an assessment of the effect of flow on the (pre-slip) fault stresses.Equation ( 14) depends on the geometrical factor J( , ) , which is discussed in the next section, and on the relative pressure drop p .When p approaches zero, the expression for the strain reduces to that for a uniform pressure field.In this case, fully analytical solutions for the displacements, strains, and stresses are available from Jansen et al. (2019) while corresponding compact expressions for the fault stresses, and a semi-analytical method to determine the resulting fault slip, can be found in Jansen and Meulenbroek (2022).

Pressure Field
Figure 3 shows an example of a typical incremental pressure field in a reservoir with a displaced normal fault.For this example, we used a fault with dip angle = 70 • , a scaled offset t = t f h = 75 225 = 1 3 , and a reservoir width w = 20h which corresponds to a distance of 4500 m between the fixed pressure boundaries.We applied a negative unit pressure drop Δp = −1 MPa such that flow is in the "up-dip" direction (i.e., from left to right in Fig. 3).As seen in Eq. ( 1), the pressure field can be expressed as an average pressure plus (11) a deviation caused by the pressure drop.Therefore, we plotted a scaled pressure defined as (p(x, y) − p)∕|Δp| , which equals −J(x, y) in the case of a negative pressure drop.Hence, the figure shows the deviations from the average pressure for a unit decline in pressure.At large distances from the fault, the flow field is one-dimensional, i.e., the pressure drop is linear in the horizontal direction, as indicated by the vertical pressure contours and horizontal streamlines.Near the fault, however, the fluid must travel around the impermeable formations surrounding the reservoir which leads to a significant vertical flow component in this region.
The pressure deviations scale linearly with the applied pressure difference Δp , but are also dependent on the reservoir geometry through the factor J(x, y).The largest deviations occur at the edges of the reservoir at x = ±2250 m, with a magnitude of ±0.5Δp (outside the figure).The function J(x, y) thus adheres to the constraint −0.5 ≤ J ≤ 0.5 .It approaches zero near the center of the domain, where the pressure is closer to the average pressure, which coincides with the location of the fault in our setup.This means that the term inside the second integral term of Eq. (B25) will always be smaller than the first integral.However, the contribution of this second term may still be relevant, especially for large values of the relative pressure drop p.
While the stresses at the fault are affected by pressure changes in the whole reservoir, the contributions of these changes diminish with their distance from the fault.In other words, the stresses along the fault are most sensitive to pressure changes close to the fault, while their values are also affected by the fault dip, offset, and reservoir width.The function value J(x, y) at the upper external corner (e.g., at J(x, y) = J(54.6,150) for the example of Fig. 3) is plotted as a function of the fault offset and dip in Fig. 4. A smaller dip corresponds to a fault that is more aligned to the main flow direction.This leads to larger pressure differences along the fault.A larger fault offset leads to a narrower opening between the left-and right-hand sides of the reservoir, which leads to larger pressure gradients along the fault.The function J is antisymmetric along the fault, so the absolute value of J is the same in the lower external corner, but with opposite sign.Figure 5 shows an example of a typical incremental pressure field in a reservoir with a displaced reverse fault.All parameter values in this example are identical to those used to generate Fig. 3, except for the dip angle which is now taken as = 120 • .Just like in Fig. 3, the flowlines close to the fault are forced to deviate from their undisturbed parallel horizontal trajectories.However, a difference between the configurations displayed in the two figures is that a normal fault causes the height of the reservoir to be locally reduced, whereas the reverse fault causes the height to be locally increased.

Incremental Flow Resistance
For a given pressure drop Δp , a reservoir without fault will provide a resistance to flow R 0 resulting in a flow rate q 0 , whereas a reservoir with a displaced fault provides a larger resistance R resulting in a smaller flow rate q.We define the incremental resistance according to where is permeability, is viscosity, ŵ = w h is scaled reservoir width, and f is a geometry-dependent factor defined in detail in Appendix A.4.The contribution of the reservoir geometry to the resistance is bundled in the dimensionless parameter R .Figure 6 displays the dimensionless resistance R as a function of and t for a fixed value ŵ = 4500 225 = 20 .We verified with a numerical model that the influence of ŵ on the added resistance is very small for ŵ > 2 , and therefore Fig. 6 will be valid for most practical situations.The range of dip angles (60 • <  < 120 • ) in the figure is typically encountered in normal and reverse faulting regimes, and it can be seen that for these values the magnitude of the incremental resistance is primarily dependent on the scaled fault throw, and only to a lesser extent on the fault angle.The range of dip angles here does not include shallow-dipping features such as thrust faults.While the Schwarz-Christoffel method is applicable to such geometries, convergence issues arise for shallow dip angles, which we will elaborate on in the Discussion section.
Somewhat surprisingly, the incremental resistance becomes slightly negative for high dip angles (corresponding to a reverse fault) and small values of t , although this is hardly visible in the figure.An explanation for this effect is the locally increased reservoir height (see Fig. 5), which, for small values of the scaled fault throw, results in a local reduction in flow resistance that just overcompensates the increased flow resistance caused by flow line disturbance.
If we consider depletion in a hydrocarbon reservoir, the absolute value of the average incremental pressure p is small at the start of production, and the relative pressure drop, which was defined earlier as p = Δp p , is therefore relatively large.As the reservoir is increasingly depleted, the absolute value of the average incremental pressure increases and the effect of fluid flow on fault stresses will therefore become relatively less significant with increase in reservoir depletion.On the other hand, in a typical geothermal doublet setting, the average pressure drop will be close to zero, because all the produced water is re-injected (although temperature effects may cause a slight shrinkage).In that case the effect of flow on the fault stresses should be compared to the effect of thermal stresses caused by the lower temperature of the re-injected water.We will elaborate on the practical consequences of these relative effects in the Discussion section.pressure drop p = Δp p for a depletion situation (p < 0) with up-dip flow (i.e., flow from left to right).The stress is scaled by a factor C = (1−2)|p| 2(1−) .The parameters used for this figure are given in Table 1.The line p = 0 coincides with a uniform incremental pressure field p = p , for which fully analytical pressure solutions are available (Jansen et al. 2019;Wu et al. 2021).Due to differential compaction between the reservoir and the surrounding rock, large shear stresses develop at the reservoir-fault corners.In case of depletion, the incremental Coulomb stress is positive between the "internal" reservoir-fault corners (i.e., for −75 ≤ y ≤ 75 m), marked in green in Fig. 7 and thus the fault moves toward failure in this region.However, this only represents the change in Coulomb stress.Whether fault slip will occur also depends on the initial Coulomb stress distribution.

Stresses Along the Fault
For a uniform incremental pressure field ( p = 0 ), the incremental stress profile is sym- metric around the x-axis.The pressure and stress profiles become (more) asymmetric when the absolute value of the pressure drop increases.In case of depletion with up-dip flow, the pressure in the upper half of the reservoir is more negative (i.e., the reservoir is more depleted here) than in the bottom half, and this asymmetry becomes more pronounced with increasing p ; see Fig. 7 (left).However, as can be seen from Fig. 7 (right), the corresponding effect of flow on the positive values of the pre-slip Coulomb stresses around y = ±75 m is very small.In case of depletion with down-dip flow (i.e., flow from right to left), the additional pressures change sign and the Coulomb stress profile is mirrored in the x-axis, but the effect remains negligible.
Figure 8 depicts the scaled pressure and incremental Coulomb stresses for a similar situation as Fig. 7, but now for an injection situation.The flow is still from left to right, but the sign of the average pressure is now positive, meaning that the dimensionless pressure drop p is now negative.The regions with positive incremental Coulomb stress now lie outside the reservoir (i.e., y < −75 m and y > 75 m).flow (from right to left) yields a similar result, although mirrored in the x-axis, such that the growth of the zones of positive incremental Coulomb stresses then occurs in the upper region.
Figures 7 and 8 both depict stress changes in a normal fault.Similar figures can be produced for a reverse fault, but they will qualitatively not be very much different from those for normal faults.Furthermore, the effects will be less severe in the sense that the combined Coulomb stresses in reverse faults will only become positive in very small regions.This is because reverse faults typically experience comparatively large compressive initial normal stresses as a remainder of the compressive tectonic stresses that caused the formation of the reverse faults (Jansen et al. 2019).
We note that positive incremental Coulomb stresses give an indication of where slip patches may develop, but do not tell the full story.The potential slip patches are also dependent on the initial Coulomb stress distribution.Once slip occurs, additional shear stresses occur along the fault leading to a shear stress redistribution and a corresponding modified post-slip Coulomb stress pattern (Jansen and Meulenbroek 2022).Moreover, the occurrence of dilation may introduce additional normal stresses and thus further complicate the post-slip stress pattern.Furthermore, slip may occur aseismically, i.e., gradually; or co-seismically, i.e., in the form of a seismic event (a.k.a. an earthquake).The nucleation of a seismic event typically occurs after growth of a slip patch over a certain distance or after the growth rate reaches a certain limit, depending on the details of the friction law; see, e.g., Uenishi and Rice (2003); Ampuero and Rubin (2008).A detailed analysis of flow-induced seismicity is well beyond the scope of the present paper, but from the results presented above we can safely draw the conclusion that fluid flow in displaced faults has a much larger effect on slip patch growth in the case of injection compared to the case of depletion.

Consequences for Practice
As a practical depletion example, we consider the Groningen natural gas field in the Netherlands (van Elk et al. 2021) which contains a large number of faults with a wide range of offsets.A typical pressure difference over a distance of about 5 km is Δp = −0.5 MPa (Pijpers 2018), while the average incremental pressure, of which probably only a small fraction is caused by displaced faults with t < 1 , is at the moment typically around p = −25 MPa.This results in p = 0.02 , which shows that the relative contribution of pressure drop to the incremental pressure is currently very small.Close to the wells (or well clusters) much stronger gradients may occur.In that case we should also take into account the strongly nonlinear relationship between pressure and volume in gasses, e.g., through the use of a pseudo gas formulation (Postma and Jansen 2018).However, even for these potentially higher relative pressure drops, it still holds that the effect of flow on the pre-slip Coulomb stresses in case of depletion is very small such that steady-state flow is therefore unlikely to significantly affect fault slip in the Groningen reservoir under the current conditions.Note, however, that for displaced faults with t > 1 another, more significant effect may be at play, as will be discussed below.
As a geothermal example, we consider a typical low-enthalpy geothermal doublet in the Netherlands (i.e., a pair of wells for hot water production at temperatures below 90 • C and subsequent re-injection of the cooled-down water in the same reservoir layer (Buijze et al. 2019b).This setup may have a similar reservoir height as the example described above, while the distance between the injector and producer will be between 1 and 2 km with a pressure drop between the wells up to 5 MPa of which typically only a small fraction will be caused by displaced faults (Daniilidis et al. 2021).The re-injection of cold water will cause shrinkage of the reservoir rock which has a similar effect as a pressure decrease due to production.In case of a temperature drop ΔT = 50 • , the equivalent average incremental reservoir pressure peq can be expressed as (Segall and Fitzgerald 1998)   where T is the linear thermal expansion coefficient.With the data given in Table 1 this results in an equivalent depletion peq ≈ −24 MPa, which is similar to the case of depletion in the Groningen gas field.The effect of steady-state flow on the fault stresses is therefore much smaller than the effect of cooling.However, steady-state flow may be significant when the cold front has not yet reached the fault, in particular for faults located close to the injector.Detailed site-specific studies will be necessary to quantify these potential effects with more certainty.
It should be noted that in the Groningen reservoir, and in typical geothermal doublets in the Netherlands, it has been found that faults are not critically stressed (van Wees et al. 2014;Dempsey and Suckale 2017;Muntendam-Bos et al. 2022).However, in tectonically more active regions, critically stressed faults do frequently occur in which case even the relatively small effects of steady-state flow might indeed lead to the reactivation of (near-) critically stressed faults.

Limitations
In this study, we assumed that the fault was completely conductive (i.e., that the permeability in the fault is the same as in the reservoir rock).In reality, faults may be partially conductive or completely sealing (Caine et al. 1996).A non-conductive fault would lead to two uniform pressure fields on either side of the fault, most likely of different magnitude.This was considered by Wu et al. (2021), and such a configuration can also be handled with our approach.A partially conductive fault, which would lead to a jump in the pressure between both sides of the fault, is more complex to analyze analytically.Anderson (2006) developed an analytical solution for the jump by treating the fault as an inhomogeneity.However, he used the Dupuit approximation and hence neglected vertical flow in the reservoir which may significantly affect the pressure along a displaced fault.Due to the complexity of the reservoir geometry, a (semi-)analytical solution for a partially conductive displaced fault seems to be out of reach.
Fault flow in a slightly different configuration was investigated numerically by Zbinden et al. (2017) who addressed the effects of flow in a displaced fault with a scaled throw t > 1 such that the top of the hanging wall is below the bottom of the foot wall.In case of a non-conductive fault, this would imply complete hydraulic separation between the reservoir blocks to the left and to the right of the fault.However, Zbinden et al. (2017) modeled the fault as conductive such that it provides a pathway for fluid flow between the blocks, and they concluded that for this configuration, where fluid is forced through a very narrow pathway, flow plays a major role in pore pressure and stress evolution within the fault, especially in case of differential depletion.
Another possible effect, not considered in the present study, is a change in fault strength that may result from flow across or along the fault, in particular when two-phase (gasliquid) effects play a role.Depending on lithology, fluid properties and detailed fault geometry such fault weakening may lead to a significantly accelerated transition from aseismic slip to seismicity; see, e.g., Hunfeld et al. (2019); Snell et al. (2020) and references therein.
We neglected flow outside of the reservoir.In reality, the permeability in the over-and underburden will be nonzero and there will be some pressure diffusion across the boundary of the reservoir.This would somewhat smoothen the discontinuity in the pore pressure across the reservoir boundary (e.g., Figs. 7 (left) and 8 (left)), which would also smoothen the stress peaks at the reservoir corners.However, this effect is expected to be negligible the permeability of the reservoir is much larger than the permeability of its surroundings.
Furthermore, we assumed that the elastic properties of the reservoir and its surroundings are equal and homogeneous.In reality, the over-and underburden are likely to have elastic properties different from the reservoir.Such contrasts in rock stiffness would alter the magnitude of the incremental stresses, but the effect of steady-state flow on the stresses is expected to remain similar.While theoretically there exist analytical methods for treating inhomogeneous inclusions (Mura 1987), these methods are only practical for elliptical inclusions (Ma andKorsunsky 2014, 2022).For other inclusion shapes, it is more convenient to resort to numerical methods.
While the Schwarz-Christoffel method can theoretically be applied to fault offset ranges from 0 < t < 1 and dip angles from 0 • <  < 180 • , in practice the transformation and/or its inverse cannot be solved with appropriate accuracy.This is the case for large fault throws, when the flow opening becomes very narrow, and for shallow dip angles.This limits the application of the Schwarz-Christoffel method to certain reservoir geometries.These limitations can be circumvented by using purely numerical methods, which would also allow for including more complex features, such as heterogeneity or transient flow.The semianalytical approach presented here is therefore not expected to replace existing numerical methods.However, the semi-analytical approach has the advantage of being able to accurately resolve the Coulomb stress peaks.Therefore, it is a helpful tool for assessing changes in Coulomb stresses, in this particular case due to steady-state flow.

Conclusions
We investigated the effect of steady-state flow in a reservoir with a displaced, conductive fault through conformal mapping.Similar to the problem of confined aquifer flow over a stepped base, we needed a semi-analytical approach to compute the inverse transformation (Bakker and Post 2022).Moreover, the more complex geometry of the faulted reservoir required the use of the Schwarz-Christoffel transformation.Inclusion theory was used to compute the stresses along the fault.
The additional resistance to fluid flow caused by a displaced fault, relative to the resistance in a reservoir without fault, is a function of dip angle, scaled fault throw (i.e., fault throw divided by reservoir height) and scaled reservoir width (i.e., reservoir width divided by reservoir height).For dip angles 60 • <  < 120 • , as typically encountered in faults orig- inating from tectonic events, the magnitude of the scaled resistance is primarily dependent on the scaled fault throw, to a lesser extent on the fault angle, and hardly at all on the scaled reservoir width.
The impact of fluid flow across a displaced fault on stresses along the fault depends on both the reservoir geometry and the pressure drop over the fault zone.The effect of flow is most strongly felt at the "external" reservoir-fault corners where the pressure is more influenced than at the "internal" corners, and where therefore also the fault stresses are more strongly influenced.It is known from earlier work that injection results in the development of zones of positive combined Coulomb stresses at the "external" reservoir corners, whereas depletion results in positive incremental Coulomb stresses at the "internal" corners (Jansen et al. 2019).Therefore, the effect of flow is considerably larger in case of injection than in case of depletion.Moreover, in case of injection with up-dip flow the effect is a stronger increase in Coulomb stress at the bottom of the reservoir, whereas in case of injection with down-dip flow the effect is a stronger increase at the top.
As a practical example of depletion in a natural gas field, we considered the Groningen field in the Netherlands.An order-of-magnitude estimate of the effect of flow indicates that steady-state flow across displaced faults is unlikely to significantly affect fault slip in the Groningen reservoir under the current conditions.As an example of geothermal energy production, we considered generic dimensions of typical low-enthalpy geothermal doublets.In this case, we find that steady-state flow may possibly play a small role in the development of stresses in displaced faults, in particular close to the injector, but site-specific assessments will be necessary to quantify the effect.transformation from the -plane to the -plane to be done analytically, if the vertices along the real axis of the -plane are chosen conveniently.Following Verruijt (1970), we choose A = 0 , E = 1 , and H = ∞ , which is consistent with the choice in the previous section.Choosing the origin of the -plane to coincide with = D , the transformation to the -plane is given by where p AH and p DE are the incremental pressures prescribed at the boundaries AH and DE, respectively, F(a|x) is the incomplete elliptic integral of the first kind and K(x) is the com- plete elliptic integral of the first kind, and where The parameter D is obtained from the solution of the Schwarz-Christoffel problem of the previous section.Thus, for any (x, y), we can calculate the corresponding coordinate in the -plane.Then, we compute the corresponding potential through Eq. (A11), which can be converted to the incremental pore pressure p through Eq. (A2), which yields

A.4 Discretization into Subdomains
The Schwarz-Christoffel algorithm typically has convergence issues when some edges of the domain much larger than others.We also run into this issue, as the reservoir width is much larger than the reservoir height and fault throw (i.e., w ≫ h > t) .This problem can be circumvented by only modeling the pressure field up to a distance of 2h from the fault Fig. 11 The maximum ratio of vertical pressure gradient over horizontal pressure gradient as a function of the dimensionless distance from the fault, for various dip angles with the Schwarz-Christoffel method.We separate the reservoir domain into three subdomains (see Fig. 10).Around the fault, in subdomain Ω 2 , the pressure field is computed using the Schwarz-Christoffel transformation.At distances further than 2h from the fault, in subdomains Ω 1 and Ω 3 , flow is assumed to be unidirectional.Therefore, the pressure is constant in the vertical direction and decreases linearly in the horizontal direction in these subdomains.This implies that assume that the effect of the displaced fault has vanished and the streamlines are parallel at a distance 2h from the fault.To verify this assumption, we computed the pressure field in the faulted reservoir numerically with the MATLAB Reservoir Simulation Toolbox (Lie 2019) using a mimetic finite volume scheme.We model the full domain w = 20h discretized into 1 million grid cells.To quantify whether flow is mostly horizontal, we compute the ratio of the vertical pressure gradient to the horizontal pressure gradient at various distances from the fault.At each point along the x-axis, we select 100 points regularly spaced from the bottom to the top of the reservoir, where we numerically approximate the horizontal and vertical pressure gradient and their ratio.We then define the relative vertical pressure gradient r at a chosen distance x k from the fault as the maximum value of these ratios where i denotes the spatial derivative in the x i -direction.In Fig. 11, this number is plotted as a function of the distance from the fault.The relative magnitude of the vertical pressure gradient decays rapidly with distance from the fault.At a distance of 2h from the fault, the vertical pressure gradient is smaller than 0.7% of the horizontal pressure gradient for a dip angle of 60 • and even smaller for larger dip angles.Hence, the assumption that flow is mostly horizontal at a distance of 2h from the fault is reasonable.Note that the choice of ±2h was chosen to make the approach applicable to a large range of dip angles and fault throws.For specific reservoir geometries, it is possible to increase this distance.
We prescribe constant incremental pressures at the reservoir edges, located at a distance of 10h from the fault ( p L and p R ).The pressures at a distance of 2h from the fault can be computed using the principle of continuity.We choose to write the solution as a function of the average incremental pressure p = p L +p R 2 and the prescribed incremental pressure drop Δp = p R − p L .The incremental pressure at a distance of 2h from the fault can then be expressed as where h � = h∕( 1 2 w − 2h) and where f = K(1 − D )∕K( D ) is a factor based on the reservoir geometry, with K(x) again being the complete elliptic integral of the first kind.The incremental pressure field in each subdomain can then be expressed as where J(x, y) is a spatially variable factor given by (A14) r(x k ) = max y p x p , (A17) p(x, y) = p + J(x, y)Δp, Writing this out for each of the strain tensor components yields Note that we take the derivative of Eq. (B20) outside of the integral.Due to the singularity in the Green's functions, Leibniz integral rule cannot be applied and the order of differentiation and integration is not interchangeable for points inside the inclusion (Mura 1987, p. 12).
The stresses can be obtained from the elastic strains using Hooke's law.The elastic strain e ij is defined the difference between the total strain and the eigenstrain The definition of the eigenstrain has already been given in Eq. (B19).Note that the eigenstrain vanishes outside the inclusion.This yields the following expressions for the three stress components

B.2 Integration of the Green's Functions
Due to the complex expression for the pressure field, closed-form expressions for the integrated Green's functions seem out of reach.Instead, we subdivide the reservoir domain   12 Example of a reservoir discretized into rectangular and triangular grid cells.The pressure may vary between grid cells, but the pressure inside each grid cell is assumed to be constant, such that closed-form expressions for the induced stresses may be used into rectangular and triangular grid cells (Fig. 12).If the grid cells are chosen sufficiently small, it is reasonable that the pressure inside each cell is constant, but this pressure may differ between cells.In this case, closed-form expressions are available for the integrated Green's functions in each cell.Using the principle of superposition, the contributions of each individual grid cell may be summed to obtain the displacements induced by pressure changes inside the entire reservoir.We thus approximate the displacements in Eq.B20 as where N is the total number of grid cells and is the double integral of the Green's function for the displacement in the i-direction.The strains are then approximated by taking the spatial derivatives of the displacements with Closed-form expressions for G ij are given by Jansen et al. (2019) and Cornelissen et al. (2023) for rectangular and triangular grid cells.The stresses are then obtained from Eqs. B28-B29.For a detailed derivation, see Cornelissen et al. (2023).
are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 3
Fig. 3 Contour plot of the scaled pressure (p(x, y) − p)∕|Δp| , with Δp = −1 in a reservoir intersected by a normal fault.The scaled pressure is equal to the negative of the geometrical pressure deviation factor J(x, y).The contours of the scaled pore pressure are indicated by dashed lines and the streamlines are indicated by solid lines.A scaled fault throw t = t f ∕h = 1 3 and a dip angle = 70 • have been used in this example

Fig. 4 Fig. 5
Fig. 4 Value of J in the upper external reservoir-fault corner as a function of scaled fault offset and fault dip.This value of J denotes the magnitude of the spatially dependent deviation from the average pressure.Values of  > 90 • indicate a reverse fault

Figure 7
Figure 7 depicts the scaled incremental pressure p p and the scaled incremental Coulomb stress C C acting on the normal fault ( = 70 • ) for various values of the dimensionless

Fig. 6
Fig. 6 Dimensionless additional resistance R as a function of fault dip and scaled fault offset t for a fixed value ŵ = 20

Fig. 7
Fig. 7 Scaled pressure and stresses under depletion with a fault throw t = 1 3 and dip = 70 • .The dotted lines indicate the coordinates of the reservoir corners and the gray shaded areas represent the reservoir geometry.Left: the dimensionless pressure p p along the fault for various values of the dimensionless pressure drops p .Right: the scaled incremental Coulomb stresses for various values of the dimensionless pressure drop p .The incremental Coulomb stresses have been scaled by a factor C = (1−2 ) |p|

Fig. 8
Fig. 8 Scaled pressure and stresses under injection with a fault throw t = 1 3 and dip = 70 • .The dotted lines indicate the coordinates of the reservoir corners and the gray shaded areas represent the reservoir geometry.Left: the dimensionless pressure p p along the fault for various values of the dimensionless pressure drop p .Right: the scaled incremental Coulomb stresses for various values of the dimensionless pressure drop p .The incremental Coulomb stresses have been scaled by a factor C = (1−2 ) |p| 2 (1− ) .The green patches indicate the regions of positive incremental Coulomb stresses for p = −1 where the fault moves closer to failure, while the cyan patches indicate the regions of positive incremental Coulomb stresses for p = −10 (B24) xx = D x ∬ Ω p( , )g x (x, y, , ) d d , (B25) yy = D y ∬ Ω p( , )g y (x, y, , ) d d , )g y (x, y, , ) d d + y ∬ p( , )g x (x, y, , )d d , (B27) e ij = ij − * ij (B28) xx =( + 2G) xx + yy − p Ω , (B29) yy = xx + ( + 2G) yy − p Ω , (B30) xy =2G xy .

Fig.
Fig. 12Example of a reservoir discretized into rectangular and triangular grid cells.The pressure may vary between grid cells, but the pressure inside each grid cell is assumed to be constant, such that closed-form expressions for the induced stresses may be used