Accuracy of fully coupled and sequential approaches for modeling hydro- and geomechanical processes

Subsurface flow and geomechanics are often modeled with sequential approaches. This can be computationally beneficial compared with fully coupled schemes, while it requires usually compromises in numerical accuracy, at least when the sequential scheme is non-iterative. We discuss the influence of the choice of scheme on the numerical accuracy and the expected computational effort based on a comparison of a fully coupled scheme, a scheme employing a one-way coupling, and an iterative scheme using a fixed-stress split for two subsurface injection scenarios. All these schemes were implemented in the numerical simulator DuMux. This study identifies conditions of problem settings where differences due to the choice of the model approach are as important as differences in geologic features. It is shown that in particular transient and multiphase flow, effects can be causing significant deviations between non-iterative and iterative sequential schemes, which might be in the same order of magnitude as geologic uncertainty. An iterated fixed-stress split has the same numerical accuracy as a fully coupled scheme but only for a certain number of iterations which might use up the computational advantage of solving two smaller systems of equations rather than a big monolithical one.


Introduction
Fluid injection into the subsurface is common to different scenarios, often related to the world's energy consumption which relies heavily on the utilization of the subsurface. Coal, oil, and natural gas account for 81-85% of the world's primary energy consumption in 2017 [13,25]. Often, the extraction of these raw materials requires fluid injections into the subsurface: in the context of the production of shale gas through hydraulic fracturing or in form of wastewater disposal from conventional and unconventional oil and gas production. Geological gas storage or geothermal systems also involve the injection of fluids. Injections affect the field of fluid pressures H. Class holger.class@iws.uni-stuttgart.de 1 Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany 2 Swiss Seismological Service, Swiss Federal Institute of Technology, ETH Zürich, Zürich, Switzerland as well as the stresses and deformations of aquifers and caprocks. Caprock integrity might be compromised. Thus, the need to understand the coupling between hydraulics and geomechanics is evident. Numerical simulation of coupled hydraulic and geomechanical processes requires the choice of an appropriate coupling scheme. This can be discussed by balancing numerical accuracy and computational efficiency and might depend on the temporal and spatial scales of the application, on the complexity of the relevant physics, on the heterogeneity, or on other scenario-specific aspects.
The classical approaches for flow and geomechanics fall into two categories: fully coupled and sequential methods. While different naming for these schemes is in use, as stated in White et al. [47], the differences are clear: The unknowns of flow and geomechanics are either solved simultaneously for one time step (termed fully coupled or monolithic) or sequentially by splitting up the coupled problem into sub-problems. One of the simplest sequential schemes is the one-way coupling. It starts with the fluid-flow subproblem and determines the effect on the geomechanics afterwards in a post-processing step. The flow problem might also account for the increase in pore volume as a result of increased pore pressures using an approach for pore compressibility. For such a scheme, the pressure evolution and the storage capacity for CO 2 injections have been studied using the TOUGH2 simulator (e.g., Zhou et al. [50], Birkholzer et al. [6]), as well as the geomechanical impact of injections in the context of hydraulic fracturing [41] and CO 2 storage [39] with the TOUGH-FLAC code. In this conceptional idea, the geomechanics part can give feedback to the flow calculation via an adapted permeability from the solution of the geomechanics sub-problem in the previous time step.
Similar, but more accurate, is a scheme where flow and geomechanics are iterated within the same time step. Preisig and Prévost [37] compare the results of a one-way coupling method, an iterative approach, and a fully coupled approach with the analytical solution of Mandel's problem. They arrive at the conclusion that a large number of iterations is needed to get accurate results for the sequential scheme. Thus, they only include a one-way coupling and a fully coupled approach for their investigation of the CO 2 injection at In Salah, Algeria. Sequential schemes were significantly improved by Kim et al. [27,28] and Mikelić and Wheeler [34]. They revealed that the undrained split and the fixed-stress split are preferable over the drained split and the fixed-strain split due to their convergence behavior and stability. This explains the large number of iterations reported by Preisig and Prévost [37]. Yoon and Kim [49] contributed to this analysis with a comparison of the fixedstress and the fully coupled approach with respect to the stability of the different discretization schemes.
In the context of sequential schemes, it is a special case of iterative coupling when the calculation of flow and mechanics is not repeated within the same time step but instead information of the mechanics is passed back into the fluid-flow calculation only for the subsequent time step. In contrast to the pure one-way coupling, such a scheme thus includes a feedback of the mechanics part onto the fluid-flow equation. We denote this special case here as an explicit coupling. The aforementioned coupling via permeability updates within the TOUGH-FLAC code is in that sense an explicit coupling in the same way as the recent expansion to include a correction term based on the fixed-stress approximation [7,29].
There is a variety of possibilities in developing the details of a sequential scheme; therefore, we put the focus of this study on the end members of the sequential spectrum: the one-way coupling and the fully iterated fixed-stress scheme.
This study evaluates suitability of solution schemes for flow and geomechanics by analyzing how the numerical accuracy in describing the relevant physical processes can be affected by the choice of scheme. The implementations are done in the numerical toolbox DuMu x [1,18,19]. We show briefly that DuMu x can correctly reproduce (a) the analytical solution for Mandel's problem and (b) the results of a one-way coupled version of TOUGH-FLAC for a homogeneous single-phase injection scenario. Then, only with the respective implementations in DuMu x , a comparison is made between the one-way coupling scheme, an iterative sequential scheme, and a fully coupled scheme for the same homogeneous single-phase example as before and for a CO 2 injection test case. A focus will be on singlephase versus multiphase physics, which has not been done so far in a systematic manner, and on the importance of transient effects where we will show that the schemes most notably differ. After that, there are some considerations on numerical efficiency provided before conclusions are drawn and an outlook is given.

Model and numerical solution schemes
The numerical simulator DuMu x [1,18,19] is an opensource package designed to model flow and transport processes in porous media. It is based on the Distributed and Unified Numerics Environment (DUNE) [3,4,8,9]. The DUNE-ALUGrid module [2] was used as the grid manager in this study. Except for the fully coupled model, which was already implemented within the doctoral thesis of Darcis [17], all models on top of DuMu x and DUNE were developed within the scope of this work.

Governing equations
The mass balance equation for multiple fluid phases in porous media and the momentum balance for the solid are introduced. The concept of effective porosity is key to describing the coupling of flow and mechanics.
Mass balance equation of the fluid: For a fluid phase α, the mass balance equation is: φ is the porosity, α the density, and S α the saturation of phase α. v α is the Darcy velocity of a phase α: and the sources and sinks of a phase α are denoted by q α . The indices w and n are for wetting and non-wetting phases, respectively. Inserting Darcy's Law in the mass balance equation gives: Effective porosity Let us briefly recapitulate the concept of effective porosity. In general, constitutive relations of poroelastic materials can be formulated dependent on dynamic variables, like stress or pressure, and on kinematic quantities, like strain or velocities. The recent paper of Steeb and Renner [42] provides a review of mechanics of poroelastic media. They highlight the role of porosity in linear poroelasticity and demonstrate how alternative pairs of kinematic variables can be used for formulating constitutive relations. In our case here, we follow the approach of Han and Dusseault [22] in our brief outline below, but we note that alternative formulations are found in the literature. The constituents of the porous medium show different responses to stress and pressure changes. The resulting change in porosity can be defined as: V p is the pore volume and the bulk volume is V b . These volumes deform differently; hence, different bulk moduli for the drained and undrained porous medium as well as for the solid grains, the fluid, and the pore volume are defined. The porosity change is then a function of the pore pressure p and the volumetric (or mean) stress σ v (positive for compression): This employs the drained bulk modulus K dr and the bulk modulus K s of the solid grains. The bulk modulus of the fluid is implicitly considered in our equations where the fluid density is formulated dependent on pressure. Han and Dusseault [22] further replace the stress with volumetric strain, so Eq. 5 becomes: The relation in Eq. 6 is dependent on the dynamic variable pressure and the kinematic variable strain. For rigid grains (K s → ∞, equivalent to a Biot's α = 1), this reduces to: Now, the effective porosity φ eff can be obtained from the volumetric strain v and the initial porosity φ 0 (assuming v,0 = 0): or as a function of the displacements by referring to the volumetric strain v as div u: In case the denominator in Eqs. 8a and 8b is considered negligible, it is simply: An alternative to deriving the porosity change from the interacting components of the porous material and their respective bulk moduli is the definition of a pore compressibility: dependent only on the changes in pore pressures and ignoring the stress contribution. The pore compressibility is related to K p , the bulk modulus of the pores, and K dr by: In a linearized form, the effective porosity is then obtained from the difference between the initial pressure p 0 and the current pressure p: As indicated further above, there are many more possibilities to formulate the incremental effects of changing fluid pressures and stresses on porosity. Later on, we use e.g. also, in order to match the analytical solution of Mandel's problem Phillips and Wheeler [35], an approach for the porosity that is expressed in terms of the Skempton coefficient, which in turn is related to the bulk modulus of the fluid, the porosity, and the undrained modulus. We shall briefly come back to this discussion then.
The change in porosity affects also the permeability tensor, so the intrinsic permeability tensor K becomes an effective permeability tensor K eff . A classical and commonly employed approach is the Kozeny-Carman relationship [14,30].
The balance equation is then written as: If the effective porosity is considered to be only a function of pressure, it can be calculated from Eq. 12, while a dependence of the porosity on the volumetric strain (8b) introduces the displacements as additional unknowns. These are provided by the momentum balance equation for the rock matrix.

Momentum balance equation for the rock matrix
Assuming quasi-static conditions, the conservation of momentum is div σ + g = 0 . (14) g is the gravity vector here. The bulk density b of the porous medium can be calculated from the weighted densities of the fluid phases w , n and the rock matrix rm : The effective stress tensor σ accounts for the contribution of the pore pressure p: The pore pressure is calculated as the effective pore pressure p eff : With this, Eq. 14 becomes: This momentum balance equation can be linearized by subtracting the initial state (subscript 0) for effective stress, effective pressure, and the bulk density: This yields: Under the assumption of small porosity changes ( φ ≈ 0, (1 − φ) ≈ 0 ) and a constant density of the rock matrix ( rm ≈ 0), the bulk density change becomes: and further simplifies to assuming a compressible non-wetting phase and an incompressible wetting phase in the linearization of b . With these simplifications, the momentum balance of the rock matrix is finally:

Constitutive equations and supplementary constraints
The system of equations is completed by employing the following constitutive equations and constraints: The strain is derived from the displacement vector u as: Linear elastic behavior of the rock obeying Hooke's Law is assumed. Young's modulus E and Poisson's ratio ν are used as the two required elastic constants.
The sum of the fluid saturations adds up to 1: The pressures for the different phases are connected via the capillary pressure, calculated here after van Genuchten [44]): with with the parameters α VG , m VG , and n VG characterizing porous media properties and the residual wetting-phase saturation S rw . The Van-Genuchten formulation is also used for the relative permeability:

Discretization and solution schemes
Using an equal-order approximation of the pressure and the solid displacement in the monolithic setting is known to create stability problems and is addressed by a staggered-grid approach by [26]. Darcis [17] avoided the encountered spurious pressure oscillations by discretizing the balance equation of the fluid phases with the Box method (cf. [23]) and the momentum balance equation with the Standard-Galerkin finite-element method. By applying different weighting functions, this leads to a quasi-staggered approach that allows for a nodal-based approximation of all primary variables. In this way, we achieve stability for the monolithic setting although we cannot prove this in a rigorous manner. Many well-established software codes such as Eclipse by Schlumberger or TOUGH2 by the Lawrence Berkeley National Laboratory [38] use the cellcentered finite-volume method, e.g., [15]. The fixed-stress approach for this study is implemented in DuMu x for both the Box method and a cell-centered finite-volume method for the two-phase flow equations. A fully implicit Euler scheme is used to for time discretization. The discretized mass balance equation for the fluid phases as in Eq. 13 contains the vector of pressure values of the wetting phase and the saturation values of the non-wetting phase as primary unknowns. For convenience, the suffix "2p" will be used when referring to the twophase equations while "el" refers to the momentum balance equation of the solid in Eq. 25.
Combining the two balance equations yields their lefthand sides as a residual vector r(w n−1 , w n ) with the components r 2p and r el : Here, the residual vector r(w n−1 , w n ) is a function of the solution vectors w n−1 and w n , which contain the discrete values of two subsequent time steps t n−1 and t n for all nodes. Both solution vectors (exemplarily written here only for w n ) include a 2p part w n 2p with p n w and S n n and an elastic part w n el with u n : In the following, we omit the time-step index n for simplicity.
As this system of residual equations is strongly nonlinear, a Newton scheme is employed. For the k-th Newton iteration, w k is the k-th estimate of the solution. To determine the new solution for k + 1, the Jacobian matrix J k is used. Here, J k is the matrix of all first-order derivatives of the residual vector at the k-th iteration: The new solution vector w k+1 is then determined by solving and by adding the update w to the previous solution vector w k Equation 32 describes the block-wise structure of the residual vector r. In the same sense, the Jacobian system solved in each Newton update can be expressed as: where J a,b denotes balance equation a differentiated with respect to solution vector b. As an example, J 2p,2p is the derivative of the fluid's mass balance equation with respect to p w and S n . The components J 2p,el and J el,2p account for the flow-geomechanics coupling.
Fully coupled scheme In the fully coupled scheme, the primary unknowns p w and S n for the two-phase flow and u for the geomechanics are solved simultaneously for each time step. This means that every Newton update is computed from the linear system in Eq. 38. The use of the full matrix J explains why this scheme is also often referred to as monolithic.
Fixed-stress scheme A scheme where the unknowns for flow and geomechanics are solved sequentially requires to split up the full problem into a flow and a geomechanics sub-problem. We denote here a coupling step as including first the flow problem, then the geomechanics sub-problem. The simplest version of an iterative coupling would be a transfer of p w and S n to the geomechanics, then inserting calculated stress and strain back into the flow problem for the next iteration. Such a drained-split scheme is only conditionally stable. A scheme with unconditional stability is the fixed-stress split [28,34]. The scheme requires that the difference between the volumetric stress of the flow problem σ n,i v, 2p and the volumetric stress σ n, i−1 v, el of the previous geomechanical solution is 0: Here, i is the index of the coupling step. Since a coupling starts with the flow problem, the last geomechanical solution was obtained in the previous coupling step and thus is denoted with the index i − 1. For the geomechanical problem, the pressure of the previous flow problem is prescribed, so p n, i el and p n, i 2p are equal within a coupling step i: Thus, one can express Eq. 39 by the pressures p n, i 2p and p n, i−1 2p and the volumetric strains n, i v, 2p and n, i−1 v, el , Here, K dr is the drained bulk modulus and the Biot coefficient α is assumed to be 1. Let us briefly note that for consolidated rocks where further pore compression is inhibited, Biot's α maybe less than 1, e.g., shown in Ingraham et al. [24] as dependent on reservoir stress conditions. In such cases, a more consistent account would require the Biot coefficient to be considered throughout the above derivation. The coupling within the flow part arises from the dependence of the porosity on the volumetric strain (8b). With Eq. 41, the flow problem can be formulated independent of the current displacement vector u by calculating the volumetric strain n, i v,2p of the flow problem in the current coupling step i from the pressure difference of the current and the previous coupling step of the flow problem and from the volumetric strain n, i−1 v,el determined in the previous coupling step for the geomechanics: For the two-phase system, the mass and momentum balances then depend on the primary variables p n, i w and S n, i n , as the pressure p becomes the effective pressure p eff : p eff = S w p w +S n p n , ( 17 revisited) neglecting the interfacial energy term in the equivalent pressure concept by Coussy et al. [16]. After solving the flow problem altered by the fixed-stress assumption, p n, i w and S n, i n are inserted into the momentum balance of the geomechanics, where only the displacement vector u n, i remains as the primary variable. With this, two separate non-linear problems are formulated and yield the following linear systems to be solved in each step of the corresponding Newton scheme: Here,J 2p,2p andJ el,el denote the derivatives of the modified balance equations:J 2p,2p contains the fixed-stress assumption, so the volumetric strain of the current coupling step i and the current Newton iteration k for a time step n is calculated from: and the effective pressure values from the flow problem are prescribed withinJ el,el . When the porosity update is calculated using Eq. 9, it is worth noting that a fixed-stress split with just one single coupling step is identical to a calculation of the porosity change from the pore compressibility (cf. Eq. 12). The previous coupling step becomes the previous time step for zero iterations. Thus, i−1 v is equal to 0 for the first time step and dependent only on the pressure difference between the previous and the current time steps for all following time steps. As these collapse to the difference between the initial and the current pressures, this simplified sequential scheme becomes equivalent to using the pore compressibility or the drained bulk modulus, respectively: Figure 1 further illustrates the differences between the discussed numerical coupling schemes. The top panel shows a numerical routine as it can be found, for example, in the original version of TOUGH-FLAC (a sequential combination of the TOUGH2 multiphase fluid and heat transport code with the commercial FLAC3D geomechanical simulator as described in [40]), where the pore compressibility is used to calculate the porosity changes within the flow problem. At time t n , the porosity is a function of the current effective pore pressure p n eff . Next, p n w and S n n are transferred to the momentum balance equation to calculate the solution at time t n . A feedback of the mechanics to the flow can be achieved via the transfer of a changed permeability k, which is used in the next flow time step. If the permeability is constant, the mechanics are just a post-processing step and this one-way coupling scheme will be referred to as a zero-iteration case here.
A more recent version of TOUGH-FLAC [7] transfers a correction term for the porosity calculated by using the fixed-stress assumption back into the flow problem (see Fig. 1, middle). Again, it is important to note that this correction becomes active only in the next time step of the flow problem. Thus, the porosity within each flow time step is a function of the current effective pore pressure p n eff and the volumetric strain of the previous time step n−1 v, el . As explained earlier, this is denoted here as an explicit coupling approach.
An iterative fixed-stress scheme repeats the coupling steps several times until convergence (see Fig. 1, bottom). Within the flow problem, the volumetric strain n v, 2p is calculated from Eq. 42. This implies (i) that the porosity depends on the values of p n eff and n v, 2p of the current time step, given that at least one iteration is performed; and (ii) that the approximated n v, 2p becomes n v, el of the geomechanics when the iteration is continued until convergence.
In summary, the use of the pore compressibility and the subsequent solution of the mechanics represents the simplest scheme. The coupling can be strengthened by a feedback of the mechanics back to the flow for the next time step via an update of permeability or porosity. A sophisticated iterative scheme requires to repeat the coupling step several times within each time step. In this study, we examine the end members of the spectrum, i.e., the first and the last schemes, and show for selected scenarios the circumstances under which significant deviation between the schemes occurs.

Verification of DuMu x implementations
Before we elaborate on the occurrence of differences for the investigated schemes, we choose to verify our implementations in DuMu x against the analytical solution Fig. 1 Numerical routine for a one-way coupling, such as the original version of TOUGH-FLAC using the pore compressibility, (top) for a more recent version using the fixed-stress assumption for an explicit coupling (middle) and for the iterative fixed-stress scheme as implemented in DuMu x (bottom) of Mandel's problem and a simplified one-way coupling mode of the well-established codes TOUGH-FLAC for a homogeneous single-phase scenario. Later on, we analyze then the schemes by using only the respective implementations in DuMu x .

Mandel's problem: comparison with the analytical solution
Mandel's problem was originally described in Mandel [32] and is a widely used setup to verify an implemented model  [10,29,35]). It consists of a poroelastic specimen with rectangular cross section of 2a in x-direction, 2b in y-direction, and infinite length in the z-direction. It is sandwiched in between two rigid plates and compressed by forces normal to those plates. For symmetry reasons, the computational domain can be limited to one-quarter of the physical domain (see Fig. 2).
Regarding our implementation of the effective porosity, we use here a formulation from Phillips and Wheeler [35] who summarize the analytical solution of Mandel's problem and use a fixed value of B = 0.833 for Skempton's coefficient. We have elaborated on the effective porosity before in Section 2.1. We use here the incremental fluid content ζ which relates to the effective porosity via the difference in volumetric strains of fluids and porous matrix. For an initial reference pressure p = 0, as we apply it in Mandel's problem, the time derivative of ζ is balanced with the divergence of the Darcy velocities and is thus used instead of porosity in the mass balance equation. For the necessary algebra and the details of transforming the formulations of poroelasticity dependent on different dynamic or kinematic variables and the various coefficients of stress-strain relations, we refer again to Steeb and Renner [42]. This formulation to model Mandel's problem then expresses incremental fluid content ζ as a function of Skempton's parameter B and the undrained bulk modulus: 1−B B K b can also be interpreted as the specific storage capacity, which is the inverse of the storage modulus M used by Biot [5].
The pressure evolution given by the analytical solution and the numerical results obtained with the fully coupled scheme implemented in DuMu x are shown in Fig. 3 for different points in time. The numerical pressure evolution is in excellent agreement with the analytical solution. This verifies the implementation and allows to proceed with a more application-based scenario.

Homogeneous test case: comparison with TOUGH-FLAC
Scenario description Water is injected into a homogeneous, two-dimensional domain fully saturated with water. The domain is 2 km × 2 km in size and located at a depth of 500 to 2500 m. The pressure distribution is hydrostatic (9.81 MPa/km) and atmospheric pressure of 0.1 MPa is assigned to the ground surface. This results in a pressure of 5 MPa at the top and 24.6 MPa at the bottom. More details can be found in Fig. 4. The grid for this domain is the same as in Fig. 8b and is later on used to include a fault in the scenario (Section 4.2); thus, this feature and the corresponding refinement is included already here. In total, 0.2 kg/s/m is injected into the two cells with the center coordinates (6.25, 2.5) and (6.25, −2.5). The rate is high enough to produce a noticeable pressure response throughout the system.
The time steps were chosen to be the same for all schemes, also for the comparison study discussed later on, to achieve a better comparability. They are small at the beginning (starting with 8.64 s) and increase over the course of the simulation.

DuMu x versus TOUGH-FLAC in one-way coupling mode
For the verification, we compare here simplified coupling strategies in DuMu x and TOUGH-FLAC [40]. Both are Fig. 4 Setup of the homogeneous single-phase test case run here in one-way coupling mode, which means that permeability remains constant and porosity changes are simply dependent on the pore compressibility. For DuMu x , we use also the denotation "fixed-stress, zero iteration" to distinguish from the iterated fixed-stress scheme.
The pressure evolution at the injection is shown in Fig. 5. Panel (a) displays the entire simulation time of 1800 days. The pressure increases rapidly, then slows down and reaches a plateau value when a state of equilibrium for the injection and the outflow through the model boundaries is achieved. There is hardly any difference visible between the different datasets. Both one-way coupled schemes (TOUGH2 using the pore compressibility and the fixed-stress implementation with zero iterations in DuMu x ) are in excellent agreement.
The results of the zero-iteration fixed-stress scheme in DuMu x and TOUGH-FLAC match perfectly also for the displacement and the volumetric strain as Fig. 6 shows.
In summary, the results of a zero-iteration fixedstress scheme implemented in DuMu x can be verified with those obtained with TOUGH-FLAC using the pore compressibility.

Homogeneous test case
We still consider the homogeneous test case as in the verification example above.  5 Full pressure evolution (a) and details of the increase (b) at the injection simulated using TOUGH-FLAC, and the cell-centered (CC) discretization in DuMu x for a fixed-stress scheme with just one coupling step (equal to zero iterations) and a fully converged simulation A detailed look during the first 80 days in Fig. 5b proves that an iterated fixed-stress scheme shows a small deviation during the transient period of the pressure increase. Both one-way coupled models underestimate the pressure increase during this period compared with the iterated fixed-stress scheme.
The underestimated pressure of the zero-iteration fixedstress scheme and the pore compressibility approach, in comparison with the iterated scheme, also affects the results of the geomechanics as seen in the displacement in xdirection and the volumetric strain in Fig. 6. Figure 7a shows the evolution of the pressure over the iterations of the fixed-stress scheme for the first time step. The zeroth iteration underestimates the pressure in comparison with the fully coupled value. In other words, this corresponds to an overestimated porosity change due to the pressure change. As soon as this is corrected in the first iteration by the volumetric strain value obtained from the previous solution of the geomechanics, the pressure estimate improves significantly. After the fourth iteration, the values are almost identical. After a sufficient number of iterations, the results of the iterated fixed-stress scheme are identical to those of the fully coupled model. This confirms that the underestimation observed for the one-way couplings, i.e., pore compressibility or zero fixed-stress iterations, is a result of the simplified representation of the flowgeomechanics coupling. A variation of this plot in Fig. 7 (right) displays exemplarily the results of the drainedsplit scheme where the flow problem is solved using the strain and stress fields of the previous geomechanics iteration, while the pressure is kept constant during the geomechanical step. It is obvious why Kim et al. [28] refer to this scheme as being, at best, conditionally stable.
In summary, for the transient phase of the simulation, a small but noticeable underestimation of the one-way coupling schemes was observed relative to a fixed-stress scheme with iterations.

CO 2 injection heterogeneous test case
Scenario description This scenario involves the injection of CO 2 into a geological formation. This adds complexity in two aspects compared with the previous scenario: first, it is a multiphase flow system, and second, it includes heterogeneities in the form of different layers and a fault zone (see Fig. 8). Similar setups have been used in previous studies, e.g., [33] and [39].
The point of injection is located in the center of a sandstone reservoir with a permeability of 1 · 10 −13 m 2 and a thickness of 100 m. Above and underneath the target reservoir is a 150-m-thick shale layer serving as a caprock with a permeability of 1 · 10 −19 m 2 . A 1-km-long fault zone dipping 80 • cuts through the reservoir and the shale layers. Its midpoint is at a depth of 1500 m and at a horizontal distance of 500 m from the injection point. With a permeability of 1 · 10 −15 m 2 , it is less permeable than the target reservoir but more permeable than the shales. Table 1 lists the parameters assigned to the different layers. The temperature is increasing with depth from 22.5 to 72.5 • C. The values for the viscosity and the density are temperature dependent, but the simulation itself solves no thermal energy balance equation.
For the capillary pressure-saturation relationship and the relative permeability, the formulations by van Genuchten [44] are used (see Eqs. 28, 30, and 31).
The injection rate is here at 0.02 kg/s/m. The grid is the same as before (see Fig. 8) and the Box method is used for the discretization of the two-phase flow equations.
Time steps are again the same for all scenarios and increase over the course of the simulation. Figure 9 illustrates how the CO 2 saturation S n and the water pressure p w evolve in the domain. Their  distribution is shown after 200 days and 1800 days. To highlight the location of the fault and the shale layers, a gray filter with different opacities was added. The less permeable a cell is, the more opaque it is. Driven by buoyancy, a large portion of the CO 2 accumulates right underneath the shale caprock. The maximum values of CO 2 saturation at the end of the simulation after 1800 days are S n = 0.56 above the point of injection and S n = 0.54 within the fault zone. Upon CO 2 reaching the fault, it moves also upward there as can be seen in Fig. 9, bottom left. The right column of Fig. 9 shows the corresponding pressure distribution in the domain. The injection pressurizes the reservoir section bounded by shale layers and the fault and also the fault itself. The wetting pressure p w reaches a peak value between 100 and 200 days and then declines.

Results
The results as presented up to now were obtained with the fully coupled approach. We have shown before in the homogeneous setup that the results of the fully coupled scheme and the iterated fixed-stress scheme are in total agreement. The same holds true for this case here as Fig. 10 shows.
The plot shows again that the pressure at the point of injection is underestimated when no iterations are performed. Dependent on the desired accuracy, it might be sufficient to use three or four iteration steps.
In the homogeneous single-phase scenario, we observed the most distinguishable differences between the zeroiteration scheme and the fully coupled scheme during the transient evolution of the system. This behavior is even more pronounced in this scenario here with layers of different hydraulic and mechanical properties, with a fault-zone feature and a second fluid phase as can be seen in Fig. 11a.
The CO 2 arrives roughly 50 days later at the fault in the zero-iteration case (see Fig. 11, right). In addition to the temporal evolution in Fig. 11, the spatial distribution of Fig. 9 Spatial distribution of the CO 2 saturation S n (left column) and the wetting pressure p w (right column) at t = 200 days (top) and t = 1800 days (bottom) in time during the pressure decrease. Again, the gray colors in the saturation distribution chart correspond with the permeability (the more gray, the lower the permeability) the differences in pressure and saturation between the fully coupled and the zero-iteration fixed-stress schemes is shown in Fig. 12. After 100 days, the relative differences in the wetting-phase pressure p w are the highest in the vicinity of the reservoir/shale layer contact (see Fig. 12, upper).
The absolute differences in the non-wetting phase saturation S n after 200 days displayed in Fig. 12 (lower) are concentrated at the CO 2 front.

Scenario variations
The pressure peak, and the following decline, is remarkably different from the behavior observed Fig. 10 Pressure at the injection over sequential iterations for fully coupled and fixed-stress schemes in the homogeneous case. We discuss now some variations of the scenario aiming (i) at elaborating a better understanding of the different effects introduced by adding the layers and the fault as well as a second fluid phase and (ii) at getting a feeling for the order of uncertainty that is introduced by the numerical scheme, fully coupled/ iterated fixed-stress versus non-iterated scheme, compared with the changes due to added physical complexities (geologic uncertainty Walter et al. [45]). Figure 13a shows pressure over time for a single-phase water injection compared with the pressure during the CO 2 injection with everything else in the setup identical. As in the homogeneous case, single-phase injection leads to a pressure plateau also in this heterogeneous setup. The fault has an effect comparable with the constantpressure boundary in the homogeneous setup. It leads to a quasi-constant pressure in the single-phase case as soon as an equilibrium between the injection and the flow leaking through the fault is established. The two-phase case is distinctly different with an increase and subsequent decline caused by relative permeability. Initially, with no CO 2 in the reservoir, the CO 2 relative permeability is close to zero and a high resistance to CO 2 flow is given. Accordingly, the injection pressure rises sharply. Subsequently, CO 2 saturation and relative permeability increase; thus, resistance becomes less and pressure required to uphold the injection flux declines.
Two further scenario variations aimed at investigating the role of the fault: (i) an increased fault permeability of  Fig. 13b. The observed pressures are even lower for the high permeable fault, which lets us conclude that the fault's breaching of the caprock outweighs the barrier effect of the fault itself since it makes a larger volume accessible for the fluids to give away for the injection.
In summary, comparing the pressure deviations in the order of 0.1-1.0 MPa between the cases of zero iteration and fully coupled with those deviations we observed above in varying the geological features, we find them in the same order of magnitude. Thus, the point we want to make here is that not iterating between flow and geomechanics may lead to comparable compromises in the accuracy of results as does a disregard of a geological feature like a fault.

Efficiency considerations
Let us extend the comparison of the different schemes by adding some comments on the efficiency of the fully coupled and the sequential schemes. We will only treat the case of employing direct solvers as an academic exercise to represent a reference. Concerning iterative solvers, the situation becomes more involved and is therefore beyond the scope of this study. We refer here to related studies on iterative solvers and preconditioners, e.g., in Gaspar and Rodrigo [20], White et al. [47,48]. The key difference between the approaches with respect to efficiency is the size of the linear systems to be solved. According to Golub and van Loan [21], it requires 2Npq floating-point operations to solve a linear system of N unknowns with a direct solver. p and q are upper and lower limits of the bandwidth. In our two-dimensional case, the ideal bandwidth would be √ N. The required operations O for a linear system of N unknowns should then be proportional to 2N 2 in an asymptotic consideration: with C as the unknown proportionality constant including also the factor 2.
Ideally, it should be expected that for a direct solver in a two-dimensional case, the computational cost of the fully coupled scheme with twice the number of unknowns is fourfold compared with just solving for p w and S n using the fixed-stress scheme. The same is true for solving only the geomechanics for the two components of the displacement. But since a coupling step requires solving both flow and geomechanics, its cost is the sum of solving both parts, being then computation-wise half as expensive as the fully coupled scheme. Theoretically, this means that two coupling steps, i.e., one iteration, can be performed until this advantage is used up. Taking into account the results obtained for the convergence of the fixed-stress scheme over the iterations, as in Figs. 7 and 10, one iteration might not be sufficient in terms of accuracy.
In three dimensions, the ideal bandwidth is N 2 3 . The difference in unknowns, i.e., 5 for the fully coupled scheme and two plus three for the flow and geomechanics subproblems of the fixed-stress scheme, results in a speedup factor of 2.37 for the fixed-stress scheme (5 7/3 / (2 7/3 + 3 7/3 )). Thus, if more than two iterations are performed, this gain is again used up.
The sequential scheme's main advantage is the possibility of combining different simulation codes. Nonetheless, its structure of nested iterations, i.e., iterations to solve each sub-problem, e.g., via Newton's scheme, as well as iterations between the sub-problems, leaves room for creativity and possible efficiency gains. For instance, Borregales et al. [10] propose a partially parallel-in-time fixed-stress splitting method for better parallelization. Solving flow and mechanics using an inexact Newton scheme for a certain number of iterations and only solving them exactly for later iterations between flow an mechanics could be an option, too.

Conclusions and outlook
This comparative study has put the focus on the end members of the spectrum of approaches to solve flow and geomechanics represented by the one-way coupling using pore compressibility (or zero-iteration fixed stress) on the one hand, and the iterated fixed-stress scheme which has been shown to achieve identical accuracy as the fully coupled scheme on the other hand.
A comparison with the analytical solution of Mandel's problem has verified this new implementation of a model for hydro-and geomechanical processes in the simulator DuMu x . The simulation of a homogeneous numerical test case with a one-way coupled version of DuMu x and the established TOUGH-FLAC code serves the same purpose.
As should be expected, for the homogeneous scenario, this implementation of the iterated fixed-stress scheme and the fully coupled scheme can produce identical results after a certain number of iterations.
The heterogeneous setup introduced different layers, a fault zone, and a second fluid phase, thus increasing the physical complexity and allowing quantitatively comparing errors introduced by not-considered geologic features and physics with errors due to inaccurate coupling. The oneway coupling with the pore compressibility and without iterations can lead to considerable deviations from the fully coupled solution in particular since transient effects are not captured accurately. The effect of relative permeability in this setup with two fluid phases has a particular pressure curve over time featuring a peak value of pressure and subsequent decline, which is notably distinct from a singlephase system. This study demonstrates the importance of accurate coupling schemes to capture this transient period properly.
Of course, often uncertainty in the geologic parameters, like permeability, porosity, elastic constants, etc., may dominate over the errors introduced by an "inaccurate" coupling scheme. However, the heterogeneous scenario with CO 2 injection scenario revealed pressure differences between the schemes were at times in the same order as the geologic uncertainty. A difference in pressure prediction of 1 MPa might be considered small relative to the total pressure of several tens of megapascal. However, if it concerns, for example, the assessment of the potential for induced earthquakes, 1 MPa can be decisive for the failure criterion to be surpassed or not.
Theoretical considerations on computational efficiency suggest that a potential advantage of the sequential scheme is used up by the iterations required to reach better accuracy. However, the gain in accuracy with only one iteration is already significant; it seems a good trade-off between computational costs and numerical accuracy to run the fixed-stress scheme with one iteration.
Explicit coupling schemes might also be an option. They update porosity and permeability to be effective only in the next time step by using the result of the geomechanics from the previous one. Thus, they are situated in-between the zero-iteration and one-iteration fixed-stress schemes.
As an outlook: there are further options that might be studied. Implementing an algebraic multigrid solver [43] or applying the fixed-stress as a pre-conditioner as suggested by White and Borja [46] could increase the efficiency of the fully coupled model even more. Both et al. [12] applied the fixed-stress scheme to a sequential scheme coupling linear elasticity and flow in unsaturated porous media modeled by the Richards equation. As they encountered difficulties when using Newton's method, they propose using the L-scheme for the Richards equation [31,36] which can be described as a standard Picard iteration with additional diagonal stabilization, resulting in an efficient and robust decoupling of the equations for geomechanics and flow. The same group also developed an optimized tuning parameter, which depends on all mechanical parameters and replaces the classically used drained bulk modulus in the fixed-stress scheme [11]. These are promising new developments that could be incorporated into the existing approaches. This study could show that it is important to put effort into developing accurate and efficient numerical schemes for the coupled solution of flow and geomechanics.
Code availability The numerical simulator DuMu X used in this study can be obtained at http://www.dumux.org/ download.php. The specific code used is available at https:// git.iws.uni-stuttgart.de/dumux-pub/beck2018a Acknowledgments Open Access funding provided by Projekt DEAL. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article 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/.