1 Introduction
In a variety of applications, it is useful to model the hydromechanical behavior of porous media infiltrated by one or more fluids—e.g. in geotechnical engineering Teatini et al. (2006); Borja and White (2010); Borja et al. (2012); Camargo et al. (2016); Cao et al. (2018); Fávero Neto and Borja (2018); Ghaffaripour et al. ; Gholami Korzani et al. (2018); Mikaeili and Schrefler (2018); Navas et al. (2018); Oka et al. (2019); Prassetyo and Gutierrez (2018); Song et al. (2018); Wang et al. (2018, 2019); Yan et al. (2019); Zhou et al. (2019), hydrocarbon recovery Zoback (2007); Feng and Gray (2019); Mashhadian et al. (2018); Peng et al. (2019); Semnani et al. (2016); Shiozawa et al. (2019); Yan and Jiao (2019); Zhang et al. (2019); Zhao et al. (2018); Zhao and Borja (2019), and geologic carbon storage Verdon et al. (2013); White et al. (2014); Jin and Arson (2019). Precise models should account for the tight interaction between solid deformation and fluid flow. The conceptual framework for modeling this coupled behavior is well established de Boer (2012); Coussy (2005), with Biot’s work A. (1941) providing a sound theoretical foundation. Computational methods for poromechanics, however, still pose many interesting challenges. In particular, this work focuses on numerical instabilities that may arise due to the discretization spaces chosen for the coupled fields.
As a representative problem, we consider a model in which an elastic solid skeleton is saturated with two immiscible fluids. We present a multiphase formulation for its relevance in many geoscience applications, though the singlephase formulation is a straightforward subcase. The behavior of the porous system is governed by a momentum balance equation for the mixture and mass balance equations for each of the fluids. A fullyimplicit time integration strategy is adopted, where all unknown fields are updated simultaneously in a monolithic manner White et al. (2016, 2018). A variety of finiteelement and finitevolume based discretization strategies may be applied to these equations, each with certain advantages Settari and Mourits (1998); Kim et al. (2011); Prévost (2014); Castelletto et al. (2015); Nordbotten (2014, 2016); Honório et al. (2018); Choo and Lee (2018); Sokolova et al. (2019)
. Of specific interest here, however, is a frequent choice: continuous trilinear interpolation for the displacement unknowns and elementwise constant fields for the pressure and saturation unknowns. Such an interpolation results, for example, when applying a continuous Galerkin finite element scheme to the momentum balance equation, and a finite volume scheme to the mass balance equations
Settari and Mourits (1998); Kim et al. (2011); Prévost (2014); White et al. (2018).This discretization works well in a variety of practical cases. The chosen interpolation spaces, however, can be problematic. In the limit of small time steps or low permeabilities, undrained deformation can occur. The fluid mass balance equations impose an incompressibility constraint on the deformation field. Like many other constrained problems—e.g. Stokes flow or incompressible elasticity—these divergence constraints can create numerical instabilities if the discrete approximations for the field variables do not satisfy the LadyzhenskayaBabuškaBrezzi (LBB or infsup) condition Babuška (1973); Brezzi (1974). Unfortunately, the combination explored here is not intrinsically LBBstable. As a result spurious oscillations, i.e. checkerboarding, may be observed in the pressure field. A less obvious, but equally important, symptom is a degradation in the convergence rate of iterative linear solvers. Nearzero eigenvalues associated with spurious modes will control the conditioning of the system matrices, leading to poorly conditioned systems and increased iteration counts. This latter issue can persist in regimes where checkerboarding is not visually apparent, and may thus go unnoticed by practitioners.
One way of circumventing these instabilities is to introduce a carefullydesigned perturbation to the constraint equations. The goal is to remove instabilities while maintaining an accurate approximation of the underlying PDEs. This is the basic rationale behind many stabilization techniques, including the Brezzi and Pitkäranta scheme Brezzi and Pitkäranta (1984), the Galerkin LeastSquares approach Hughes et al. (1986), and the Polynomial Pressure Projection technique Dohrmann and Bochev (2004). Various stabilization schemes have been proposed that devote particular attention to constant pressure elements Fortin and Boivin (1990); Bochev et al. (2006); Burman and Hansbo (2007); Hughes and Franca (1987); Silvester and Kechkar (1990), starting with early work in Pitkäranta and Saarinen (1985). In Hughes and Franca (1987) the idea of penalizing the pressure jump across interelements boundaries was introduced. An important modification of this method, the Local Pressure Jump (LPJ) stabilization, was developed in Silvester and Kechkar (1990) based on the macroelement concept.
The schemes above were primarily developed for fluid mechanics problems. Since then, many of these stabilization schemes have been successfully applied to poromechanics with singlephase flow White and Borja (2008); Preisig and Prévost (2011); Choo and Borja (2015); Sun et al. (2013); Berger et al. (2015); Rodrigo et al. (2018); Honório et al. (2018). However, the study of stabilization procedures addressing multiphase problems is still incipient, with just a few studies available Wan,Jing (2002); Truty and Zimmermann (2006).
This paper proposes a new stabilization technique in which the balance of mass equations are supplemented with stabilizing flux terms. The resulting technique mimics the LPJ stabilization Silvester and Kechkar (1990); Silvester (1994); Elman et al. (2005) in its basic design, but with suitable extensions to handle the poromechanical system of interest here. The additional stabilization terms are dependent on a stabilization parameter that must be well chosen to suppress instabilities while not compromising solution accuracy. We identify an optimal value for this parameter using an analysis of the eigenvalue distribution of the macroelement Schur complement matrix. The resulting method is simple to implement and preserves the underlying sparsity pattern of the original discretization. Another appealing feature of the method is that mass is exactly conserved on macroelements.
The paper is organized as follows. The governing equations and discretization scheme are introduced in Section 2 and 3. In Section 4 we examine the behavior of this model in the undrained limit, in order to identify the source of spurious modes. To fix this deficiency, our stabilization scheme is detailed in Section 5. The resulting approach both treats spurious pore pressure oscillations and improves the conditioning of the system matrices. This is demonstrated through numerical examples presented in Section 6. Finally, concluding remarks are given in Section 7.
2 Governing Equations
We consider a multiphase poroelastic problem in which two immiscible fluids fill the voids of the porous, deformable solid skeleton. We focus on a displacementsaturationpressure formulation, ignoring dynamic and nonisothermal effects. We further neglect capillary forces, meaning, the wetting and nonwetting fluid phases have equal pressure inside the pore. This simplification is common in many reservoirscale simulations Ertekin et al. (2001). The methods described here may be readily extended to include more sophisticated formulations.
The porous medium occupies a domain over time interval . The unknown fields are the displacement of the solid , the wetting fluid saturation , and the fluid pressure . The initial/boundary value problem is governed by a linear momentum balance for the mixture and two mass balance equations for the wetting ( and nonwetting () fluids, respectively:
(1a)  
(1b)  
(1c) 
In Eq. (1a), the effective Cauchy stress depends on the symmetric gradient of the displacement field as
(2) 
where
is the drained elasticity tensor. Biot’s coefficient
may be calculate from , the drained skeleton bulk modulus, and , the intrinsic bulk modulus of the solid phase. The mixture density is related to individual the phase densities—denoted by for and —via the relationship(3) 
where is the porosity. Porosity changes are related to solid deformation and fluid pressure changes as
(4) 
which introduces a coupling between the momentum and mass balance equations. Each fluid phase requires a density model , such as the simple linear model
(5) 
with phase bulk modulus and reference density at a reference pressure .
In Eq. (1b) and Eq. (1c), is the mass per unit volume for the fluid phase , with
(6a)  
(6b) 
The source terms are used to model well sources for injection and production of fluids, using a Peaceman well model Peaceman (1978); Chen et al. (2006). The mass flux is linked to the pore pressure field via the generalized Darcy’s law as
(7) 
The constitutive relation in (7) defines the volumetric flux using the phase mobility , the viscosity , and the relative permeability . Specific relationships for viscosity and relative permeability must be defined for the fluids and porous medium under consideration. Additionally, represents the absolute permeability tensor, the gravitational acceleration, and the elevation above a datum.
The domain boundary is decomposed into regions where Dirichlet and Neumann boundary conditions are specified, denoted by for the momentum balance and for the mass balances. These divisions follow the overlap restriction . Specifically,
(8a)  
(8b)  
(8c)  
(8d)  
(8e)  
(8f) 
where the boundary conditions prescribing displacement (8a), total traction (8b), pore pressure (8c), wetting phase saturation (8d), wetting phase mass flux (8e) and nonwetting phase mass flux (8f) are given. Here,
denotes the outer normal vector. Homogeneous conditions on the displacement and external fluxes were chosen here to simplify some notations below, but these can be easily relaxed.
Initial conditions are specified as
(9a)  
(9b)  
(9c) 
Note that the singlephase poromechanics model arises as a subcase of these general equations if one fixes either or . In this case, only one mass balance equation is required.
Clearly a number of modeling and constitutive assumptions are embedded in the formulation described above, but it remains a useful approximation for many subsurface applications. This formulation also contains many of the salient mathematical features that may be encountered in other models used in practice.
3 Discrete Formulation
Figure 1 illustrates the basic geometry under consideration. The domain is partitioned into a computational mesh made of nonoverlapping elements such that . Every element face is assigned a unique unit normal vector . For our stabilization procedure, we further assume that these elements may be grouped into macroelements consisting of patches of eight hexahedra in 3D or four quadrilaterals in 2D. This configuration is readily achieved by beginning with a coarse version of the mesh and applying one level of structured refinement. The discretization of the governing equations (1) is then obtained using a combined finiteelement/finitevolume approach. An extensive description of this formulation is presented in White et al. (2018). Here, we briefly summarize the key components, but refer the interested reader there for a more complete exposition.
Time integration relies on a fully implicit backward Euler scheme, with the time interval divided into subintervals of length . We will use the notation for other timedifferenced quantities as well. For the spatial discretization, we introduce the specific spaces
(10a)  
(10b) 
where and are the space of continuous and square Lebesgueintegrable functions on , respectively. denotes the space of linear functions (trilinear in 3D or bilinear in 2D) and the space of constant functions on a given element .
The discrete weak form of (1) reads: Find such that for time step
(11a)  
(11b)  
(11c) 
where are discrete test functions. The symbol denotes the jump of a quantity across face in . For an internal face, , with and the restriction of on cells and sharing , respectively. For a face belonging to the domain boundary, the jump expression reduces to . The term denotes a discrete mass flux, i.e. . These are computed using a standard twopoint flux approximation scheme, with upwinding of the phase density and mobility Aziz (1979):
(12) 
Here, is the transmissibility coefficient for the face, which is computed knowing the mesh geometry and permeability. The jump in the elevation datum should be understood as the difference in the coordinate of the respective cell centroids.
The unknown fields are interpolated as
(13a)  
(13b)  
(13c) 
with and bases for and , respectively. are nodal values of the displacement components, while and are cellcentered values for the saturation and pressure fields. An identical basis is introduced for the test functions.
The fully discrete system of equations at time is then obtained by introducing these bases into the weak form Eqs. (11a)(11c) and applying the standard finite element procedure. This leads to a set of algebraic equations for the unknown degreesoffreedom , and . These degreesoffreedom are assembled in an algebraic vector , , . The nonlinear system of equations is assembled in a residual vector
(14) 
The nonlinearity of the system here results from nonlinear constitutive behavior embedded in the relative permeability relationship . Furthermore, this formulation is general enough to accommodate other nonlinear constitutive relationships for the various solid and fluid components. A Newton iteration scheme is used to solve this system, which requires the linearization of the threefield problem. The linearized problem is defined by a Jacobian system with a 3 3 block structure of the form
(15) 
where is the Jacobian matrix, with , , and the Newton search directions for each field. The superscript stands for the Newton iteration count. Full expressions for each elemental subvector of and submatrix of are reported in White et al. (2018). We note that the solution of this linear system is typically the most expensive component of a fullyimplicit code, and good solver performance is therefore essential.
4 Incompressibility
It may not be immediately apparent that the system (15) may be subject to an infsup condition on its solvability. Indeed, for many problem configurations the discrete system is perfectly wellposed. Instabilities can arise, however, when two conditions are satisfied:

during undrained loading, i.e. as either or ;

when the solid and fluid phases approach incompressibility, i.e. for .
Note that it is sufficient to merely approach these limits, a situation that occurs frequently in practice. This is particularly true at early simulation times, when small time steps are often required to resolve rapidly evolving solution fields. Liquid and solid compressibilities are also often small for many geologic systems.
To highlight the origin of difficulties, we first revisit the continuum governing equations assuming the conditions above are exactly satisfied. In this case, several relationships simplify, in particular
(16a)  
(16b)  
(16c)  
(16d)  
(16e) 
We set here under the assumption that these source terms represent wells, which cannot inject when permeability goes to zero. The mass balance equations (1b)(1c) reduce to
(17a)  
(17b) 
Adding these two equations implies , and therefore . The reduced system of governing equations is therefore
(18a)  
(18b) 
with . We observe that under these conditions the solid deformation field must satisfy a divergence constraint condition, while the saturation field becomes fixed in time at its initial conditions. This result is physically intuitive. If the fluids can neither flow nor compress, they will not allow the solid skeleton to deform volumetrically, nor is there a mechanism for saturations to evolve. The result is a twofield problem only in displacement and pressure. With some additional manipulations one can show that this set of governing equations is equivalent to Stokes’ equations.
It is also instructive to perform the same exercise for the algebraic system (15). When phase compressibility is zero and undrained conditions are reached, the Jacobian matrix becomes
(19) 
with the following expressions for the individual blocks:
(20a)  
(20b)  
(20c)  
(20d)  
(20e)  
(20f) 
Now imagine that we add the third block row to the second, scaled by their (constant) densities. We observe that
(21a)  
(21b) 
We may therefore identify a subsystem for displacements and pressures that is uncoupled from the saturation field,
(22) 
One would arrive at the same system through a direct discretization of the reduced governing equations above. It is clear that this matrix is in saddlepoint form, and that the spaces chosen for the pressure and displacement fields must satisfy an infsup compatibility condition to ensure has full rank. For our chosen interpolation, however, this is not the case. As the incompressible limit is approached, —and equivalently, —may contain nearsingular modes that can express as spurious oscillations in the pressure field.
5 Stabilized Formulation
As a fix for this difficulty, we propose a simple modification to the way the discrete fluxes are treated in Eqs. (11b)–(11c). As described earlier, the mesh is decomposed into macroelements. Let denotes the union of all faces that lie interior to any macroelement. That is, any two cells connected across a face are members of the same parent macroelement.
For any face in , we augment the physical flux with an additional stabilization flux . For a given time increment , we replace
(23) 
where for each phase the stabilization flux is a function of an interelement jump across the face, scaled by particular constants,
(24a)  
(24b) 
These artificial fluxes will be used to control spurious pressure modes associated with nonphysical pressure jumps across faces. Here, is a stabilization parameter and is the average volume of the child elements in the macroelement. The remaining terms are the upwinded density and phase saturation for the respective phases. Note that these are lagged in time to simplify the linearization, as a lagged approximation of these quantities is usually sufficient for stabilization purposes. We will discuss the choice of stabilization constant below, which is critical to success.
This flux form is quite similar to the physical flux computation Eq. (12), so it may be readily added to an existing facebased assembly loop. Any addition of artificial fluxes, however, will break the elementwise mass conservation property of the underlying finite volume scheme. Because these artificial fluxes are only added to internal faces of the macroelement, however, exact mass conservation is still preserved on the macroelement level.
When assembled, these flux terms add additional entries to two blocks of the system matrix,
(25) 
where the stabilizing entries are assembled facewise for any as
(26a)  
(26b) 
In the incompressible limit, these contributions will not vanish, so that
(27) 
In practice, we always solve the threefield problem. However, it is instructive to apply the same reduction procedure as before for the incompressible limit. This leads to a reduced system
(28) 
where
(29) 
Thus, the original saddlepoint system is modified so that new entries appear in the lowerrighthand block. For each face , this matrix contains contributions
(30) 
Because of the macroelement construction, the resulting matrix is extremely sparse. In 3D, it is blockdiagonal with one block for each macroelement in the mesh, with entries
(31) 
where is the average volume of the child elements in the macroelement. Note that in our mesh geometry, for any face interior to the macroelement, with the face area and the distance between the centroids of the neighboring cells. By construction, this pattern preserves the underlying twopoint flux approximation (TPFA) stencil adopted in the original finite volume scheme, and will not cause any fillin. Note that the matrices and will have the same sparsity pattern as , though their entries are weighted by the local saturation and densities at the faces.
For the reduced system, the stabilization contribution has a very similar form to the Local Pressure Jump (LPJ) stabilization as originally formulated for solving the Stokes equation Silvester and Kechkar (1990). Indeed, this basic idea of using interelement pressure jumps to control spurious modes inspired the method proposed here. There are two key differences for the multiphase poromechanics application, however:

In the threefield formulation, a separate contribution is made to each mass balance equation, weighted by appropriate phase density and saturation;

The optimal stabilization constant will differ due to the nature of the underlying equations.
We still have to address this question of what is an appropriate value for the stabilization parameter . As proposed in Silvester (1994), good candidates for can be determined by examining the spectrum of the Schur complement matrix,
(32) 
which corresponds to a further block reduction of to a pressureonly system. We focus on a patch test involving a single macroelement, with rigid and impermeable boundary conditions (Figure 2). If stability can be demonstrated for a single macroelement, theoretical results in Boland and Nicolaides (1983); Stenberg (1984) prove stability for discretizations on arbitrary grids constructed by “gluing together” stable macroelements.
The resulting Schurcomplement is rank deficient without stabilization. In 3D, there are eight pressure unknowns in the cells but only three displacement components at the central node. Similar to Elman et al. (2005)
, the eigenvalues and eigenvectors of this system may be readily computed. Let the cells have edge lengths
, , and . The element volume is , and each face has area , , or . The resulting eigenvalues are(33a)  
(33b)  
(33c)  
(33d)  
(33e)  
(33f) 
where and are the two Lamé parameters characterizing the elastic mechanical response.
When no stabilization is applied to the macroelement, that is when , five out of eight eigenvalues are zero. The null eigenvectors include one constant pressure mode (associated to ) and four spurious pressure modes (associated to –). The constant mode is expected here since the boundary conditions only determine the pressure solution up to an arbitrary constant. We see that for , the stabilization will remove the spurious pressure modes from the null space of . Unfortunately, the choice of will also impact the physical modes associated with –.
Symbol  Parameter  Drained  Undrained  Modified  Units 

Time step  s  
Final time  s  
Biot coefficient  1  1  1    
Young’s modulus  2.5  Pa  
Poisson ratio  0.1  0.1  0.25  
Permeability  
Viscosity  Pas  
Density 
For stability, all eigenvalues must remain bounded away from zero except for the constant mode . Taking too large, however, will corrupt the physical solution and compromise the approximation. A good choice to balance these competing priorities is to choose a that minimizes the condition number . Considering a regular cube with equal sides , the minimal condition number is retrieved for any stabilization parameter lying within the range
(34) 
Figure 3 illustrates how the condition number varies depending on the stabilization parameter. We can see that the condition number attains its minimal value of within the range prescribed in Eq.(34). Within this range, neither extremal eigenvalue actually depends on . Therefore, a reasonable choice for the stabilization parameter is
(35) 
In order to explore the sensitivity of numerical solutions to this constant, it is convenient to present results in terms of the ratio , i.e. the ratio of a given stabilization to the recommended value , with being “optimal”.
We emphasize that the recommended value here is based on the analysis of a fullyincompressible, homogeneous, cubic macroelement. It therefore only depends on two elastic constants. A more precise analysis may lead to a stabilization constant that additionally depends on the Biot coefficient, fluid compressibilities, macroelement heterogeneity, and so forth. For example, when
a better estimate would be
.6 Numerical Examples
We begin with a few singlephase examples () to demonstrate the performance of the method in a simpler setting. We then conclude with a multiphase demonstration for a benchmark reservoir simulation problem. These numerical experiments were implemented using Geocentric, a simulation framework for computational geomechanics that relies heavily on finite element infrastructure from the the deal.ii library Bangerth et al. (2007).
6.1 Singlephase Examples


For the singlephase examples, we consider several variants of Barry and Mercer’s problem for a twodimensional, poroelastic medium Barry and Mercer (1999). Parameter values are summarized in Table 1. The first test illustrates that the proposed stabilization scheme does not compromise solution accuracy under drained conditions, when no instabilities are expected. Two subsequent examples confirm the effectiveness of the scheme in suppressing spurious pressure oscillations under undrained conditions.
6.1.1 Drained BarryMercer
Barry and Mercer Barry and Mercer (1999) provide an analytical solution for a two dimensional problem. The problem setup consists of a square domain subjected to a periodic point source given as
(36) 
Here, and denote the two elastic Lamé parameters, while and are the isotropic absolute permeability and viscosity, respectively. The point source is located at , and indicates a Dirac function. All sides of the computational domain are constrained with zero pressure and zero tangential displacement boundary conditions, as depicted in Figure 4. The simulation parameters provided in Table 1 (drained conditions) are the same as those used in Fu (2019); Boffi et al. (2016); Rodrigo et al. (2016); Phillips and Wheeler (2007). Note that the time step and final time correspond to a normalized time of and , respectively.
Figure 5 shows the resulting pressure profile at the final time along a vertical line through the source point. Both the analytical solution and the numerical solution for different mesh refinements are shown. Good agreement between the exact and computed results is also indicated by Figure 6, which shows convergence behavior of the error for the pressure solution for both the stabilized and unstabilized formulations. One observes a linear and essentially identical error behavior for both, indicating that the macroelement stabilization does not compromise solution accuracy in regimes where it is not strictly needed.
6.1.2 Undrained BarryMercer
The goal of this example is to show the effectiveness of proposed stabilization scheme in treating nonphysical pressure oscillations. These spurious pressure oscillations appear in the limit of low permeability or fast loading rates. As in Rodrigo et al. (2016); Fu (2019); Boffi et al. (2016), we use the same simulation parameters as the previous section, but we decrease the value of the permeability to and perform only one time step of s.
Figure 7 shows the pressure contour plot for the uniformly discretized domain with . The pressure field exhibits mild oscillations close to the sourcepoint. These oscillations are eliminated when using the proposed stabilization technique.
It is also interesting to examine the conditioning of the stabilized system and its impact on iterative solver performance. To do so, we define a scaled Schurcomplement matrix,
(37) 
where is the mass matrix on the pressure space. For the space, this is a diagonal matrix with entries corresponding to the element volumes. The inverse of this diagonal matrix introduces a volume scaling that allows eigenvalues to be properly compared across different mesh refinement levels.
Figure 8 presents the condition number of the for various choices of stabilization constant. The minimum is achieved close to the recommended value that was inferred from the single macroelement analysis. Furthermore, the removal of nearsingular modes from the system matrix has a dramatic impact on iterative solver performance. Figure 9 presents the number of Krylov iterations to convergence needed for different mesh refinements. For low values of the stabilization constant, the nearsingular modes cause a dramatic degradation in the linear solver performance.
6.1.3 Modified Undrained BarryMercer
This last variant of the Barry Mercer’s problem tests the efficacy of the stabilization when even more severe pressure oscillations are present. The setup is based on Haga et al. (2012), where the only difference with the previous setup is that the source point is switched from ratecontrolled to pressurecontrolled. The applied pressure varies according to
(38) 
with . The other simulation parameters follow Haga et al. (2012) and are listed in Table 1. An exact solution is not available in this case. However, at early times ( s) we can infer that the whole domain should have zero pressure except for the cell where the pressure is enforced. Figure 10 illustrates the solution obtained from the unstabilized and the stabilized schemes for a domain discretized with cell size . As expected, the stabilization considering eliminates the wild oscillations. Table 2 reports the extremal eigenvalues and condition number of the for solutions with and without stabilization. One observes that as the mesh is refined the minimal eigenvalue and the condition number converge to an asymptotic value different than zero and infinity, respectively, only when using the stabilized scheme. Figure 11 presents the condition number of the matrix as a function of the stabilization constant. Once again, the minimimal condition number is attained when . Table Figure 12 shows the resulting improvement in Krylov convergence.
6.2 Multiphase Example
Parameter  Units  Value 

Porosity:  
Highperm zone  –  0.20 
Lowperm zone  –  0.05 
Permeability:  
Highperm zone  mD  1000 
Lowperm zone  mD  1 
Relative perm:  
Residual wetting sat.  –  0.2 
Residual nonwetting sat.  –  0.2 
Wetting Fluid:  
Reference density  kg/m  1035 
Bulk modulus  MPa  
Viscosity  cP  0.3 
Nonwetting Fluid:  
Reference density  kg/m  863 
Bulk modulus  MPa  
Viscosity  cP  3.0 
Rock:  
Young’s modulus  MPa  5000 
Biot coefficient  –  1 
Grain density  kg/m  2650 
Well control:  
Injection BHP  MPa  
Production BHP  MPa  
Ramp time  day  1 
Well radius  m  0.1524 
Skin factor  –  0 
Timestepping:  
Initial  day  0.0001 
Maximum  day  1 
End time  day  100 
Solver tolerances:  
Newton  –  10 
Krylov  –  10 



Lastly, we consider a multiphase poromechanics example. The test problem is based on the staircase benchmark problem originally presented in White et al. (2018). Figure 13 provides an illustration of the problem geometry. The domain contains two regions, a highpermeability channel and a lowpermeability host rock. The highpermeability channel winds its way in a spiral, staircase fashion from an upper injection well to a lower production well. This spiral geometry is obviously artificial, but it introduces very strong coupling between the displacement, pressure, and saturation fields. Water is injected at the upper well, while both fluids may be produced from the lower well. These wells have a bottomhole pressure (BHP) control. The injector (producer) is ramped up to MPa ( MPa) overpressure over one day, and then held at a constant pressure. All problem parameters are given in Table 3. We have set the fluids to be incompressible to accentuate any instabilities in the formulation. Specifics regarding the well control, relative permeability model, and mechanical boundary conditions may be found in White et al. (2018). At the beginning of the simulation, the initial time step is day. This time step is then doubled every step until a maximum time step of day is reached. The whole simulation is run for 100 days. Note that small time steps are the most problematic from a stability point of view.
Figure 13(c) presents pressure and saturation snapshots for this simulation, using both an unstabilized and stabilized formulation (with ). At the first time step ( day) checkerboard oscillations are apparent in the pressure field for the unstabilized formulations. Note that we have truncated the colorbar, cutting off the peak pressures, in order to accentuate these oscillations visually. The stabilization successfully suppresses this checkerboarding. At the end of the simulation ( day), we see that the stabilized and unstabilized formulations produce essentially identical results. The addition of the artificial flux terms does not compromise overall solution quality, with the saturation field being advected the same distance along the highperm channel in both cases.
It is interesting to examine the linear solver behavior at early times in the simulation (Figure 15). At each Newton step of the nonlinear solver, a preconditioned GMRES iteration is used to solve the Jacobian system. The preconditioner we use is the multistage preconditioner described in White et al. (2018), which in general provides excellent convergence behavior for this class of problem. In the first day of simulation time, however, we see a substantial degradation in solver performance using the unstabilized formulation. This is a direct results of the presence of nearsingular modes, to which Krylovbased solvers are extremely sensitive. With the addition of stabilization, however, this problem is completely removed and GMRES once again exhibits excellent convergence. Later in the simulation, as grows, the physical fluxes between elements grow and even the unstabilized formulation becomes intrinsically stable. As a result, we see the solver convergence behavior merge at later times for the two formulations.
7 Conclusion
In this work, we have presented a stabilized formulation for discretizations of single and multiphase poromechanics. The stabilization is achieved by adding artificial flux terms to faces interior to macroelements. We have also identified an appropriate value for the stabilization parameter based on an eigenvalue analysis of an incompressible macroelement patch test. The stabilization is easy to implement in existing codes, and does not change the underlying sparsity pattern of the finite volume stencil. While exact mass conservation on individual elements is sacrificed, exact mass conservation on macroelements is retained. We have demonstrated, through a number of single and multiphase examples, that the method is effective in practice. It can suppress spurious oscillations and also prevent unwanted degradation in iterative solver convergence in the presence of nearsingular modes. The latter is a critical issue for largescale simulations of geosystems.
While the discussion here has been limited to discretizations, a similar approach can likely be used for other unstable interpolation pairs involving piecewise constant pressure approximations–e.g. on more general hexahedral or tetrahedral meshes.
Acknowledgements
Funding for JTC and JAW was provided by Total S.A. through the FCMAELSTROM Project. JTC also acknowledges financial support provided by the Brazilian National Council for Scientific and Technological Development (CNPq) and the the John A. Blume Earthquake Engineering Center. RIB was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Geosciences Research Program, under Award Number DEFG0203ER15454. The authors wish to thank Nicola Castelletto for helpful discussions. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DEAC5207NA27344.
References
 General theory of threedimensional consolidation. Journal of Applied Physics 12 (2), pp. 155–164. External Links: Document Cited by: §1.
 Petroleum reservoir simulation. Vol. 476, Applied Science Publishers, London. Cited by: §3.
 The finite element method with Lagrangian multipliers. Numerische Mathematik 20 (3), pp. 179–192. External Links: ISSN 09453245, Document Cited by: §1.
 deal.II – a General Purpose Object Oriented Finite Element Library. ACM Trans. Math. Softw. 33 (4), pp. 24/1–24/27. Cited by: §6.
 Exact solutions for twodimensional timedependent flow and deformation within a poroelastic medium. Journal of applied mechanics 66 (2), pp. 536–540. External Links: Document Cited by: §6.1.1, §6.1.
 Stabilized lowestorder finite element approximation for linear threefield poroelasticity. SIAM Journal on Scientific Computing 37 (5), pp. A2222–A2245. External Links: Document Cited by: §1.
 Stabilization of loworder mixed finite elements for the Stokes equations. SIAM Journal on Numerical Analysis 44 (1), pp. 82–101. External Links: Document Cited by: §1.
 A nonconforming highorder method for the Biot problem on general meshes. SIAM Journal on Scientific Computing 38 (3), pp. A1508–A1537. External Links: Document Cited by: §6.1.1, §6.1.2.
 Stability of finite elements under divergence constraints. SIAM Journal on Numerical Analysis 20 (4), pp. 722–731. External Links: ISSN 00361429, Document Cited by: §5.
 Multiphysics hillslope processes triggering landslides. Acta Geotechnica 7 (4), pp. 261–269. External Links: ISSN 18611133, Document Cited by: §1.
 Continuum deformation and stability analyses of a steep hillside slope under rainfall infiltration. Acta Geotechnica 5 (1), pp. 1–14. External Links: ISSN 18611133, Document Cited by: §1.
 On the stabilization of finite element approximations of the Stokes equations. In Efficient Solutions of Elliptic Systems: Proceedings of a GAMMSeminar Kiel, January 27 to 29, 1984, W. Hackbusch (Ed.), pp. 11 – 19. External Links: ISBN 9783663141693, Document Cited by: §1.
 On the existence, uniqueness and approximation of saddlepoint problems arising from Lagrangian multipliers. R.A.I.R.O. Analyse Numérique 8 (R2), pp. 129–151. External Links: Document Cited by: §1.
 A unified stabilized method for Stokes and Darcy’s equations. Journal of Computational and Applied Mathematics 198 (1), pp. 35 – 51. External Links: ISSN 03770427, Document Cited by: §1.
 Numerical limit analysis of threedimensional slope stability problems in catchment areas. Acta Geotechnica 11 (6), pp. 1369–1383. External Links: ISSN 18611133, Document Cited by: §1.
 On the soil water characteristic curves of poorly graded granular materials in aqueous polymer solutions. Acta Geotechnica 13 (1), pp. 103–116. External Links: ISSN 18611133, Document Cited by: §1.
 Accuracy and convergence properties of the fixedstress iterative solution of twoway coupled poromechanics. International Journal for Numerical and Analytical Methods in Geomechanics 39 (14), pp. 1593–1618. External Links: Document Cited by: §1.
 Computational Methods for Multiphase Flows in Porous. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. Cited by: §2.
 Stabilized mixed finite elements for deformable porous media with double porosity. Computer Methods in Applied Mechanics and Engineering 293, pp. 131 – 154. External Links: ISSN 00457825, Document Cited by: §1.
 Enriched Galerkin finite elements for coupled poromechanics with local mass conservation. Computer Methods in Applied Mechanics and Engineering 341, pp. 311–332. External Links: Document Cited by: §1.
 Poromechanics. John Wiley & Sons, Ltd. External Links: Document Cited by: §1.
 Theory of porous media: highlights in historical development and current state. Springer Science & Business Media. External Links: Document Cited by: §1.
 A stabilized finite element method for the stokes problem based on polynomial pressure projections. International Journal for Numerical Methods in Fluids 46 (2), pp. 183–201. External Links: Document Cited by: §1.
 Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Oxford University Press. Cited by: §1, §5.
 Basic applied reservoir simulation. Society of Petroleum Engineers, Richardson. Cited by: §2.
 Continuum hydrodynamics of dry granular flows employing multiplicative elastoplasticity. Acta Geotechnica 13 (5), pp. 1027–1040. External Links: ISSN 18611133, Document Cited by: §1.
 XFEMbased cohesive zone approach for modeling nearwellbore hydraulic fracture complexity. Acta Geotechnica 14 (2), pp. 377–402. External Links: ISSN 18611133, Document Cited by: §1.
 Iterative stabilization of the bilinear velocity–constant pressure element. International Journal for Numerical Methods in Fluids 10 (2), pp. 125–140. External Links: Document Cited by: §1.
 A highorder HDG method for the Biot’s consolidation model. Computers & Mathematics with Applications 77 (1), pp. 237 – 252. External Links: ISSN 08981221, Document Cited by: §6.1.1, §6.1.2.
 [30] Fully coupled elastoplastic hydromechanical analysis of unsaturated porous media using a meshfree method. International Journal for Numerical and Analytical Methods in Geomechanics. External Links: Document Cited by: §1.
 SPH approach for simulating hydromechanical processes with large deformations and variable permeabilities. Acta Geotechnica 13 (2), pp. 303–316. External Links: ISSN 18611133, Document Cited by: §1.
 On the causes of pressure oscillations in lowpermeable and lowcompressible porous media. International Journal for Numerical and Analytical Methods in Geomechanics 36 (12), pp. 1507–1522. External Links: Document Cited by: §6.1.3.
 A stabilized elementbased finite volume method for poroelastic problems. Journal of Computational Physics 364, pp. 49 – 72. External Links: ISSN 00219991, Document Cited by: §1, §1.
 A new finite element formulation for computational fluid dynamics: V. Circumventing the BabuškaBrezzi condition: a stable PetrovGalerkin formulation of the Stokes problem accommodating equalorder interpolations. Computer Methods in Applied Mechanics and Engineering 59 (1), pp. 85 – 99. External Links: ISSN 00457825, Document Cited by: §1.
 A new finite element formulation for computational fluid dynamics: VII. the Stokes problem with various wellposed boundary conditions: symmetric formulations that converge for all velocity/pressure spaces. Computer Methods in Applied Mechanics and Engineering 65 (1), pp. 85 – 96. External Links: ISSN 00457825, Document Cited by: §1.
 Fluiddriven transition from damage to fracture in anisotropic porous media: a multiscale XFEM approach. Acta Geotechnica. External Links: ISSN 18611133, Document Cited by: §1.
 Stability, accuracy, and efficiency of sequential methods for coupled flow and geomechanics. SPE Journal 16 (02), pp. 249–262. External Links: Document, ISBN 1086055X Cited by: §1.
 Probabilistic multiscale characterization and modeling of organicrich shale poroelastic properties. Acta Geotechnica 13 (4), pp. 781–800. External Links: ISSN 18611133, Document Cited by: §1.
 XFEM, strong discontinuities and secondorder work in shear band modeling of saturated porous media. Acta Geotechnica 13 (6), pp. 1249–1264. External Links: ISSN 18611133, Document Cited by: §1.
 Explicit meshfree solution for large deformation dynamic problems in saturated porous media. Acta Geotechnica 13 (2), pp. 227–242. External Links: ISSN 18611133, Document Cited by: §1.
 Cellcentered finite volume discretizations for deformable porous media. International Journal for Numerical Methods in Engineering 100 (6), pp. 399–418. External Links: Document Cited by: §1.
 Stable cellcentered finite volume discretization for Biot equations. SIAM Journal on Numerical Analysis 54 (2), pp. 942–968. External Links: Document Cited by: §1.
 A computational model for dynamic strain localization in unsaturated elastoviscoplastic soils. International Journal for Numerical and Analytical Methods in Geomechanics 43 (1), pp. 138–165. External Links: Document Cited by: §1.
 Interpretation of wellblock pressures in numerical reservoir simulation. Soc. Pet. Eng., AIME J. 18 (3), pp. 183–194. External Links: Document Cited by: §2.
 Fully coupled hydraulic fracture simulation using the improved element partition method. International Journal for Numerical and Analytical Methods in Geomechanics 43 (1), pp. 441–460. External Links: Document Cited by: §1.
 A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discreteintime case. Computational Geosciences 11 (2), pp. 145–158. External Links: ISSN 15731499, Document Cited by: §6.1.1.
 A multigrid version of a simple finite element method for the Stokes problem. Mathematics of Computation 45 (171), pp. 1–14. External Links: ISSN 00255718, 10886842 Cited by: §1.
 Highorder ADE scheme for solving the fluid diffusion equation in nonuniform grids and its application in coupled hydromechanical simulation. International Journal for Numerical and Analytical Methods in Geomechanics 42 (16), pp. 1976–2000. External Links: Document Cited by: §1.
 Stabilization procedures in coupled poromechanics problems: a critical assessment. International Journal for Numerical and Analytical Methods in Geomechanics 35 (11), pp. 1207–1225. External Links: Document Cited by: §1.
 Twoway coupling in reservoir geomechanical models: vertexcentered Galerkin geomechanical model cellcentered and vertexcentered finite volume reservoir models. International Journal for Numerical Methods in Engineering 98 (8), pp. 612–624. External Links: Document Cited by: §1.
 Stability and monotonicity for some discretizations of the Biot’s consolidation model. Computer Methods in Applied Mechanics and Engineering 298, pp. 183–204. External Links: Document Cited by: §6.1.1, §6.1.2.
 New stabilized discretizations for poroelasticity and the Stokes’ equations. Computer Methods in Applied Mechanics and Engineering 341, pp. 467–484. Cited by: §1.
 Thermoplasticity and strain localization in transversely isotropic materials based on anisotropic critical state plasticity. International Journal for Numerical and Analytical Methods in Geomechanics 40 (18), pp. 2423–2449. External Links: Document Cited by: §1.
 A coupled reservoir and geomechanical simulation system. SPE Journal 3 (03), pp. 219–226. External Links: Document Cited by: §1.
 The effect of stress boundary conditions on fluiddriven fracture propagation in porous media using a phasefield modeling approach. International Journal for Numerical and Analytical Methods in Geomechanics 43 (6), pp. 1316–1340. External Links: Document Cited by: §1.
 Stabilised bilinearconstant velocitypressure finite elements for the conjugate gradient solution of the Stokes problem. Computer Methods in Applied Mechanics and Engineering 79 (1), pp. 71 – 86. External Links: ISSN 00457825, Document Cited by: §1, §1, §5.
 Optimal low order finite element methods for incompressible flow. Computer methods in applied mechanics and engineering 111 (34), pp. 357–368. External Links: ISSN 00457825, Document Cited by: §1, §5.
 Multiscale finite volume method for finitevolumebased simulation of poroelasticity. Journal of Computational Physics 379, pp. 309 – 324. External Links: ISSN 00219991, Document Cited by: §1.
 Localized failure in unsaturated soils under nonisothermal conditions. Acta Geotechnica 13 (1), pp. 73–85. External Links: ISSN 18611133, Document Cited by: §1.
 Analysis of mixed finite element methods for the Stokes problem: a unified approach. Mathematics of Computation 42 (165), pp. 9–23. External Links: ISSN 00255718, 10886842, Document Cited by: §5.
 A stabilized assumed deformation gradient finite element formulation for strongly coupled poromechanical simulations at finite strain. International Journal for Numerical and Analytical Methods in Geomechanics 37 (16), pp. 2755–2788. External Links: Document Cited by: §1.
 Groundwater pumping and land subsidence in the EmiliaRomagna coastland, Italy: modeling the past occurrence and the future trend. Water Resources Research 42 (1), pp. . External Links: Document Cited by: §1.
 Stabilized mixed finite element formulations for materially nonlinear partially saturated twophase media. Computer Methods in Applied Mechanics and Engineering 195 (13), pp. 1517 – 1546. External Links: ISSN 00457825, Document Cited by: §1.
 Comparison of geomechanical deformation induced by megatonnescale CO2 storage at Sleipner, Weyburn, and In Salah. Proceedings of the National Academy of Sciences 110 (30), pp. E2762–E2771. External Links: Document, ISSN 00278424 Cited by: §1.
 Stabilized finite element methods for coupled geomechanics and multiphase flow. Ph.D. Thesis, Stanford University. Cited by: §1.
 A coupled 3dimensional bonded discrete element and lattice Boltzmann method for fluidsolid coupling in cohesive geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics 42 (12), pp. 1405–1424. External Links: Document Cited by: §1.
 A dual mortar contact method for porous media and its application to claycore rockfill dams. International Journal for Numerical and Analytical Methods in Geomechanics 43 (9), pp. 1744–1769. External Links: Document Cited by: §1.
 A twostage preconditioner for multiphase poromechanics in reservoir simulation. arXiv preprint arXiv:1812.05540. Cited by: §1, §3, §3, §6.2, §6.2.
 Blockpartitioned solvers for coupled poromechanics: a unified framework. Computer Methods in Applied Mechanics and Engineering 303, pp. 55–74. External Links: Document Cited by: §1.
 Geomechanical behavior of the reservoir and caprock system at the In Salah CO2 storage project. Proceedings of the National Academy of Sciences 111 (24), pp. 8747–8752. External Links: Document Cited by: §1.
 Stabilized loworder finite elements for coupled soliddeformation/fluiddiffusion and their application to fault zone transients. Computer Methods in Applied Mechanics and Engineering 197 (49), pp. 4353 – 4366. External Links: ISSN 00457825, Document Cited by: §1.
 A 2D coupled hydrothermal model for the combined finitediscrete element method. Acta Geotechnica 14 (2), pp. 403–416. External Links: ISSN 18611133, Document Cited by: §1.
 FDEMTH3D: A threedimensional coupled hydrothermal model for fractured rock. International Journal for Numerical and Analytical Methods in Geomechanics 43 (1), pp. 415–440. External Links: Document Cited by: §1.
 On the preferential flow patterns induced by transverse isotropy and nonDarcy flow in double porosity media. Computer Methods in Applied Mechanics and Engineering 353, pp. 570 – 592. External Links: ISSN 00457825, Document Cited by: §1.
 Deformation and strength of transversely isotropic rocks. In Desiderata Geotechnica, W. Wu (Ed.), Cham, pp. 237–241. External Links: ISBN 9783030149871, Document Cited by: §1.
 On the strength of transversely isotropic rocks. International Journal for Numerical and Analytical Methods in Geomechanics 42 (16), pp. 1917–1934. External Links: Document Cited by: §1.
 Micromechanism of the diffusion of cementbased grouts in porous media under two hydraulic operating conditions: constant flow rate and constant pressure. Acta Geotechnica 14 (3), pp. 825–841. External Links: ISSN 18611133, Document Cited by: §1.
 Reservoir geomechanics. Cambridge University Press. External Links: Document Cited by: §1.
Comments
There are no comments yet.