Finite volume hydromechanical simulation in porous media.

Cell-centered finite volume methods are prevailing in numerical simulation of flow in porous media. However, due to the lack of cell-centered finite volume methods for mechanics, coupled flow and deformation is usually treated either by coupled finite-volume-finite element discretizations, or within a finite element setting. The former approach is unfavorable as it introduces two separate grid structures, while the latter approach loses the advantages of finite volume methods for the flow equation. Recently, we proposed a cell-centered finite volume method for elasticity. Herein, we explore the applicability of this novel method to provide a compatible finite volume discretization for coupled hydromechanic flows in porous media. We detail in particular the issue of coupling terms, and show how this is naturally handled. Furthermore, we observe how the cell-centered finite volume framework naturally allows for modeling fractured and fracturing porous media through internal boundary conditions. We support the discussion with a set of numerical examples: the convergence properties of the coupled scheme are first investigated; second, we illustrate the practical applicability of the method both for fractured and heterogeneous media.

Formulation of consistent discretization of fluid flow and mechanical deformation within a finite volume framework

Numerical validation of convergence and application to model problems

Application to fractured and fracturing porous media

1 Introduction

Geomechanical response to subsurface exploration has long been accepted as important as evidenced through subsidence associated with both petroleum and ground water extraction. Coupled hydromechanical processes are furthermore emerging as increasingly important in a range of rapidly expanding subsurface engineering applications. These include, but are not limited to, CO2 storage, unconventional oil and gas production, and geothermal energy recovery. Particular cases have gained recent attention: surface uplift due CO2 injection has emerged as an appealing monitoring mechanism [Rutquist et al., 2010]; the shale-gas revolution is enabled due to improved understanding of hydraulic fracturing [Cipolla et al., 2013]; and finally unlocking geothermal resources heavily relies on accurate understanding of fluid flow, heat transfer, and mechanical response [Tester, 2006].

Common for most subsurface applications is that the fluids are of primary economic or environmental interest. Thus, the simplest analysis, in terms of screening, research, and decision purposes, involves only flow calculations [see, e.g., Celia and Nordbotten, 2009]. As additional physical complexity is required, we therefore frequently confront a situation wherein the mechanical response is desired in a context for which the flow calculations are already established. This motivates us to seek ways to seamlessly integrated mechanical response within the context of flow in porous media. We note that this is in contrast to problems related to earthquakes, wave loading, and consolidation, where it is the material deformation which is of primary interest, and the fluid flow is of secondary importance [Lewis and Schrefler, 1998].

Two strategies can be readily identified in order to couple mechanical responses and flow simulations. The classical approach is to consider a separate framework for mechanical deformations, and couple it to the flow calculations by passing appropriate variables. This can be implemented in both a sequential or iterative setting. There has been recent work providing rigorous analysis of coupling schemes for hydromechanical systems [see, e.g., Kim et al., 2011; Mikelic and Wheeler, 2012]. This approach is also standard for industrial simulation, as evidenced by the software packages Visage and Eclipse, respectively, for mechanics and flow, both marketed by Schlumberger. This classical approach has the advantage that established numerical methods and software packages can be used for the two subproblems. However, the coupling of different software packages and grid structures may lead to both excessive numerical diffusion, as well as slow convergence of the coupling schemes [Pettersen, 2012].

Herein, we advocate a different approach. Our perspective is that the mechanical response needs to be within the same numerical and software frameworks as the flow calculations. This will allow equivalent grids and parameter representations to be used in the two parts of the problem, moreover it avoids the use of iterative coupling schemes, thus facilitating simultaneous and fully coupled solution strategies [Haga, 2011]. Fully coupled simulation has previously been recognized to have superior accuracy [see, e.g., Lewis and Schrefler, 1998]; however, it has to the best of our knowledge primarily been implemented with node-centered or face-centered variables for the mechanical deformation (e.g., finite elements [Lewis and Schrefler, 1998; Kim et al., 2011; Haga, 2011] or face-centered finite volume [Lemaire, 2013]). An example of node-centered variables for deformation and cell-centered variables for flow is also analyzed in Wheeler et al. [2013]. Since the majority of flow and transport simulation in porous media relies on a cell-centered finite volume structure, we advocate that coupled hydromechanical simulations must be resolved within a consistent framework utilizing cell-centered finite volumes also for the deformation problem. Utilizing finite volume methods also for elasticity gives the practical advantage of equivalent data structures for flow and deformation, which simplifies both the design of linear and nonlinear solvers, as well as the inclusion of irregular domain features such as fractures (as we will discuss later). A similar idea was briefly considered in the petroleum literature [Shaw and Stone, 2005; Stone et al., 2000], using finite volume methods similar to those proposed in Jasak and Weller [2000]. Their approach is however not suited to accurately handle discontinuities [Liu et al., 2004].

A cell-centered discretization for mechanics was recently introduced by the author for the purpose of geomechanical applications (J. M. Nordbotten, 2013, Cell-centered finite volume discretizations for deformable porous media, submitted to International Journal of Numerical Methods in Engineering). Herein, we explore the coupling of this novel discretization with a compatible discretization for flow. This involves in particular identifying appropriate and consistent treatment of the coupling terms between mechanics and flow. A particular feature of cell-centered variables is furthermore that they represent the mean solution within the grid cell, and as such do not explicitly provide a spatially continuous approximation. While this may be a deficiency in some applications, we show that this allows for a particularly simple approach to modeling fractured materials.

As a comment, we note that it is not appropriate in the context of fluid flow and deformation in heterogeneous porous media to consider finite volume methods as simply an interpretation of finite element methods. Indeed, the notion of conservative postprocessing of finite elements [Hughes et al., 2000] does not align with the discontinuous material coefficients prevailing in porous media, and arithmetic means are recovered where harmonic means are desired [Babuska and Osborn, 1983]. The treatment of fractures discussed below would also not adapt naturally to such an approach. On the other hand, deriving finite volume methods by exploiting reduced integration for mixed finite element methods [Russell and Wheeler, 1983] is complicated by the lack of simple mixed finite element methods for elasticity [Arnold and Winther, 2002].

We structure the manuscript as follows. In the next section, we state the model equations, including the description of fractures as internal boundary conditions. This is followed in section 3 by the identification of appropriate coupling between flow and deformation in a cell-centered finite volume framework. We validate our approach in section 4 by comparison to analytical solutions. These validations are supplemented in section 5 by examples containing prototypical problems from applications. Finally, in section 6, we summarize and conclude the paper.

2 Governing Equations

We will model the porous medium as a poroelastic material according to Biot's equations [Biot, 1941]. Furthermore, we emphasize the inclusion of fracture faces in the form of internal boundaries of the material. In general, fluid flow equations can be considered within the fractures, we will herein, for simplicity, make the approximation that the potential drop within a (connected) fracture system is low, such that the fluid potential is constant within any connected fracture network. This conceptualization leads to the model equations as given below. For a more complete description, see, e.g., Lewis and Schrefler [1998].

2.1 Equations for the Porous Material

For concreteness, we will denote our model domain as Ω. For this, and other domains, we denote the boundary using the symbol ∂, such that the boundary of our model domain is given as ∂Ω.

Acceleration is given by the imbalance of mechanical forces according to Newton's law for any subdomain ω of our model domain which is fixed in time [Temam and Miramville, 2000],

Here we denote the material deformation as , the force exerted on a surface with normal vector as , the density of the material as ρ, and finally the forces acting on the bulk are given by the gravitational vector .

Similarly, the accumulation of fluid mass for any subdomain is given by the imbalance of fluid mass flowing through the boundary and the internal sources or sinks,

Here fluid mass is denoted m, the fluid flux through a surface with normal vector as , the density of the fluid ρ, and finally sources (or sinks) are denoted by q.

For both the fluid flux and the surface forces, a geometric argument (Cauchy's first law) allows us to write the dependency of the surface terms on the normal vector as a linear function, such that we have both

The variables and are known as the Cauchy stress tensor and the Darcy flux, respectively.

Equations (1) and (2) represent the fundamental balance laws from which the finite volume methods are derived. Both for the analysis of the equations, as well as to derive the finite element discretizations, it is common to consider the differential representation of the system. For balance of momentum, the differential equation is while for conservation of mass

In the continuation, we will work with equations (1-3) directly, without explicitly considering the differential equations (4) and (5).

The conservation laws given above are complemented by constitutive laws. Here we consider only linear constitutive laws, corresponding to small deformations of an elastic porous material with slow flow of the fluid. Then for the solid, Biot introduced the fluid pressure as an addition to Hooke's law, to obtain that the stress is linearly dependent on both the displacement gradient and the fluid pressure p

The fourth-order tensor is the stiffness tensor of the material, and α is Biot's coupling coefficient.

For the fluid, the constitutive law is given by Darcy as:

We will abuse convention, and refer to the linear coefficient in Darcy's law as the permeability, suppressing the dependence on the fluid viscosity. Finally, we note that in a porous medium the fluid mass is proportional to the density and porosity, where the former is modeled as a function of pressure and the latter as a function of the compression of the solid,

Equations (1-3) together with equations (6-8) form a complete system of equations for flow and deformation of porous media. Our exposition has neglected the nuances which appear between Eularian and Lagrangian coordinates, and it is therefore implicit that this formulation is only valid in the limit of small deformations. For increasing deformations, the transformation of variables induced by the deformation itself becomes of importance; this is beyond the scope of the current paper.

2.2 Boundary Conditions and Inclusion of Fractures

Fractures physically represent internal discontinuities in the porous material, and may be subject to other governing equations than those governing the flow and deformation within the material. Here we will only consider two types of fractures: closed fractures and open fractures. These will both be included as internal boundary conditions within the domain, but are distinguished by different boundary conditions.

To be precise, we can without loss of generality assume that the domain contains a single connected fracture system (multiple fractures systems are handled equivalently). A fracture system may contain multiple fractures, but by definition each fracture system is connected. Let the external boundary of the domain Ωbe denoted Γ, and the boundary between the fracture and the domain be denoted γ. For an open fracture, the boundary γ will encompass a volume, while for a closed fracture γ will be a simple surface.

To solve the equations given in section 2.1, we require boundary conditions on the full boundary of the domain. For the outer boundary Γ, we assume that suitable boundary conditions have been chosen. Out interest is the boundary conditions for the fracture system γ.

2.2.1 Closed Fractures

We define closed fractures as closed surfaces γ, which are either permeable or impermeable to flow, but are in full contact from a mechanical perspective. Thus, denoting by superscripts + and − the variables on the two sides of the internal boundary, with flux continuity (Neumann) boundary conditions for the fluid:

Where for the impermeable fractures we have, while for permeable fractures, we have the additional constraint of continuity of pressure (although note that for slightly permeable fractures there will in general be a jump in pressure proportional to the flux w). We complement the fluid conditions by continuity in displacement and surface forces for the solid:

Considering the closed fractures in the context of the continuity equations (1) and (2), we note that permeable closed fractures are themselves simply expressions of conservations (and thus do not affect the system), while impermeable closed fractures represent an impermeable internal surface to the fluid flow. Our main interest will therefore be in open fractures.

2.2.2 Open Fractures

We define open fractures as surfaces γ enclosing a fluid volume. We will assume for simplicity that the variation in fluid potential within this volume is negligible. This assumption is valid in the limit of fluids or gases of low viscosity—the extension to more complex models of flow in the fracture is possible [see, e.g., Sousa et al., 1993]; however, since the focus of this contribution is on the hydromechanical coupling in the porous medium, we will prefer to keep the simplest model for fluid flow in the fractures. Denote thus the constant fluid potential within the fracture as ψ. The fluid flow equation then assumes a Dirichlet boundary condition on γ

Furthermore, the solid stress is carried purely by the fluid pressure at the boundary

In applications where a pressure-controlled well injects directly into a fracture network (which will sometimes be the case for hydraulic fracturing), the fracture potential ψ will be known. Otherwise, the fracture itself will satisfy the conservation law (4), which provides a closure relationship for ψ.

3 Finite Volume Approximation of Coupling Terms

We seek to develop a finite volume approximation for the equations of hydromechanics as outlined in section 1. We will achieve this goal by exploiting standard finite volume discretizations for the mass balance equations, coupled with a recent finite volume discretization for the deformation. The key issue in this section will be to properly understand the discrete representation of the coupling terms.

3.1 Finite Volume Setting

To be precise, let ω represent the cells of a nonoverlapping partition of the domain. Then equations (1) and (2) must hold for very cell, and we obtain the balance of mechanical forces for each cell i as

Here we have used Biot's stress tensor, and split the surface integral to emphasize the structure. The left-hand side of (13) is thus the standard terms appearing with the stress tensor from Hook's law, while the surface integral on the right-hand side of the equation represents the coupling between the mechanical deformation and the fluid.

Before we state the conservation of mass, we note that from equation (8) we have as a consequence of the chain rule

Here is the compressibility of the fluid, while the porosity is assumed to only be a function of the trace of the strain, hence. For the conservation of fluid mass, we then obtain

In the quasi-linear regime (or within an iterative solver for the fully nonlinear system), this simplifies to the model equation for mass balance

Again we have reorganized the equation to emphasize the structure. The left-hand side terms of equation (16) are the standard terms appearing in discretizing flow in porous media. The volume integral on the right-hand side represents the coupling with the mechanical deformation arising from the change in volume available for the fluid.

Finite volume approximations to equation (13) and (16) are usually written in terms of cell-average and face-average quantities. We define p, , and q as the cell-average values of p, , and q, respectively. For the face separating cells i and j, we further define the face-averaged fluid fluxes and mechanical stresses as

Note carefully that only represents the mechanical component of the face-average force.

With these definitions, the model equations (13) and (16) can be stated as and

It is important to note that equations (18) and (19) are exact expressions of conservation.

Local finite volume approximations to the fluxes are well known in literature. The inconsistent two-point flux (using only the values of p in the cells neighboring the face) is most widely used [Aziz and Settari, 1979], while consistent approximations have received significant attention the last two decades [see, e.g., Aavatsmark, 2002; Aavatsmark et al., 1996; Klausen and Winther, 2006]. Herein we chose to apply the multipoint flux approximations (MPFA) as described in Aavatsmark [2002]. In contrast, local finite volume approximations to normal stress have not received similar attention, since finite element approximations or node-centered finite volume methods have generally been preferred for deformations [Lewis and Schrefler, 1998]. As noted in the introduction, the MPFA methods have recently been extended to normal stress, where they are referred to as multipoint stress approximations (MPSA; J. M. Nordbotten, submitted manuscript, 2013).

We will adopt the MPFA and MPSA stencils for the internal fluxes and surface forces appearing in equations (18) and (19). This approach is also applicable for Dirichlet boundary conditions on both internal and external boundaries. For the fluxes and surface forces which align with internal or external Neumann boundaries, we will explicitly use the boundary conditions directly in equations (18) and (19). Note that the choice of discrete flux and surface force representations inherit the limitations associated with these methods. For the MPFA methods, limitations are well understood, with the most general argument established by Agelas and coworkers [Agelas et al., 2010]. For the MPSA method, similar analysis is not yet complete due to the added challenges associated with the nontrivial kernel for the elasticity operator.

It remains to obtain the coupling terms on the right-hand side of these equations. Fortunately, both pressures and displacements are available at cell faces as part of the MPFA and MPSA calculations. In both instances, these expressions have received little attention, since the MPFA and MPSA stencils have not been considered for coupled hydromechanics before. We will therefore review them below.

3.2 Approximation of Coupling Terms

We will approximate the coupling terms in the framework of the MPFA and MPSA methods. To this aim, let us first briefly review the main idea behind these methods, while referring to other publications for details (J. M. Nordbotten, submitted manuscript) [Aavatsmark, 2002]. We will in the following few paragraphs for simplicity only refer to the scalar discretization (MPFA), as the situation for the vector equation is analogous.

Considering a grid as illustrated by solid black lines in Figure 1, where cell centers are denoted x. The MPFA method then utilizes a concept of a dual grid, whose cells are termed “interaction regions,” as illustrated by gray solid lines. Note that by considering the intersection of the primal and dual grid each cell will be subdivided into multiple subcells, while each face will be subdivided into two (in 2-D) or multiple (in 3-D) subfaces. The key purpose of the MPFA methods is now to obtain explicit expressions of the fluxes for each subface.

Figure 1

Sketch of key concepts for MPFA and MPSA discretizations.

For concreteness consider face, which is in Figure 1 further subdivided into two subfaces. To calculate the flux (normal stress) across the upper subface, the MPFA O(η) methods utilize the following local calculation. Within each subcell, the pressure (displacement) is allowed to vary linearly. In 2-D, this gives rise to 3 degrees of freedom in each subcell. These degrees of freedom are constrained by matching the pressure in the cell center , and at one point on each subface, denoted. Furthermore, using Darcy's law, also flux is forced to be continuous across each subface. For the scalar equation, this leads to a solvable local system on most geometries [Aavatsmark, 2002], while for the vector equation, additional constraints are imposed (J. M. Nordbotten, submitted manuscript). We comment that the variable η refers to the position of the continuity points for pressure (displacement), and can in general be different from cell to cell [see, e.g., Nordbotten et al., 2007] or between the scalar and vector discretizations. The choice of η can in special cases lead to a higher-order methods [Edwards and Rogers, 1998], and more generally impacts the monotonicity of the discretization [Nordbotten et al., 2007].

Solving the local systems as outlined above give explicit expressions for the subface fluxes (normal stresses) as functions of the nearby cell-center pressures (displacements). This is the traditional use of the MPFA methods, which then determines the linear systems arising on the left-hand side of equations (18) and (19).

The above local calculation yields as an intermediate calculation the local expressions for the subface pressures at the continuity points. Note that since Darcy's law was used in the local calculation, this local expression for the subface pressure is not a simple geometric average of the nearby cell-center values, but also honors the effect of anisotropy and local heterogeneity in permeability. Using these subface pressures, the surface integral appearing in the right-hand side of equation (18) can readily be approximated with the continuity points taking the role of quadrature points for the integral.

Similarly, for the vector case, the displacement at continuity points is a consequence of the local calculations. As in the scalar case, the resulting expressions are not simply dependent on the geometry of the grid, but also honor the physics as expected from Hooke's law. As in the scalar case, the surface integral appearing in the right-hand side of equation (19) can readily be approximated.

The fact that the calculation of the coupling terms are dependent on the underlying constitutive laws, implies that even though the coupling terms are formally adjoint operators for the continuous problem, their matrix representations in the discrete setting will in general not be negative transposes. This forms a contrast to pure finite element discretizations for hydromechanics, but is similar to the situation for some of the coupled schemes as described in the introduction.

3.3 Application to Fracturing Materials

As discussed in section 2, we consider fractures as boundaries of the poroelastic domain. The finite volume methods of the preceding section naturally accept both Dirichlet and Neumann boundary conditions, and as such static fractures pose no particular problems [Aavatsmark, 2002].

Of interest are also materials where the fractures evolve over time. We consider the setting where the fracture evolution is governed by a prescribed fracture mechanical model at a local scale [see, e.g., Boone et al., 1986], and are concerned with how this fracture evolution can be combined with the overall poroelastic discretization. In this setting, there are two main classes of models. In applications where the precise structure of the fracture is important, significant (adaptive) grid refinement will be desired around the fracture tip. As the grid is locally updated for each time step, this is from a spatial discretization perspective equivalent to the problem with static fractures. The finite volume methods utilized herein are well adapted to grid refinement as required in this class of models [Aavatsmark et al., 2010]. However, the novelty of this contribution does not lie in adaptive grid refinement.

The second class of models concerns applications where the fractures are part of a much larger domain of interest, the grid is typically not refined around the individual fractures, and thus the grid is conceptually static despite the changing fractures. This is the case where the present methodology has particular advantages.

For standard methods for solid mechanics, such as finite elements, the displacement is by choice of function spaces inherently continuous. To accommodate fractures, either nodes must be “doubled,” or special elements are introduced, with fractures propagating internally to the elements of the grid. One such concept is the extended finite elements [Sukumar et al., 2000]. While this is conceptually straight-forward, practical problems arise when multiple fractures are present within the same element, especially for general grids in 3-D. Furthermore, in both the above approaches the number of degrees of freedom increases during the simulation, leading to additional work in updating the data structures for the linear solvers in time-dependent problems [Sukumar et al., 2000]. While extended finite elements have been extended to fluid flow in fractures [e.g., Fumagalli, 2012], this author is not aware that full hydromechanical simulation has been undertaken with this methodology.

In the finite volume setting, fractures naturally coincide with the cell faces, since no explicit continuity is required by the discretization. This allows us to consider two possible models for time-evolving fractures. The simplest approach is to constrain fractures to follow existing faces of the original discretization. This has the advantage that the number of unknowns in the discretization remains unchanged during the simulation, and only the discretization where the fracture has moved needs to be updated. The second approach is to allow fractures to intersect a grid cell, effectively splitting it into two parts. This approach is conceptually similar to the extended finite elements, and thus provides limited benefit in terms of using a finite volume discretization.

In the results section, we will therefore emphasize with numerical examples simulations wherein the fractures that are allowed to evolve along the existing faces of the grid, allowing the grid structure to remain unchanged for the full duration of the simulation.

3.4 Prototypical Discretization

We close the section with recalling prototypical discretization of poroelasticity, based on the spatial stencils discussed in sections 3.1 and 3.2 [see also Lewis and Schrefler, 1998]. We refer to the standard finite volume discretization matrices as and for stress and fluid flux, respectively. Similarly, we refer to the coupling terms as grad and div for the coupling terms on the right-hand sides of equations (18) and (19), recognizing that the continuum limits these terms, represent the gradient and divergence, respectively. We consider a single-step discretization (so-called θ-Euler) of the temporal term, as is standard in industrial applications, after which the linear system can be expressed as

Here we have used, and represented all right-hand side entries by r, which for the purpose of time integration is defined as

The superscripts on the solution variable indicate the time level, and the discretization parameter θ is related to the time discretization. The customary choices are, which yield the so-called explicit, midpoint, and implicit discretizations, respectively.

In the general setting of heterogeneous parameters, the general structure of the system will remain the same. As such we will for simplicity frequently refer to equation (20), while keeping in mind the more general system given by equations (18) and (19).

Finally, we point out again that when arising from a finite element discretization (20) will be symmetric whenever the stress tensor is symmetric, as both and are symmetric and also. As noted previously in the section, none of these hold for finite volume discretizations, which lead to a (weakly) nonsymmetric system (20). This has three immediate implications. First, the stability of the discretization with respect to spatial grid refinement cannot be immediately established, and we will therefore verify the convergence of the coupled discretization numerically in the next section. Second, the linear solvers chosen for (20) must be designed for nonsymmetric systems (e.g., Generalized Minimal Residual (GMRES)). Finally, the time discretization may for certain problems have a lower bound on the time step in order to obtain stability [Lewis and Schrefler, 1998]. Unconditionally stable finite volume discretizations for poroelasticity is a topic of ongoing research [see, e.g., Lemaire, 2013]; at present, no results are known for cell-centered methods on general grids.

4 Numerical Validation

Finite volume discretizations in general are well established and validated. Of particular relevance for the current manuscript are theoretical results proving, in various settings, the convergence of the finite volume discretizations for the fluid problem [see, e.g., Klausen and Winther, 2006]. Both the fluid and mechanical discretizations have also been validated for a wide range of heterogeneities and grid distortions in, e.g., Aavatsmark et al. [2010], Eigestad and Klausen [2005], and J. M. Nordbotten (submitted manuscript). The purpose of this section is therefore restricted to novel issues, which may arise from the coupling. In particular, we will verify that (a) the discretization of the coupling terms found on the right-hand sides of equations (18) and (19) are convergent, and (b) that the fully coupled system on the form given in equation (20) is stable and convergent. Since our primary interest is in the spatial discretization, we will in both cases (a) and (b) restrict ourselves to considering refinements of the spatial grids. Here and in the examples in section 5, we specify the MPFA method as the so-called O(0) method, and the MPSA method as the generalized O-method.

Recognizing the importance of unstructured grids in order to capture the complex geometries associated with geological porous media, all convergence and stability tests are conducted with both structured (quadrilateral) and unstructured grids (triangles and general polygons).

We utilize the unit square where suitable analytical solutions can be defined by use of elementary differentiable functions. In particular, we will for these tests use the analytical solution and. These functions satisfy zero Dirichlet boundary conditions, and the problem is thus driven by internal source-sink terms (for the flow equation) and body forces (for the momentum equation).

We consider three types of grids, square, triangular, and general polygons. The square grid forms a regular lattice. The triangular grid is formed by subdividing each square of a perturbed regular lattice along the shortest diagonal. Finally, we obtain the unstructured grid by taking the dual grid of the triangular grid. Examples of triangular and unstructured grids at the third level of grid refinement are given in Figure 2. For each grid, we assign the characteristic grid density as the square root of the number of grid cells. This generalizes the notion of number of grid cells in each dimension. To clearly identify the behavior of the method, we also considered the case where the coupling terms are omitted (by setting α = 0 in equation (20)); however, the results were indistinguishable from the coupled situation and are therefore not displayed in separate figures.

Figure 2

Example of (left) triangular and (right) unstructured grids used in convergence study. Note that the unstructured grid is the dual of the triangular in the sense that cell centers and cell corners are interchanged. Furthermore, each edge in the triangular grid crosses exactly one edge of the unstructured grid.

The convergence results for the coupled problem are shown in Figures 3. Here two layers of information are provided. The black lines give the approximation error for the two problems, and show that the discretizations are convergent in both the primary variables, as well as the coupling terms. We note that for square grids, the approximation error for all terms converge with second-order accuracy. For triangular and unstructured grids, the pressures and displacements still converge as second order; however, the coupling terms are only first-order convergent. These results are comparable to what could be expected from lowest-order finite elements, where for the uncoupled problem the discretization is typically second-order accurate in primary variables but only first-order accurate in the coupling terms. The second layer of information is given with gray lines on the figure, where the coupling terms are applied to the pressures and displacements of the analytical solution. We see that only minor differences are observed, indicating that the observed accuracy of the coupling terms is inherent to their definition, and not a result of pollution from the accuracy of the discrete solutions.

Figure 3

Convergence of discretization errors on (top left) square, (top right) triangular, and (bottom left) unstructured grids. The figures show with black lines the normalized error in the discrete solution and the coupling terms as compared to an analytical solution on the unit square. For comparison, the error in the coupling terms assessed directly by application to the analytical solution is shown in gray.

5 Example Applications

We consider three example applications to assess the flexibility of the coupled finite volume discretizations presented herein. The first application is to a heterogeneous material subject to external loading, where we assess the time evolution of the compression of the poroelastic material. The second application is to a fractured porous media, where two fractures intersect at a nonorthogonal angle in the interior of the domain. Here we consider the steady state of the system. Finally, we apply the coupled discretization to the tensile fracturing of a heterogeneous rock. Throughout the three applications we utilize a range of grids, giving us further opportunity to assess the flexibility and robustness of the finite volume framework.

5.1 Loading of a Layered Material

This example has a geometry that is motivated by loading of a layered material at the laboratory scale.

We consider a 1 m by 1 m sample, wherein softer, more permeable, material is embedded into a harder, less permeable material, as illustrated in Figure 4. We are interested in the response of the material to a loading on the top left-hand half of the sample. To be specific, we take the two materials separated by 2 orders of magnitudes in their parameters. Thus, we consider the materials to have permeabilities of 10−12 m2 and 10−14 m2, respectively, and Lamé coefficients of 1 and 100 GPa. This is within the range of values that may be observed for rocks in the Earth's crust [Ji et al., 2010]. Furthermore, we take the width of the soft material to be 20 cm, and the external loading to be 106 N. For this case, we consider mixed boundary conditions wherein the lower boundary is fixed, while the remaining boundaries have zero stress. The boundaries are furthermore open to fluid flow. We have chosen a quadrilateral grid, where the corners are perturbed randomly only in the horizontal direction to not interfere with the layered heterogeneity.

Figure 4

Illustration of the geometry for the loading and fracture problems for sections (left) 5.1 and (right) 5.3. Dark gray indicates the harder material, arrow indicates region of loading, while the heavy straight line at the base of the figure indicates the side of the material which is clamped. See text for details.

Since there is no fluid forcing function, the steady state solution for this problem is equivalent to the elastic problem, and thus we are interested in the time-dependent solution. The solution after 10 s is shown in Figure 5 (left). Note that within the high-permeable region, the fluid has equilibrated to a near-constant pressure, while significant pressure buildup remains in the low-permeable regions. The solution also conforms to our expectation that the majority of the compaction happens in the soft layer, while the hard layers are either at rest (the lower layer) or exposed to near-rigid motion (upper layer). For this example, we also show the maximum pressure and the deformation as a function of time in Figure 5 (right). The figure clearly illustrates both how the heterogeneous material introduces two time scales into the problem, and also how the fluid sustains loading during the transient phase.

Figure 5

Solution to loading problem. (left) The solution after 10 s, with pressure as gray scale (5.5 MPa maximum variation) and red arrows deformation (0.018 cm maximum deformation). The deformation arrows are given for only a subset of the grid cells. (right) The time evolution of the maximum deformation (top-left corner) and maximum pressure buildup are shown. Note that units are chosen to make both graphs fit on the same axis.

This example thus verifies the applicability of the coupled finite volume approach to standard problems of fluid flow and loading with parameters typical of those encountered in subsurface applications.

5.2 Injection Into Fractured Material

This example has a geometry that is motivated by geothermal energy production. In this application, injection typically occurs into a fracture network, and it is therefore important to understand the response of the material due to such operations.

For concreteness, we consider a 1 km by 1 km domain, wherein two fractures are placed horizontally and at a 62° angle (two-to-one slope of the inclined fracture). We resolve the domain using a grid that conforms to the fractures, and choose triangular grid cells for this example. Notice that the grid is constructed so that the aspect ratio of the triangles degenerates toward two of the corners of the grid. This apparently presents no numerical problems, which is consistent with previous studies [Klausen et al., 2008], (J. M. Nordbotten, submitted manuscript). Zero Neumann and Dirichlet boundary conditions are used on the external flow and deformation boundaries, respectively.

As described in section 2.2.1, the fractures are assumed to be highly conductive relative to the solid, and the pressure drop within the fracture is neglected. The solid is taken with a permeability of 10−12 m2, and with Lamé parameters of 1 GPa. A steady state solution is then calculated based on the injection of 10−3 m3/s of fluid into the fracture system, which is withdrawn from a well placed at the location (250 m).

The resulting solution is shown in Figure 6, which also illustrates the domain geometry and grid. Here the fractures are shown in solid, bold lines. The pressure profile is illustrated in gray scale, while arrows indicate the deformation of the solid. As expected, a smooth pressure field drives flow from the fracture to the production well. The observed pressure variation between the fracture and the well is 6.4 MPa for this calculation. Simultaneously, the medium deforms in response to the injection. On the far side of the pumping well, relatively little deformation is seen, while on the side of the pumping well the medium deforms toward the well, with a concomitant opening of the fracture network. The maximum deformation observed in this example is 13 cm.

Figure 6

Simulation on fractured domain. Note the two intersecting fractures in bold. Gray scale indicates pressure (6.4 MPa variation), while the red arrows indicated deformation. The length of the arrows is proportional to the deformation, with the largest deformation 13 cm. To increase clarity, the deformation arrows are given for only a subset of the grid cells.

This example justifies the claim that the finite volume methodology adapts naturally to fractured porous media, and emphasizes that no particular treatment is needed for intersecting fractures.

5.3 Fracturing of a Heterogeneous Material

This final example considers the evolution of a hydraulically stimulated fracture within a heterogeneous rock. To stay within the framework of quasi-static poroelastic deformation, we will also assume that the fracture evolves in a quasi-static sense [Davy et al., 2013]. Thus, we let the fluid pressure and mechanical response equilibrate between each evolution step of the fracture. This is numerically advantageous, since the explicit treatment of the fracture would otherwise impose a Courant-Friedrichs-Lewy (CFL)-type condition associated with the speed of the fracture propagation.

The model geometry is given in Figure 5 (right), and the material properties are identical to example in section 5.1. However, in this case, there is no external mechanical loading. Instead the pressure boundary conditions are linear in the y coordinate, varying by 1 MPa from top to bottom of the sample. The fracture is seeded in the middle of the top boundary, and the pressure in the evolving fracture is set equal the pressure at the top boundary, such that fluid enters the rock from the fracture.

The fracture evolves along the grid edges, and we evolve the fracture one edge at a time by converting the internal edge with largest tensile force into a new fracture edge. This simplistic fracture evolution is meant to illustrate the flexibility of the current discretization, more advanced fracture evolution models are possible [Boone et al., 1986]. As the fracture evolves, no new unknowns are needed, and no new connectivity appears in the discrete equations—thus the data structure remains identical. Only the local discretization at the fracture tip needs to be updated to account for the appearance of a new internal boundary. This has significant advantages in terms of implementation and computation.

Due to the fractures being restricted to edges of the grid, the computation is prone to grid orientation effects. To investigate this issue, we use two random grids as benchmarks (one unstructured and one a Delaunay triangulation—see Figure 2), and compare the results to that obtained using a square lattice grid. The comparison is presented in Figure 7, using between 30,000 and 45,000 degrees of freedom in the calculation. We see that the calculations on both triangular and unstructured grids predict that the fracture deviates from a straight fracture as it approaches the heterogeneity, and after intersecting the relatively harder material deviates toward the edge of the sample. As the setup is symmetric, it is arbitrary what side of the sample the crack deviates toward, and the runs have been mirrored so that the crack always deviates toward the right boundary to facilitate comparison. The calculation on a uniform grid supports this general description of the fracture evolution, although significantly stronger grid orientation effects are seen for this grid. The calculations have been verified with several realizations of the grid, resulting in qualitatively similar behavior. Note that with both the uniform grid as well as the unstructured grid, minor secondary fractures are seen deviating from the main fracture path.

Figure 7

Solution to the fracturing problem. The problem is represented on (left) randomized triangular, (middle) randomized unstructured, and (right) uniform Cartesian grids, respectively. The grids have comparable diameters of the grid cells, with the two random grids have roughly 15,000 grid cells, while the uniform grid has 10,000 grid cells.

5.4 Comments on Implementation and Computational Issues

The numerical examples have all been computed using a Matlab implementation of the equations given above. As the governing partial differential equations (1-8) are linear, the only nonlinearity in the problem enters for the case of an evolving fracture as discussed in section 5.3.

In order to keep the emphasis on the spatial discretization, we have kept the temporal integration as simple as possible. To this end, we have employed a standard implicit (backward Euler) time discretization for the time derivative in the mass conservation equation for section 5.1. Section 5.2 is a steady state calculation. Finally, as noted in section 5.3, the fracture evolves in a quasi-static sense. This is realized by calculating the equilibrium pressure and displacement field between each update of the fracture tip. It is not realistic to evaluate the fracture evolution and the equilibrium fields in a coupled sense, since these are coupled through the (changing) discretization itself. Thus for all three sections, a fully coupled hydromechanical system needs to be solved for each time step. For the smaller examples, this system can be solved using a direct solver. For larger examples, standard iterative approaches are efficient [Mikelic and Wheeler, 2012]. While we have not conducted a comparison to symmetric discretizations herein, our experience from other work is that the nonsymmetric discretizations are mainly penalized due to the added memory usage of nonlinear iterative solvers.

The spatial discretization uses previous code for MPFA and MPSA discretizations on unstructured grids. These codes have been modified to provide the coupling terms appearing on the right-hand side of equations (18) and (19). For the case of the evolving fracture, the equations are rediscretized only around the fracture tip for each evolution of the fracture.

While it may not be appropriate to comment on computational performance of nonptimized research code, we note that resolving the MPFA and MPSA discretizations involve a local linear system associated with each vertex in the grid. This is in contrast to finite element methods where it is sufficient to evaluate a quadrature. In our implementation, calculating the discretization itself represents the main computational cost in the above examples.

6 Conclusions

We present a novel cell-centered finite volume approach to discretizing coupled flow and deformation of porous media. Through this contribution, we establish that it is possible to model hydromechanical coupling consistently within a single unified finite volume approach. The methodology is particularly adapted to problems where the flow is of primary importance, justifying the choice of finite volume methods for both flow and deformation. The compatible choice of discretization has several advantages from a practical perspective: (a) the data structures are equivalent for flow and deformation, simplifying the design of efficient solvers. (b) The discretization can naturally be developed within a single software framework, enabling fully coupled simulation. (c) Fractures aligned with cell boundaries can be naturally included through boundary conditions on internal surfaces.