A scalable collocated finite volume scheme for simulation of induced fault slip

We present a scalable collocated Finite Volume Method (FVM) to simulate induced seismicity as a result of pore pressure changes. A discrete system is obtained based on a fully-implicit fully-coupled description of ﬂow, elastic deformation, and contact mechanics at fault surfaces on a ﬂexible unstructured mesh. The cell-centered collocated scheme leads to a convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a generic multi-point ﬂux approximation is formulated to treat heterogeneity, anisotropy, and cross-derivative terms for both ﬂow and mechanics equations. The resulting system, though ﬂexible and accurate, can lead to excessive computational costs for ﬁeld-relevant applications. To resolve this limitation, a scalable processing algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate and eﬃcient; and provides a promising framework for accurate and eﬃcient simulation of induced seismicity in various geoscientiﬁc applications. © 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the


Introduction
Quantification of the mechanical deformation of subsurface geological reservoirs as a result of pore pressure change plays a key role in safe and optimal operations of many geoscience and geo-energy applications [1].As an example, rock mechanics stays at the core of hydraulic fracturing operation and design.Gas production, as another example, often causes subsidence and can initiate induced seismicity and serious damage to surface infrastructures [2][3][4][5][6].Moreover, in geothermal operations, the re-injection of cooler fluid causes stress and strain changes that can potentially activate faults.Successful exploitation of these geo-energy resources, therefore, depends highly on the development of robust and efficient computational simulators for coupled hydro-mechanical processes.
Several studies have been developed in the literature for simulation of geomechanics and induced seismicity.More precisely, Biot's consolidation model has been simulated with a finite difference scheme on staggered 1D [7] and 2D [8] grids.Galerkin finite element methods have been also employed by many researchers to investigate induced seismicity in different geoscience applications [4,9,10].Mixed finite elements methods were also extensively developed to allow for more accurate treatment of geomechanics [11,12].To enhance convergence properties of the underlying linear systems, weaklyimposed symmetry of the stress tensor in mixed finite element methods has been also developed [13] and employed in combination with a multi-point stress approximation [14].
Recently, the Finite Volume Method (FVM) has gained considerable interest in the computational geoscience community.This is despite the fact that a locally conservative stress field for geomechanics does not seem to be as crucial as the locally conservative mass flux typically required in multiphase flow simulations.Nonetheless, the FVM is still an attractive choice because it represents an integral form of conservation laws.Recent literature includes the development of the FVM for geomechanical simulations with both staggered [15,16] and collocated grids [17][18][19][20].The resulting linear systems are non-symmetric, yet can be made weakly symmetric [21].
Different approaches have been developed to address the approximation of fluid flow in faults.The Discrete Fracture Model (DFM) [22] allows faults to be resolved explicitly at the interfaces of the grid cells.It respects fault geometry, material contrasts in the vicinity of the fault.In contrast to DFM, the Embedded Discrete Fracture Model (EDFM) [23] and projectionbased EDFM [24] do not require conformal meshing of faults, hence, independent grids for faults and rock matrix can be employed.A comprehensive review and benchmarking study have been presented in reference [25].Both concepts have been successfully applied for mechanics and poromechanics modeling [9,10,26,27,18,28].In the embedded models special discontinuity basis functions are used to resolve contact mechanics [27,28] whereas the use of vertex-based [9,10] grids allows for natural treatment of the discontinuities within a DFM unstructured mesh.The present study follows the latter approach, where the displacement discontinuities are modeled as extra unknowns within a cell-centered collocated FVM grid.
Some authors [27,29] use fixed-stress splitting algorithms [30,31] to decouple mechanics and flow equations.These are a form of sequential implicit (SI) solution schemes and often lead to more efficient simulations than fully implicit (FI) simulations.However, sequential schemes introduce certain restrictions on time step sizes.On the other hand, FI schemes [16,9,29,10,18] are unconditionally stable and are often more robust and convenient approaches for the investigation of complex multiphysics problems.When comparing FI and SI approaches designed using the most efficient numerical techniques, the FI often outperforms the SI approach in the context of coupled thermo-compositional-mechanics simulation [29].
Although the FI approach does not imply any restriction on time step size, it requires efficient nonlinear and linear solution strategies for high-resolution models.One such strategy is to construct a preconditioner based on the idea of the SI approach.In [32], the authors employ a fixed-stress splitting concept in a sparse approximation of the Schur complement in order to obtain a block-preconditioned solution strategy.Later this approach was combined with a constrained pressure residual (CPR) preconditioner to construct a robust and effective solution strategy for coupled multiphase flow and mechanics [33].
In this study, we develop a collocated FI multi-point FVM scheme for poromechanics simulation of faulted reservoirs, extending the work of [19,20].The scheme can be used to solve arbitrarily anisotropic poromechanics problems on unstructured polyhedral grids with a minimum number of degrees of freedom per cell.It is also capable to take into account material heterogeneity while preserving mass and momentum balances.To this end, we also extend the state-of-the-art collocated FVM to take into account discontinuities in displacements at faults consistent with the DFM framework.More precisely, in contrast to the state-of-the-art existing collocated FVM schemes [18,26], in this work, the displacement discontinuities are considered as primary unknowns.This allows for the convenient development of a coupled multi-point approximation with contact interfaces.Such an approach simplifies the incorporation of friction laws and the general formulation of contact mechanics.To further enhance the linear scalability of our method, a block-partitioned preconditioning strategy is also implemented.
The developed methods are implemented in the open-source Delft Advanced Research Terra Simulator (DARTS).DARTS is a scalable parallel simulator, which has been successfully applied for modeling of hydrocarbon [34,35], geothermal [36][37][38] and CO 2 sequestration [39,40] applications.This study adds geomechanical modeling to the existing advanced hydrothermal capacity of DARTS, making it a fully hydro-thermo-mechanical simulator for complex geoscience and geo-energy applications.
The rest of the paper is organized as follows.First, we present the governing equations for coupled flow and geomechanics.Then we develop the discrete system of equations, the coupled multi-point approximation, and the strategy to solve this system.Next, we perform validation against analytical and numerical approaches.Thereafter, we study slip and stress profiles over a fault in different configurations.Through these proof-of-concept test cases, we assess the performance and applicability of the proposed simulation approach.

Governing equations
Single-phase flow in poroelastic media is governed by a coupled system of stationary momentum balance and fluid mass balance equations [41], i.e., subject to the relations as defined in [41,42]: where is the rank-two total stress tensor, h is a vector of volumetric forces, φ is porosity, ρ is fluid density, v is Darcy's velocity, h is a source (sink) of fluid mass, C is the rank-four drained stiffness tensor, E is the rank-two strain tensor, u = [u x u y u z ] T is a vector of displacements, ∇u and (∇u) T stand for the Jacobian of the displacements and its transpose respectively, B is the rank-two Biot (tangent) tensor [41][42][43], p is pore pressure, β = (B : I) / 3 is one-third of the trace of B, I denotes the identity matrix, K s is the bulk modulus of the solid phase, K is the rank-two permeability tensor, μ is fluid viscosity, and g is acceleration of gravity.Throughout the paper we follow the sign convention that tension stresses are positive.Furthermore, a subscript "0" denotes the reference state of a variable, i.e., Eqs. ( 2) and ( 3) represent stress and porosity changes within anisotropic poroelasticity, Eq. ( 4) describes infinitesimal strains, and Eq. ( 5) is Darcy's equation.The fluid properties density and viscosity are functions of pressure, as stated in Eq. ( 6).We note that the full Biot tensor results from general thermodynamic assumptions to poroelasticity but that in most practical applications it is only possible to experimentally determine a single scalar Biot coefficient.However, we maintain the full Biot tensor formulation in view of its potential relevance to upscaling.All variables have been listed in the Nomenclature section at the end of the paper.We denote vectors with bold lowercase letters, second-order tensors and matrices with bold capital letters, and the tensors of rank higher than two with script font.We also define the total traction vector f over an interface with unit normal n as Boundary conditions for the system in Eq. ( 1) can be written as where f b , u b and p b are traction, displacement and pore pressure at the boundary, respectively.In addition, a n , b n , a t , b t , a p , b p are coefficients that determine the magnitude of their corresponding boundary conditions, while r n , r t , r p are the corresponding condition values.Eqs. ( 9) can describe a broad range of possible boundary conditions including The rank-two Terzaghi effective stress tensor [41] and effective traction vector f are defined as At the fault interfaces we consider a gap vector g that is equal to the jump of displacements over the contact g = u + − u − , where +/− denote a particular side of the fault.The conditions that are applied to the contact read ġN = 0, (11) ġT = 0, < 0, (stick), (13) where f N = n T f and f T = (I − nn T )f are the scalar normal and vectorial tangential projections of f on the fault; g N and g T are the equivalent normal and tangential projections of g on the fault; ġ stands for the time derivative of the gap vector; η is the friction coefficient, and = f T − ηf N is the Coulomb friction function.Eq. ( 11) represents a non-penetration condition which prohibits opening as well.This implies that the gap vector only allows for tangential displacements and therefore becomes a slip vector.Eq. ( 12) governs relaxation of tangential traction once slip occurs, and Eq. ( 13) sets the change of tangential gap (i.e. the slip) to zero if the slip criterion is not exceeded.

Numerical scheme
In this section, we briefly describe the discrete schemes used for numerical solution of the governing equations.First, we present FV schemes for the equations in matrix and fault cells.The next parts provide the derivation of flux approximations at both continuous cell interfaces and faults.

Discrete equations for matrix
According to the FVM, the momentum balance and fluid mass balance in Eq. ( 1) in integral form for the i-th cell can be stated as where superscripts n, and n + 1 denote the variables taken from the current and next time step, respectively.Also, subscript j enumerates the interfaces of the i-th cell, t is the time step size, V i is the volume of i-th cell, δ ij denotes the area of the interface between cells i and j, f is defined in Eq. ( 8), and φ, q and qij , are defined as where the following relations are used Here B is assumed to be symmetric, porosity in Eq. ( 3) is treated as φ = φ + B : ∇(u − u 0 ) and the last term is approximated using Gauss' formula as a sum of fluxes qij over cell interfaces.The term (ρ/μ) ij is calculated using an upwind approximation while the approximations of f, q and q are discussed in the following sections.

Discrete equations on faults
To satisfy Eqs. ( 11)-( 13), we use a penalty regularization [44,9,45] which leads to the return-mapping algorithm of g n+1,k+1 N where indices n and k denote time step and nonlinear iteration respectively.Moreover, f T denotes trial traction, which represents the penalized effective tangential traction [44], the Coulomb friction function used as a slipping criterion is evaluated at the trial state = ( f ) that accounts for the change of slip g T over the time step.Macaulay brackets a are equal to a if a ≥ 0 and otherwise equal to zero.Thus, in the slip state = 0 Eq. ( 19) requires contact to remain at the slipping surface defined by = 0 where the direction of forces is defined by the trial traction.Contact reaches the stick state once the slip increment in Eq. (18) becomes negligible compared to the previous traction ( ġT = 0).In this case Eq. (19) claims the traction to be equal to the trial one.Although the model is capable to take into account a friction coefficient η as governed by various friction laws, we only consider a constant friction coefficient in this paper.In our experience, the return-mapping algorithm described in Eqs. ( 18)- (19) does not exhibit significant convergence problems, except the cases with severe inf-sup instability (pressure oscillations) and when the slip direction inverses.We may also expect convergence issues in the presence of intersecting faults or in the case of hydraulically active fault when its volume and transmissibilities depend on the aperture g N .
Although Eqs. ( 17)-( 19) can be treated as boundary conditions with no extra degrees of freedom [26], it might be convenient to treat the gap vector g as an unknown assigned to particular fault cells, especially in induced seismicity applications.Moreover, it allows for maintaining the block structure of the Jacobian.
The fluid mass balance in fault cells can be written as where (ρ/μ) ij is taken from the upwind direction, fluxes of displacements qn ij , qn+1 ij are evaluated only between fault and matrix cells, and where the approximations of f, q, and q are discussed in the following sections.It is worth to be pointed out that we expect discontinuous displacements over the fault interface whereas pressure remains continuous there according to the assumptions of the DFM approach.
The influence of mechanical stresses on the conductivity of hydraulically active faults can be of high interest in the modeling of hydraulic fracturing.Although we do not consider fault opening in this paper, this effect was addressed by other researchers using FEM with penalty regularization [9], Nitsche method [10] or Lagrange multipliers [46]; and in the FVM framework with using Lagrange multipliers [18].The latter development within FVM also demonstrates the applicability of our developed method to model the hydraulically active faults.
Below we will use the following decomposition where subscripts i = 1, 2 refer to the cells neighboring an interface, r 1 and r 2 are distances between the cell centers and the interface, y 1 , y 2 are projections of the cell centers on the interface, T 1 and T 2 are 3 × 3 matrices, 1 and 2 are 3 × 9 matrices, while scalars λ 1 and λ 2 and vectors γ γ γ 1 and γ γ γ 2 provide co-normal decompositions of K 1 and K 2 .The 4 × 1 vectors ξ 1 and ξ 2 , and the 12 × 1 vectors ξ τ 1 and ξ τ 2 represent normal and tangential projections respectively of the gradients of the unknowns.
Using the notation introduced above, the continuity of fluxes, represented in Eqs. ( 22) and ( 23), can be written in the coupled form of where Here A i , Âi and Q i are 4 × 4 matrices, i is a 4 × 12 matrix and ζ i is a 4 × 1 vector and "±" in According to Eq. ( 21) the tangential projections of the gradients are ξ τ 1 = ξ τ 2 = ξ τ .Deriving ξ 2 from Eq. ( 21) and substituting the result into Eq.( 29) we obtain the following expression for ξ 1 : Substituting Eq. ( 33) into the left-hand side of Eq. ( 29), one obtains the multi-point approximation for the traction f, as given in Eq. ( 8), and the total fluid flux q t = q + (q n+1 − qn )/ t as The coupled multi-point approximation Eq. ( 34) was presented in [20].Approximation of the boundary conditions defined in Eq. ( 9) requires a specific treatment and we refer to [20] for the details.The approximate Eq. ( 34) involves fluid properties which may depend on pore pressure or composition in the case of multicomponent flow.In such cases Eq. (34) can not be applied directly and further steps are required, as will be discussed below.

Gradient reconstruction
The approximation in Eq. ( 34) requires the tangential projection of the gradients of the unknowns to be reconstructed.
One can derive ξ 2 from Eq. ( 29) and substitute it into Eq.( 21) to obtain the interpolation equation as (35) It is necessary to consider at least 3 interfaces (in 3D) of the first cell to enclose the system with respect to the 12 components of ∇ ⊗ w 1 .
Bringing together the results of Eq. (35) for N considered interfaces of the i-th cell, we build up the system where M i is a 4N × 12 matrix and D i is a 4N × 4(N + 1) matrix of coefficients in front of the corresponding unknowns at the right-hand side of Eq. ( 35), while ψ i is a 4(N + 1) × 1 vector of N + 1 unknowns (or free terms in right-hand side of boundary conditions).The solution of Eq. ( 36) can be obtained in a least-squares sense as Fig. 1.Cells that contribute to the approximation of fluxes over the interfaces of cell i.Index j denotes the nearest neighbors of cell i.Index k denotes farther neighbors that contribute to the gradients reconstructed in cells j.
For the approximation of ξ τ in Eq. ( 34) the following combination of gradients is used in [20] A set of cells that contribute to the approximation Eq. ( 38) for each interface of some cell i is illustrated in Fig. 1.
It is worth to be mentioned that the least squares solution in Eq. ( 37) allows computing the gradients of unknowns locally and independently for every cell.Note, however, that it does not guarantee the local conservation property for the scheme.In order to maintain it, individual gradients for every interface that respects the corresponding flux balance should be employed.

Introducing discontinuity
DFM is a well-known technique of fault representation extensively used for the description of flow dynamics in fractured subsurface reservoirs [22,25].In combination with contact mechanics, it has been proved for practical engineering subsurface problems [45,10] and was widely validated against semi-analytical poromechanics solutions [9,18].
Although faults are usually represented as 2D surfaces in computational grids, some aperture is assigned so that they have a volume greater than zero in the discrete fluid mass balance.Displacements can be discontinuous at faults while pressure remains continuous according to the assumptions of the DFM (see Fig. 2).In that case Eq. ( 21) does not hold anymore and the tangential gradients of the displacements may be different at different sides of the fault.We introduce an additional degree of freedom (d.o.f.) per faulted interface, namely gap vector g = u + 2 − u − 1 , which is equal to the jump of displacements over the interface.The following continuity condition will be used instead of the top three rows in Eq. ( 21) (39) where "±" is positive when the first and second cells are located at negative and positive sides of the fault respectively.Assuming a zero normal pressure gradient in a fault cell, the continuity of momentum fluxes, as given in the top three rows in Eq. ( 29), can be rewritten as where we redefine g ← {g, 0}.Eq. ( 33) changes to (41) and Eq. ( 34) becomes Also, the following expression can be used for the reconstruction of displacement gradients (top 9 rows) instead of Eq. ( 35) (43) where gap gradients are reconstructed by simply using We use a Discrete Fracture Model (DFM) approach [9] for the flow which implies having discrete equations for fracture segments of non-zero volume.The particular scheme was shown in Eq. (20).We assume continuous pressure in Eq. ( 21), and we have two matrix-fracture connections per fracture segment in the fluid mass balance for that segment.The following two-point approximation of Darcy flow is used for these connections where ρ 12 = (r 1 ρ 2 +r 2 ρ 1 )/(r 1 +r 2 ) is linearly interpolated between cells.We assume that the flow in fractures is dominated by fracture conductivity and use a two-point flux approximation Eq. (45) for the fluid fluxes between fracture segments.
Since we consider continuity conditions for fluid mass balance in Eq. ( 22) and for pressure in Eq. ( 21) between matrixfracture and fracture-fracture cells, we use the following equation in the reconstruction of pressure gradients for these connections Although Eq. ( 42) approximates the total traction f at contact, we found that taking pressure p δ from fault cell in the approximation of the Biot term p δ Bn leads to more accurate results.This can be achieved assuming Â1 , Â2 = 0 in Eq. (42) and adding p δ B δ n in order to calculate total traction f, where B δ denotes matrix of Biot coefficients defined in a fracture cell.Note, that in this case neither the approximation of Biot term in total traction, nor two-point flow approximation in Eq. ( 45) depends on pressure gradients.However, the condition in Eq. ( 46) is taken into account in other fluxes of the contacting cells.We expect that approximation of fluxes over the adjacent faces that belong to different fractures (in the case of crossing fractures) should not bring additional difficulties.However, this aspect has not been tested yet and the consideration of complex fracture networks is outside the scope of this study.

Reconstruction of stresses
Solving the system of discrete equations Eqs.(14) provides the vector of unknowns in the cell centers.To reconstruct effective stresses at the same locations we use the algorithm described in [20].For the jth interface of the cell let us construct the matrices where n x , n y and n z are components of the unit normal to the jth interface; x j , y j and z j are components of the jth interface center position; and x V , y V and z V are components of the cell center position.
Collecting these matrices for N interfaces of the cell we will obtain the following 3N × 6 and 3N × 6 matrices, and 3N × 1 vector , (48) where vector b represents the Biot tensor B in Voigt notation.The stresses at the cell center can then be reconstructed using the least squares solution where the Biot effective stress tensor [41] σ = [σ xx , σ yy , σ zz , σ yz , σ xz , σ xy ] T is written in Voigt notation.

Discretization
In real-world subsurface applications, the fluid is usually represented as a mixture composed of multiple chemical components and phases.In the case of multiphase multicomponent fluid flow, its governing properties like density, viscosity, and others depend on mole fractions, pressure, and temperature.Eqs. ( 42) and ( 43) involve fluid viscosity and density (in the gravity term) which means that they may depend at least on pressure.These approximations become nonlinear with respect to pressures and demand nonlinear iterations.It is therefore required to repeat the gradient reconstruction and the approximations at each iteration throughout the nonlinear solution.
Assuming that the coefficient of the Biot tensor remains constant, the terms that include these coefficients can be omitted in the local problem for poromechanics.This allows fluid density and viscosity to be excluded from these expressions and the approximations in Eq. ( 42) to be calculated once before iteration over time.Density and viscosity can be taken into account having the equations for fluid flow and flux of displacement in Eq. ( 15) already approximated.
This approach requires a separate assembly of Darcy, Biot, and gravity terms in Eq. ( 15).The Biot term does not affect the reconstruction of the gradients because the tensor of Biot coefficients is constant over the domain whereas the gradients are linear with respect to the gravity contribution.In turn, approximations of tractions and the terms in the fluid flow expression are linear with respect to gradients and gravity.The operation can be performed by omitting all left-bottom terms in the matrices in Eqs.(30), (31) and the Biot terms in vectors ζ 1 and ζ 2 and setting t = 1 in the other terms.This technique, as introduced in [20], allows us to retrieve the approximation of the Darcy flux q separately from the Biot contribution (q n+1 − qn )/ t, and to decouple fluid properties from approximation of fluxes.It also allows us to calculate the coefficients of the approximations once in the pre-processing stage.The flux of displacement q can then be calculated as (50) where Ã1 = A 1 − Â1 .At the boundaries, this flux can be calculated separately as well; see [20] for details.
The first stage of discretization consists of the reconstruction of gradients according to Eqs. ( 35), ( 43) and (44).Once the gradients are reconstructed, they can be used to approximate fluxes according to Eqs. ( 34), ( 42) and (50).For non-faulted interfaces, the averaged gradient is used as defined in Eq. (38).

Linear solvers
In order to improve the efficiency of the solution of the system given in Eq. ( 1), an advanced linear solver strategy is required.In [32], a fixed-stress split concept is used to construct a block-partitioned preconditioner for poromechanics, which was extended for coupled multiphase flow and mechanics in [33].Here we employ the first stage of this approach to construct an efficient preconditioner.
The idea of the method is to consider the block-partitioned linear system where J, r u and r p are the Jacobian and the residuals for momentum and mass balance equations produced by the numerical scheme, δu and δ p are unknown increments of the vector of displacements and pressure, J uu , J up , J pu , J pp are contributions to momentum balance and mass balance equations from displacement and pressure unknowns, L is the lower-triangular term in an LDU decomposition of the Jacobian, and S pp = J pp − J pu J −1 uu J up represents the Schur complement of block J uu in the Jacobian.
The concept of fixed-stress splitting is used here to provide the sparse approximation of S pp as where e = [1, 1, .., 1] T is a probing vector, and P −1 u is a preconditioner used for the elasticity system.This preconditioner is applied as follows: 1.At the beginning of every nonlinear (Newton) iteration, Spp is evaluated.A single V-cycle of an algebraic multi-grid (AMG) solver is typically used for P −1 u .
2. At every iteration of the linear solver (GMRES) we employ the approximation Spp to solve the upper triangular system in Eq. ( 51) as follows: (a) Solve the second equation in upper triangular system in Eq. ( 51): Spp δ p = −r p .(b) Use the known value of δ p to solve the first equation in Eq. ( 51): J uu δu = −r u − J up δ p. 3. Provide the linear solver with the approximate solution { δu, δ p} as initial guess.
The system for the pressure unknowns is preconditioned using a single V-cycle as well.The generalized minimal residual (GMRES) method is used as an outer solver.

Convergence study
We make a comparison with the analytical solution proposed in [19] to check the order of approximation on different grids.We consider a 3D unit cubic domain  at the right side inside the domain and apply Dirichlet boundary conditions according to the reference solution.
To assess the convergence of the numerical solutions we calculate the discrete L 2 norm of their deviation from the reference solution.Results are presented in Fig. 3 for calculations performed using a structured cubic mesh, an unstructured mesh of wedges (an x-y 2D triangle mesh extruded along z-axis), and an unstructured hexahedral mesh (an x-y 2D corner point grid extruded along the z-axis).The last two meshes are truly unstructured only in the x-y plane.The solution obtained on the structured cubic mesh demonstrates 2nd-order convergence with respect to displacements and more than 1st-order convergence with respect to stresses.On the unstructured wedge mesh, the solution also shows 2nd-order convergence with respect to displacements and nearly 1st-order convergence with respect to stresses.

Displaced fault
Faulted reservoirs can experience different mechanical conditions over their geological history.Nearly always this results in an offset such that stratigraphic layers become displaced.An analytical investigation of stresses initiated around a displaced fault in a depleted reservoir was carried out in [47].Here we compare the solution we obtained numerically with the analytical one.The model domain is illustrated in Fig. 4a.It consists of two regions which are treated as elastic bodies with Young's modulus E = 16 GPa, Poisson's ratio ν = 0.2, and Biot's coefficient b = 0.9.The reservoir region (in the middle) has permeability 100 mD whereas the permeability of the surrounding rock is equal to 10 −5 mD.
A comparison of the numerical solutions with the analytical one is displayed in Fig. 4b.A build-up of pore pressure by p = 10 MPa inside the reservoir (width 800 m, height 300 m and offset 100 m) perturbs the system.The calculated profile of the shear component of the stress tensor over the vertical centerline of the reservoir (white line) is compared against the analytical result.The calculated solution matches the analytical one quite well.

Contact problem
The next case concerns a fault with Coulomb friction in the center of a square domain of size a = 1 m; see Fig.    Eq. ( 12) is exceeded and the fault is prohibited from opening as stated in Eq. (11).A structured quadrilateral grid is used in the calculation.The resulting displacement and stress fields for the vertical fault are shown in Figs. 6 and 7.One can notice the small jump in horizontal displacements across the fault.
The results are compared with other ones, obtained by the PorePy simulation tool [48], for different orientations of the fault determined by a dip angle φ.The comparison is displayed in Fig. 8. Tangential tractions over the fault calculated by DARTS and PorePy match quite well.For higher dip angles (nearly vertical orientation) the resulting slip also fits quite well.
However, for lower dip angles the slip is decreasing and in the case of φ = 72 • it becomes located at the tips of the fault whereas the center of the fault displays higher normal tractions and smaller slip values.With decreasing dip angles, the slip magnitude decreases as well and a tiny mismatch in tangential traction results in a higher mismatch in slip (for φ = 72 • ).
Below some threshold angle, no slip over the fault is observed.

Accuracy on adaptive meshes
The figures displayed in the previous subsection were calculated using a structured hexahedral grid that provided stable solutions.However, an attempt to calculate contact problems using adaptive unstructured grids may fail to produce smooth traction profiles.We considered a no-slip case within the same setup.This can be achieved either by increasing the friction coefficient or by using a meshed line without displacement discontinuity.To avoid the complexity of extra terms associated with a contact, we follow the last option for a line orientation φ = 90 • (see Fig. 5).The resulting traction profiles obtained using wedge-type unstructured grids of different adaptivity are compared against the calculations made using a structured hexahedral mesh which is treated as a reference.The cell sizes near the fault and the border define the mesh adaptivity.
Tangential and normal components of traction over the fault are presented in Fig. 9.The three adaptive meshes used in calculations are depicted at the right of the figure.They have cell sizes [0.003, 0.1] , [0.003, 0.03] , [0.03, 0.03] (in meters) near fault and border respectively.The domain size is a = 1 m and the fault length is L = 0.4a.In this comparison, it is seen that the calculation made using the mesh of large cell size at the boundary of the domain (red lines) provides a poor approximation of the traction vector over the fault.At the same time, the less accurate near-the-fault mesh (purple lines) results in more accurate profiles even though it involves a lower number of cells.One can observe that the profiles of normal traction exhibit some noise and are slightly skewed whereas the tangential ones remain smooth and symmetric.This may be explained by the absolute values they have.The boundary conditions applied in this example result in a smaller normal traction than the tangential one.That means that the accuracy of calculated normal traction is poor because the whole traction vector lies almost within the plane of the contact interface.The relative tolerance of traction components is important to resolve the friction contact.

Linearity preservation
with the following linear reference solution  which satisfies the system in Eq. ( 1).Note that in this example we employ a full Biot tensor B. Although such a matrix is of no physical significance for practical applications, it serves to test computational aspects of the algorithm.We assume that density and viscosity in Eq. ( 6) remain constant and that the porosity is equal to where M is Biot's modulus, while p 0 and u 0 are fluid pressure and displacement vector at reference state.Substituting Eqs. ( 55), ( 54) and (56) into Eq.( 1), we can derive the right-hand side terms.Using these terms while specifying initial and boundary (Dirichlet) conditions according to Eq. ( 55), we expect to observe the reference solution inside the cube.The parameters used for calculation are listed in Table 1.
The results are summarized in Table 2.The calculations were performed using structured grids composed of cubes and wedges, and an unstructured tetrahedral grid up to time t = 1 day.The linear solution was observed for every time step and for each type of grid.Absolute errors at t = 1 day with respect to the reference solution in Eq. ( 55) are almost at the level of machine double floating-point precision.

Convergence
For the same and the following reference solution [20]: The source terms h and h are calculated by substituting Eqs.(57), (56) and (58) into Eq.(1) using automatic differentiation.Dirichlet boundary and initial conditions are applied according to Eq. (58).We assume that density and viscosity in Eq. ( 6) remain constant while porosity is calculated according to Eq. (56).All relevant parameters used in the calculations are listed in Table 3.
L2 norms of differences between reference and calculated solutions are depicted in Fig. 10.The calculations were done for a few cubic meshes of different resolutions and using different time step sizes.They demonstrate superlinear convergence for the vector of displacements and nearly linear convergence for the pressure.The fact that the order of convergence is below the second order can be explained by the Backward-Euler scheme used for the time integration [20].

Mandel's problem
Consider a rectangular domain illustrated in Fig. 11.Roller boundary conditions are applied to the left and bottom boundaries of the domain.The right boundary is free of both normal and tangential forces while a normal load is applied from the top.No-flow conditions are specified for all boundaries except for the right one where a Dirichlet condition p = p 0 is applied.This set-up is the so-called Mandel problem which is often used as an example to demonstrate specific aspects of poroelasticity [41].
We consider a porous homogeneous domain characterized by Young's modulus E = 1 GPa, Poisson's ratio ν = 0.25, a diagonal permeability tensor k x = k y = k z = 1 mD, saturated with a single-phase fluid with viscosity μ = 0.0981 cP, and with a Biot modulus M = 5.44 • 10 −6 bar −1 and a Biot tensor equal to the identity matrix B = I.The distributed load f = 10 MPa is applied to the initially undeformed domain with p 0 = 0. Fig. 12 depicts a comparison of pressures for the numerical and analytical solutions.Time is expressed in dimensionless [49,50], where S is the storage coefficient [50], k = k x is the scalar permeability, L is the horizontal extent of the domain, K d is the drained bulk modulus, G is the drained shear modulus, b = 1 is the scalar Biot coefficient, M is the Biot modulus, and μ is fluid viscosity.Pressure is normalized to half the applied load.The calculated pressure distribution over the entire domain is depicted in Fig. 13.
The efficiency of the fixed-stress block-partitioned preconditioning strategy is tested in the Mandel problem against ILU0 preconditioner.The results are presented in Table 4.The number of linear iterations (LI) made by the block-partitioned preconditioner is more than a magnitude lower than the number of iterations performed with ILU0 preconditioner.At the same time, we observe only a double reduction in the time the solver takes for the setup and solve calls.We attribute this to the cost of AMG setup for J uu block.

Core injection experiment
The developed framework has been applied to model fluid-induced fault slip behavior observed in a core injection experiment [51].The authors performed two fluid-induced slip experiments (SC1, SC2) conducted on permeable Bentheim sandstone samples crosscut by a fault.Fault slip was triggered by fluid injection into the core at different rates.As a result, they measured slip, slip velocity, normal and shear stresses, and the apparent friction coefficient as a function of time.
In order to approximately simulate the experiment the following setup is considered.A rectangular 2D domain (Fig. 14) of 50 × 100 mm size is loaded from the left and the right with a confining normal stress |f| = 350 bar and is loaded from The normal displacements applied from the top are chosen to match the initial vertical stress (σ 1 in [51]) observed in the experiment.Using a Neumann condition for the top boundary directly is not possible because once the fault starts sliding, the problem becomes purely Neumann so that the displacements in the top piece become defined with respect to a constant.
Water of constant viscosity μ = 10 −3 Pa • s and compressibility 10 −9 Pa −1 is injected into the domain through the bottom boundary at a specified pressure p = p(t) whereas the other boundaries are impermeable.Evolution of injection pressure and friction coefficient are taken precisely from the experimental data and are shown in Fig. 15.The injection pressure is increased step-wise in SC1 while in SC2 it is increased gradually while keeping the pressure constant during short periods of time.
Fig. 15 depicts the dynamics of experimental slip, slip velocity, normal and shear stresses compared to calculated values.Although we used a 2D approximation of a 3D setup, there is a good match in terms of slip, both stresses, and slip velocity in both tests.It is worth to be pointed out that we used the friction coefficient calculated from the experiment; interpretation of friction behavior according to friction laws is out of the scope of this paper.

SPE 10 with mechanical extension
In this field-scale test case, we use data from the SPE10 benchmark for flow supplemented by mechanical parameters, in particular spatial distribution of Young's modulus that depends linearly on the porosity [29].The original dataset represents a reservoir characterized by a channelized permeability and by a permeability field that has a Gaussian spatial covariance as shown in Fig. 16.The dataset was coarsened using a volume-averaging approach [29].Although the original SPE10 benchmark was designed as a two-phase flow problem, here we consider single-phase flow.The reservoir is produced by a single doublet of an injector and producer.No-flow boundary conditions are prescribed for all the boundaries.Normal displacements and tangential tractions are set to zero at all boundaries except for the top boundary where a uniform distributed load of 900 bar is applied.Poisson's ratio is taken as constant ν = 0.2, and Young's modulus and lateral permeability fields are depicted in Fig. 16.  17.They demonstrate the applicability of block-partitioned preconditioning for the solution of a discrete system produced by a coupled FVM multi-point scheme.
In Fig. 17a it can be seen that the number of linear iterations (LI) increases as the problem size increases while the number of nonlinear iterations (NI) remains equal to 20.Similar behavior was observed in the nonlinear problem of twophase flow in a deformable poroelastic medium [33] and was explained by the effect of upscaling.The upscaled grid tends to be more homogeneous and less anisotropic, which may also improve the convergence of the linear solver.To distinct the    effect of upscaling, we also present calculations for a homogeneous reservoir in Fig. 17b).In the homogeneous case, the number of linear iterations flattens as the problem size grows which demonstrates the scalability of the solution strategy.

LBB-instability
In the original test case reported by [20], the FV approximation showed a divergent solution for very small time steps.In our implementation, we observed unstable behavior for small time steps as well.In Fig. 18 the evolution of the discrepancy between true and computed solutions is depicted.The setup was taken from test case 1 and was calculated using an 8 ×8 ×8 cubic grid for a time period t ∈ [0, 0.01] days.For calculations made with t = 0.01 and t = 0.001 the discrepancy remains at the level of machine precision whereas for smaller time steps it rises drastically.
Such behavior is explained by the violation of the inf-sup stability condition which is also known as Ladyzhenskaya-Babuska-Brezzi (LBB) condition which is a sufficient condition for a saddle point problem to have a unique solution [52].
The saddle point appears in our numerical scheme because matrices A i and A j defined in Eq. ( 30) have eigenvalues of different signs.The occurrence of the saddle point can be avoided using the eigen-splitting of the indefinite matrices [26] which we will address in our future work.

Conclusion
We have developed a collocated Finite Volume Method with a multi-point approximation of fluxes for geomechanics and poromechanics.A collocated grid allows assigning unknowns and material properties in the same regions such that the block structure of the Jacobian is maintained.A fully coupled approximation provides a robust and unified framework to approach hydromechanical problems.The method can cope with discontinuities in displacements, as occur in faults, on the level of discretization.The multi-point approach requires to perform discretization only once in the pre-processing stage.These discrete terms remain valid for all different contact states.We introduced gap degrees of freedom over the fault, which significantly simplified the formulation of contact conditions using a minimal number of degrees of freedom.We validated the method against several analytical and numerical solutions in different test cases.In order to speed up the linear solution, we used a block-preconditioning solution strategy based on the fixed-stress splitting concept in the approximation of a sparse Schur complement.Based on the literature [20] we expect that our method can be extended to more complex grid geometries.The developed numerical framework demonstrates convergence with timestep and grid refinement.The multi-point stress discretization behaves robustly on either structured or unstructured simplex and hexahedral grids.The developed solution strategy is capable of providing efficient solutions even for highly resolved models at different scales of interest.In the limit of an incompressible or impermeable medium or at small timesteps, the coupled collocated scheme may exhibit unstable behavior, which can be explained by the occurrence of a saddle-point problem.

Numerical variables
A i , Âi -4 × 4 matrices, γ γ γ i -tangential projections of co-normal permeability vector, i -3 × 9 matrices, t -time step size, δ j -area of jth interface, ζ i -4 × 1 vector, i -4 × 12 matrix, λ i -tangential projection of co-normal permeability vector, J -Jacobian matrix, J uu , J up , J pu , J pp -blocks of globally-partitioned Jacobian matrix, ξ i , ξ τ i -4 × 1 and 12 × 1 vectors of normal and tangential projections of the gradients of unknowns in ith cell, P −1 u , P −1 p -preconditioners used for displacement and pressure blocks of Jacobian matrix, Q i -4 × 4 matrices, r i -distance between ith cell center and interface, S i -9 × 9 matrix, S pp -Schur's complement of J uu in Jacobian matrix,

Fig. 2 .
Fig. 2. Fault representations in the approximation of momentum fluxes (a) and the approximation of fluid mass fluxes (b).

Fig. 3 .
Fig. 3. L 2 norms of errors in displacement (a) and stress tensors (b) computed on structured cubic, unstructured extruded wedge and unstructured extruded hexahedral meshes.The order of convergence for each mesh mentioned in the legend.

Fig. 4 .Fig. 5 .
Fig. 4. The domain with a vertical displaced fault and boundary conditions (a).Only inside the reservoir, the pore pressure was subtracted from the stress tensor.Tangential stress profile over the vertical centerline calculated using analytics and numerical simulation (b).(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) 5. The right boundary is fixed, the top and bottom boundaries are free of any forces, while displacements u le f t = {0.001,0.01} m are prescribed at the left.A plane strain setup is considered.The stiffness matrix is determined by Lame coefficients λ = G = 1 Pa.A fault of length L = 0.4a with friction coefficient η = 0.85 is allowed to slip once the Coulomb criterion in

Fig. 6 .
Fig. 6.Horizontal (left) and vertical (right) displacements for the case of a vertical fault (φ = 90 • ).Note the jump in vertical displacements at the fault.

Fig. 8 .
Fig. 8.Comparison of slip and tangential traction along the fault for different dip angles.

Fig. 9 .
Fig.9.Tangential and normal components of traction vector over the fault.Calculations were made using three adaptive wedge-type meshes depicted at the bottom and defined by cell sizes near fault and border.The reference profiles are calculated using a structured hexahedral grid.

Fig. 10 .
Fig. 10.Error norms of pressure and displacements calculated using cubic meshes of various resolutions.

Fig. 11 .
Fig. 11.The domain and boundary conditions for Mandel problem.

Fig. 12 .
Fig. 12.Comparison of calculated pressure (DARTS) with analytical solution.Normalized pressure at the left bottom cell over time (left) and normalized pressure along horizontal centerline at different moments of time.

Fig. 13 .
Fig. 13.Pressure over domain in bars at various moments of time.

Fig. 15 .
Fig. 15.Numerical simulation of core injection experiments from[51].Two injection scenarios were considered: step-wise increase (SC1) and gradual increase (SC2) of injection pressure from 50 to 350 bar.Time-dependent friction coefficient was taken from experiment.

Fig. 16 .
Fig. 16.Young modulus (in GPa) and lateral permeability (in mD) fields shown from the top (top row) and from the bottom (bottom row) of reservoir.Young's modulus is calculated as a linear function of porosity.The top 24 layers represent a channelized structure whereas the bottom 16 layers correspond to a Gaussian distribution.

Fig. 17 .
Fig. 17.Cumulative number of linear iterations and computational time taken by the linear solver to calculate 20 time steps and 20 nonlinear iterations for five different resolutions.For every resolution results are depicted for an upscaled heterogeneous reservoir (a) and for a homogeneous one (b).The rise in the number of linear iterations in the heterogeneous case can be explained by the effects of upscaling, while the flattening in the homogeneous case demonstrates the scalability of the solution strategy.

Fig. 18 .
Fig. 18.L2 norm of the difference between computed and true solutions over time.The setup is the taken from Test Case 1. Calculations were made with different time steps with magnitude t = 0.01, 0.001, 0.0001 and 0.00001 d.Below a certain time step size the error increases sharply.
9 × 9 matrices where C denotes a 6 × 6 symmetric stiffness matrix in Voigt notation and where

Table 1
Parameters used for calculations in linearity preservation test case.

Table 2
Absolute errors of solutions calculated on different mesh types for Test Case 1.

Table 3
Parameters used for calculations in convergence test case.

Table 4
Speed up obtained using fixed-stress preconditioning compared with an ILU0 preconditioner for the full system.TS -time steps, NI -nonlinear iterations, LI -linear iterations.