A macroelement stabilization for multiphase poromechanics

Strong coupling between geomechanical deformation and multiphase fluid flow appears in a variety of geoscience applications. A common discretization strategy for these problems is a continuous Galerkin finite element scheme for the momentum balance equations and a finite volume scheme for the mass balance equations. When applied within a fully-implicit solution strategy, however, this discretization is not intrinsically stable. In the limit of small time steps or low permeabilities, spurious oscillations in the pressure field, i.e. checkerboarding, may be observed. Further, eigenvalues associated with the spurious modes will control the conditioning of the matrices and can dramatically degrade the convergence rate of iterative linear solvers. Here, we propose a stabilization technique in which the balance of mass equations are supplemented with stabilizing flux terms on a macroelement basis. The additional stabilization terms are dependent on a stabilization parameter. 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, despite the addition of artificial fluxes. The efficacy of the proposed technique is demonstrated with several numerical examples.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

12/05/2020

A finite-element framework for a mimetic finite-difference discretization of Maxwell's equations

Maxwell's equations are a system of partial differential equations that ...
06/18/2020

Well-posedness, discretization and preconditioners for a class of models for mixed-dimensional problems with high dimensional gap

In this work, we illustrate the underlying mathematical structure of mix...
09/29/2021

A Mass Conserving Mixed hp-FEM Scheme for Stokes Flow. Part III: Implementation and Preconditioning

This is the third part in a series on a mass conserving, high order, mix...
10/13/2019

A Conservative Finite Element Method for the Incompressible Euler Equations with Variable Density

We construct a finite element discretization and time-stepping scheme fo...
03/12/2020

A Class of Collocated Finite Volume Schemes for Incompressible Flow Problems

In this paper, we present a class of finite volume schemes for incompres...
12/27/2021

Robust approximation of generalized Biot-Brinkman problems

The generalized Biot-Brinkman equations describe the displacement, press...
04/27/2020

Clustering via torque balance with mass and distance

Grouping similar objects is a fundamental tool of scientific analysis, u...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 single-phase formulation is a straightforward sub-case. 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 fully-implicit time integration strategy is adopted, where all unknown fields are updated simultaneously in a monolithic manner White et al. (2016, 2018). A variety of finite-element and finite-volume 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 element-wise 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 Ladyzhenskaya-Babuška-Brezzi (LBB or inf-sup) condition Babuška (1973); Brezzi (1974). Unfortunately, the combination explored here is not intrinsically LBB-stable. 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. Near-zero 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 carefully-designed 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 Least-Squares 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 inter-elements 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 single-phase 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 displacement-saturation-pressure formulation, ignoring dynamic and non-isothermal effects. We further neglect capillary forces, meaning, the wetting and non-wetting fluid phases have equal pressure inside the pore. This simplification is common in many reservoir-scale 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 non-wetting () 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 non-wetting 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 single-phase 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:

Example mesh with nodal and cell-centered degrees-of-freedom. Each element is assigned to a parent macroelement.

Figure 1 illustrates the basic geometry under consideration. The domain is partitioned into a computational mesh made of non-overlapping 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 finite-element/finite-volume 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 time-differenced quantities as well. For the spatial discretization, we introduce the specific spaces

(10a)
(10b)

where and are the space of continuous and square Lebesgue-integrable 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 two-point 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 cell-centered 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 degrees-of-freedom , and . These degrees-of-freedom 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 three-field 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 sub-vector of and sub-matrix 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 fully-implicit 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 inf-sup condition on its solvability. Indeed, for many problem configurations the discrete system is perfectly well-posed. Instabilities can arise, however, when two conditions are satisfied:

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

  2. 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 two-field 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 sub-system 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 saddle-point form, and that the spaces chosen for the pressure and displacement fields must satisfy an inf-sup 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 near-singular 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 inter-element jump across the face, scaled by particular constants,

(24a)
(24b)

These artificial fluxes will be used to control spurious pressure modes associated with non-physical 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 face-based assembly loop. Any addition of artificial fluxes, however, will break the element-wise 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 face-wise for any as

(26a)
(26b)

In the incompressible limit, these contributions will not vanish, so that

(27)

In practice, we always solve the three-field 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 saddle-point system is modified so that new entries appear in the lower-right-hand block. For each face , this matrix contains contributions

(30)

Because of the macroelement construction, the resulting matrix is extremely sparse. In 3D, it is block-diagonal 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 two-point flux approximation (TPFA) stencil adopted in the original finite volume scheme, and will not cause any fill-in. 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 inter-element pressure jumps to control spurious modes inspired the method proposed here. There are two key differences for the multiphase poromechanics application, however:

  1. In the three-field formulation, a separate contribution is made to each mass balance equation, weighted by appropriate phase density and saturation;

  2. 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 pressure-only 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.

Figure 2: Single macroelement patch test geometry in 2D, with one pressure in each cell and two displacement components at the central node. The 3D macroelement is similar, involving eight pressures and three displacement components.

The resulting Schur-complement 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 .

Figure 3: Analytically computed condition number of as a function of the stabilization parameter for the single macroelement patch test. Here, with being the recommended value.
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
Table 1: Simulation parameters used in the three single-phase Barry-Mercer examples.

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 fully-incompressible, 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

.

In summary, the proposed stabilization method consists of adding the artificial flux terms in Eq. (24) to all macroelement-interior faces, weighted by the stabilization constant recommended in Eq. (35). We now explore the effectiveness of this simple method with several numerical examples.

6 Numerical Examples

We begin with a few single-phase 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 Single-phase Examples

Figure 4: Geometry and boundary conditions for the Barry-Mercer problem.
Figure 5: Pressure plot along the vertical line for the drained Barry-Mercer problem at using the stabilized formulation.
Figure 6: Convergence of the relative -error in pressure for the drained Barry-Mercer problem at , using both the unstabilized () and stabilized () formulations.
(a) Unstabilized scheme
(b) Stabilized scheme (with )
Figure 7: Pressure distribution for the undrained Barry-Mercer problem.
Figure 8: Computed condition numbers of for the undrained Barry-Mercer problem for various stabilization values .
Figure 9: Krylov iterations to convergence as a function of mesh refinement and stabilization constant for the undrained Barry-Mercer problem.
(a) Unstabilized scheme
(b) Stabilized scheme (with )
Figure 10: Pressure distribution for the modified undrained Barry and Mercer’s problem.
Number
of cells
Condition
number ()
8 x 8 2.51 e 0.330 131.56
16 x 16 5.35 e 0.333 621.85
32 x 32 1.19 e 0.333 2799.72
(a)
Number
of cells
Condition
number ()
8 x 8 0.222 0.539 2.421
16 x 16 0.224 0.546 2.437
32 x 32 0.225 0.548 2.438
(b)
Table 2: Extremal eigenvalues of the scaled Schur complement matrix for the modified undrained Barry-Mercer problem. Minimum and maximum eigenvalues as well as the condition number are presented as a function of mesh refinement. Two cases are considered: using (a) an unstabilized scheme and (b) the proposed stabilization.
Figure 11: Computed condition numbers of for the modified undrained Barry-Mercer problem for various stabilization values .
Figure 12: Krylov iterations as a function of the mesh refinement and the stabilization constant for the modified undrained Barry Mercer problem.

For the single-phase examples, we consider several variants of Barry and Mercer’s problem for a two-dimensional, 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 Barry-Mercer

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 Barry-Mercer

The goal of this example is to show the effectiveness of proposed stabilization scheme in treating non-physical 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 source-point. 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 Schur-complement 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 near-singular 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 near-singular modes cause a dramatic degradation in the linear solver performance.

6.1.3 Modified Undrained Barry-Mercer

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 rate-controlled to pressure-controlled. 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

Figure 13: Problem geometry for the multiphase poromechanics example. Grey region corresponds to a high-permeability channel, which spirals in a “staircase” fashion from the upper injection well (blue triangle) to the lower production well (red triangle). Blue region is a low-permeability zone.
Parameter Units Value
Porosity:
High-perm zone 0.20
Low-perm zone 0.05
Permeability:
High-perm zone mD 1000
Low-perm zone mD 1
Relative perm:
Residual wetting sat. 0.2
Residual non-wetting sat. 0.2
Wetting Fluid:
Reference density kg/m 1035
Bulk modulus MPa
Viscosity cP 0.3
Non-wetting 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
Time-stepping:
Initial day 0.0001
Maximum day 1
End time day 100
Solver tolerances:
Newton 10
Krylov 10
Table 3: Simulation parameters used in the multiphase example.
(a) Colorbar is pressure in MPa at day for the unstabilized (left) and stabilized (right) formulations. Note that the colorbar is truncated to visually accentuate oscillations.
(b) Colorbar is pressure in MPa at day for the unstabilized (left) and stabilized (right) formulations.
(c) Colorbar is saturation in decimals at day for the unstabilized (left) and stabilized (right) formulations.
Figure 14: Comparison of unstabilized and stabilized formulations for the multiphase example. The stabilization suppresses checkerboarding at early simulation times (a), but does not otherwise compromise solution quality at late times (b-c).
Figure 15: Comparison of average Krylov iterations per Newton step for the stabilized and unstabilized formulation applied to the multiphase poromechanics example. Without stabilization, the iterative solver convergence degrades at early simulation times when small time steps are employed.

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 high-permeability channel and a low-permeability host rock. The high-permeability 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 bottom-hole 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 high-perm 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 near-singular modes, to which Krylov-based 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 near-singular modes. The latter is a critical issue for large-scale 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 FC-MAELSTROM 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 DE-FG02-03ER15454. 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 DE-AC52-07-NA27344.

References

  • B. A. (1941) General theory of three-dimensional consolidation. Journal of Applied Physics 12 (2), pp. 155–164. External Links: Document Cited by: §1.
  • K. Aziz (1979) Petroleum reservoir simulation. Vol. 476, Applied Science Publishers, London. Cited by: §3.
  • I. Babuška (1973) The finite element method with Lagrangian multipliers. Numerische Mathematik 20 (3), pp. 179–192. External Links: ISSN 0945-3245, Document Cited by: §1.
  • W. Bangerth, R. Hartmann, and G. Kanschat (2007) deal.II – a General Purpose Object Oriented Finite Element Library. ACM Trans. Math. Softw. 33 (4), pp. 24/1–24/27. Cited by: §6.
  • S. Barry and G. Mercer (1999) Exact solutions for two-dimensional time-dependent 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.
  • L. Berger, R. Bordas, D. Kay, and S. Tavener (2015) Stabilized lowest-order finite element approximation for linear three-field poroelasticity. SIAM Journal on Scientific Computing 37 (5), pp. A2222–A2245. External Links: Document Cited by: §1.
  • P. Bochev, C. Dohrmann, and M. Gunzburger (2006) Stabilization of low-order mixed finite elements for the Stokes equations. SIAM Journal on Numerical Analysis 44 (1), pp. 82–101. External Links: Document Cited by: §1.
  • D. Boffi, M. Botti, and D. A. Di Pietro (2016) A nonconforming high-order 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.
  • J. M. Boland and R. A. Nicolaides (1983) 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.
  • R. I. Borja, X. Liu, and J. A. White (2012) Multiphysics hillslope processes triggering landslides. Acta Geotechnica 7 (4), pp. 261–269. External Links: ISSN 1861-1133, Document Cited by: §1.
  • R. I. Borja and J. A. White (2010) Continuum deformation and stability analyses of a steep hillside slope under rainfall infiltration. Acta Geotechnica 5 (1), pp. 1–14. External Links: ISSN 1861-1133, Document Cited by: §1.
  • F. Brezzi and J. Pitkäranta (1984) On the stabilization of finite element approximations of the Stokes equations. In Efficient Solutions of Elliptic Systems: Proceedings of a GAMM-Seminar Kiel, January 27 to 29, 1984, W. Hackbusch (Ed.), pp. 11 – 19. External Links: ISBN 978-3-663-14169-3, Document Cited by: §1.
  • F. Brezzi (1974) On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers. R.A.I.R.O. Analyse Numérique 8 (R2), pp. 129–151. External Links: Document Cited by: §1.
  • E. Burman and P. Hansbo (2007) A unified stabilized method for Stokes and Darcy’s equations. Journal of Computational and Applied Mathematics 198 (1), pp. 35 – 51. External Links: ISSN 0377-0427, Document Cited by: §1.
  • J. Camargo, R. Q. Velloso, and E. A. Vargas (2016) Numerical limit analysis of three-dimensional slope stability problems in catchment areas. Acta Geotechnica 11 (6), pp. 1369–1383. External Links: ISSN 1861-1133, Document Cited by: §1.
  • J. Cao, J. Jung, X. Song, and B. Bate (2018) 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 1861-1133, Document Cited by: §1.
  • N. Castelletto, J. A. White, and H. A. Tchelepi (2015) Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics. International Journal for Numerical and Analytical Methods in Geomechanics 39 (14), pp. 1593–1618. External Links: Document Cited by: §1.
  • Z. Chen, G. Huan, and Y. Ma (2006) Computational Methods for Multiphase Flows in Porous. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. Cited by: §2.
  • J. Choo and R. I. Borja (2015) 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 0045-7825, Document Cited by: §1.
  • J. Choo and S. Lee (2018) 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.
  • O. Coussy (2005) Poromechanics. John Wiley & Sons, Ltd. External Links: Document Cited by: §1.
  • R. de Boer (2012) Theory of porous media: highlights in historical development and current state. Springer Science & Business Media. External Links: Document Cited by: §1.
  • C. R. Dohrmann and P. B. Bochev (2004) 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.
  • H. C. Elman, D. J. Silvester, and A. J. Wathen (2005) Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Oxford University Press. Cited by: §1, §5.
  • T. Ertekin, J. H. Abou-Kassem, and G. R. King (2001) Basic applied reservoir simulation. Society of Petroleum Engineers, Richardson. Cited by: §2.
  • A. H. Fávero Neto and R. I. Borja (2018) Continuum hydrodynamics of dry granular flows employing multiplicative elastoplasticity. Acta Geotechnica 13 (5), pp. 1027–1040. External Links: ISSN 1861-1133, Document Cited by: §1.
  • Y. Feng and K. E. Gray (2019) XFEM-based cohesive zone approach for modeling near-wellbore hydraulic fracture complexity. Acta Geotechnica 14 (2), pp. 377–402. External Links: ISSN 1861-1133, Document Cited by: §1.
  • M. Fortin and S. Boivin (1990) 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.
  • G. Fu (2019) A high-order HDG method for the Biot’s consolidation model. Computers & Mathematics with Applications 77 (1), pp. 237 – 252. External Links: ISSN 0898-1221, Document Cited by: §6.1.1, §6.1.2.
  • [30] O. Ghaffaripour, G. A. Esgandani, A. Khoshghalb, and B. Shahbodaghkhan Fully coupled elastoplastic hydro-mechanical analysis of unsaturated porous media using a meshfree method. International Journal for Numerical and Analytical Methods in Geomechanics. External Links: Document Cited by: §1.
  • M. Gholami Korzani, S. A. Galindo-Torres, A. Scheuermann, and D. J. Williams (2018) SPH approach for simulating hydro-mechanical processes with large deformations and variable permeabilities. Acta Geotechnica 13 (2), pp. 303–316. External Links: ISSN 1861-1133, Document Cited by: §1.
  • J. B. Haga, H. Osnes, and H. P. Langtangen (2012) On the causes of pressure oscillations in low-permeable and low-compressible porous media. International Journal for Numerical and Analytical Methods in Geomechanics 36 (12), pp. 1507–1522. External Links: Document Cited by: §6.1.3.
  • H. T. Honório, C. R. Maliska, M. Ferronato, and C. Janna (2018) A stabilized element-based finite volume method for poroelastic problems. Journal of Computational Physics 364, pp. 49 – 72. External Links: ISSN 0021-9991, Document Cited by: §1, §1.
  • T. J.R. Hughes, L. P. Franca, and M. Balestra (1986) A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations. Computer Methods in Applied Mechanics and Engineering 59 (1), pp. 85 – 99. External Links: ISSN 0045-7825, Document Cited by: §1.
  • T. J.R. Hughes and L. P. Franca (1987) A new finite element formulation for computational fluid dynamics: VII. the Stokes problem with various well-posed 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 0045-7825, Document Cited by: §1.
  • W. Jin and C. Arson (2019) Fluid-driven transition from damage to fracture in anisotropic porous media: a multi-scale XFEM approach. Acta Geotechnica. External Links: ISSN 1861-1133, Document Cited by: §1.
  • J. Kim, H. A. Tchelepi, and R. Juanes (2011) Stability, accuracy, and efficiency of sequential methods for coupled flow and geomechanics. SPE Journal 16 (02), pp. 249–262. External Links: Document, ISBN 1086-055X Cited by: §1.
  • M. Mashhadian, S. Abedi, and A. Noshadravan (2018) Probabilistic multiscale characterization and modeling of organic-rich shale poroelastic properties. Acta Geotechnica 13 (4), pp. 781–800. External Links: ISSN 1861-1133, Document Cited by: §1.
  • E. Mikaeili and B. Schrefler (2018) XFEM, strong discontinuities and second-order work in shear band modeling of saturated porous media. Acta Geotechnica 13 (6), pp. 1249–1264. External Links: ISSN 1861-1133, Document Cited by: §1.
  • P. Navas, L. Sanavia, S. López-Querol, and R. C. Yu (2018) Explicit meshfree solution for large deformation dynamic problems in saturated porous media. Acta Geotechnica 13 (2), pp. 227–242. External Links: ISSN 1861-1133, Document Cited by: §1.
  • J. M. Nordbotten (2014) Cell-centered 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.
  • J. M. Nordbotten (2016) Stable cell-centered finite volume discretization for Biot equations. SIAM Journal on Numerical Analysis 54 (2), pp. 942–968. External Links: Document Cited by: §1.
  • F. Oka, B. Shahbodagh, and S. Kimoto (2019) A computational model for dynamic strain localization in unsaturated elasto-viscoplastic soils. International Journal for Numerical and Analytical Methods in Geomechanics 43 (1), pp. 138–165. External Links: Document Cited by: §1.
  • D. W. Peaceman (1978) Interpretation of well-block pressures in numerical reservoir simulation. Soc. Pet. Eng., AIME J. 18 (3), pp. 183–194. External Links: Document Cited by: §2.
  • S. Peng, Z. Zhang, J. Mou, B. Zhao, Z. Liu, and J. Wang (2019) 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.
  • P. J. Phillips and M. F. Wheeler (2007) A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discrete-in-time case. Computational Geosciences 11 (2), pp. 145–158. External Links: ISSN 1573-1499, Document Cited by: §6.1.1.
  • J. Pitkäranta and T. Saarinen (1985) 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.
  • S. H. Prassetyo and M. Gutierrez (2018) High-order ADE scheme for solving the fluid diffusion equation in non-uniform grids and its application in coupled hydro-mechanical simulation. International Journal for Numerical and Analytical Methods in Geomechanics 42 (16), pp. 1976–2000. External Links: Document Cited by: §1.
  • M. Preisig and J. H. Prévost (2011) 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.
  • J. H. Prévost (2014) Two-way coupling in reservoir geomechanical models: vertex-centered Galerkin geomechanical model cell-centered and vertex-centered finite volume reservoir models. International Journal for Numerical Methods in Engineering 98 (8), pp. 612–624. External Links: Document Cited by: §1.
  • C. Rodrigo, F. Gaspar, X. Hu, and L. Zikatanov (2016) 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.
  • C. Rodrigo, X. Hu, P. Ohm, J. H. Adler, F. J. Gaspar, and L. Zikatanov (2018) New stabilized discretizations for poroelasticity and the Stokes’ equations. Computer Methods in Applied Mechanics and Engineering 341, pp. 467–484. Cited by: §1.
  • S. J. Semnani, J. A. White, and R. I. Borja (2016) 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. Settari and F. M. Mourits (1998) A coupled reservoir and geomechanical simulation system. SPE Journal 3 (03), pp. 219–226. External Links: Document Cited by: §1.
  • S. Shiozawa, S. Lee, and M. F. Wheeler (2019) The effect of stress boundary conditions on fluid-driven fracture propagation in porous media using a phase-field modeling approach. International Journal for Numerical and Analytical Methods in Geomechanics 43 (6), pp. 1316–1340. External Links: Document Cited by: §1.
  • D.J. Silvester and N. Kechkar (1990) Stabilised bilinear-constant velocity-pressure 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 0045-7825, Document Cited by: §1, §1, §5.
  • D. Silvester (1994) Optimal low order finite element methods for incompressible flow. Computer methods in applied mechanics and engineering 111 (3-4), pp. 357–368. External Links: ISSN 0045-7825, Document Cited by: §1, §5.
  • I. Sokolova, M. G. Bastisya, and H. Hajibeygi (2019) Multiscale finite volume method for finite-volume-based simulation of poroelasticity. Journal of Computational Physics 379, pp. 309 – 324. External Links: ISSN 0021-9991, Document Cited by: §1.
  • X. Song, K. Wang, and M. Ye (2018) Localized failure in unsaturated soils under non-isothermal conditions. Acta Geotechnica 13 (1), pp. 73–85. External Links: ISSN 1861-1133, Document Cited by: §1.
  • R. Stenberg (1984) 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.
  • W. Sun, J. T. Ostien, and A. G. Salinger (2013) 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.
  • P. Teatini, M. Ferronato, G. Gambolati, and M. Gonella (2006) Groundwater pumping and land subsidence in the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend. Water Resources Research 42 (1), pp. . External Links: Document Cited by: §1.
  • A. Truty and T. Zimmermann (2006) Stabilized mixed finite element formulations for materially nonlinear partially saturated two-phase media. Computer Methods in Applied Mechanics and Engineering 195 (13), pp. 1517 – 1546. External Links: ISSN 0045-7825, Document Cited by: §1.
  • J. P. Verdon, J. Kendall, A. L. Stork, R. A. Chadwick, D. J. White, and R. C. Bissell (2013) Comparison of geomechanical deformation induced by megatonne-scale CO2 storage at Sleipner, Weyburn, and In Salah. Proceedings of the National Academy of Sciences 110 (30), pp. E2762–E2771. External Links: Document, ISSN 0027-8424 Cited by: §1.
  • Wan,Jing (2002) Stabilized finite element methods for coupled geomechanics and multiphase flow. Ph.D. Thesis, Stanford University. Cited by: §1.
  • M. Wang, Y.T. Feng, G.N. Pande, and T.T. Zhao (2018) A coupled 3-dimensional bonded discrete element and lattice Boltzmann method for fluid-solid coupling in cohesive geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics 42 (12), pp. 1405–1424. External Links: Document Cited by: §1.
  • W. Wang, M. Zhou, B. Zhang, and C. Peng (2019) A dual mortar contact method for porous media and its application to clay-core rockfill dams. International Journal for Numerical and Analytical Methods in Geomechanics 43 (9), pp. 1744–1769. External Links: Document Cited by: §1.
  • J. A. White, N. Castelletto, S. Klevtsov, Q. M. Bui, D. Osei-Kuffuor, and H. A. Tchelepi (2018) A two-stage preconditioner for multiphase poromechanics in reservoir simulation. arXiv preprint arXiv:1812.05540. Cited by: §1, §3, §3, §6.2, §6.2.
  • J. A. White, N. Castelletto, and H. A. Tchelepi (2016) Block-partitioned solvers for coupled poromechanics: a unified framework. Computer Methods in Applied Mechanics and Engineering 303, pp. 55–74. External Links: Document Cited by: §1.
  • J. A. White, L. Chiaramonte, S. Ezzedine, W. Foxall, Y. Hao, A. Ramirez, and W. McNab (2014) 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.
  • J. A. White and R. I. Borja (2008) Stabilized low-order finite elements for coupled solid-deformation/fluid-diffusion and their application to fault zone transients. Computer Methods in Applied Mechanics and Engineering 197 (49), pp. 4353 – 4366. External Links: ISSN 0045-7825, Document Cited by: §1.
  • C. Yan, Y. Jiao, and S. Yang (2019) A 2D coupled hydro-thermal model for the combined finite-discrete element method. Acta Geotechnica 14 (2), pp. 403–416. External Links: ISSN 1861-1133, Document Cited by: §1.
  • C. Yan and Y. Jiao (2019) FDEM-TH3D: A three-dimensional 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.
  • Q. Zhang, J. Choo, and R. I. Borja (2019) On the preferential flow patterns induced by transverse isotropy and non-Darcy flow in double porosity media. Computer Methods in Applied Mechanics and Engineering 353, pp. 570 – 592. External Links: ISSN 0045-7825, Document Cited by: §1.
  • Y. Zhao and R. I. Borja (2019) Deformation and strength of transversely isotropic rocks. In Desiderata Geotechnica, W. Wu (Ed.), Cham, pp. 237–241. External Links: ISBN 978-3-030-14987-1, Document Cited by: §1.
  • Y. Zhao, S. J. Semnani, Q. Yin, and R. I. Borja (2018) 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.
  • Z. Zhou, X. Du, S. Wang, X. Cai, and L. Chen (2019) Micromechanism of the diffusion of cement-based 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 1861-1133, Document Cited by: §1.
  • M. D. Zoback (2007) Reservoir geomechanics. Cambridge University Press. External Links: Document Cited by: §1.