Cell-centered finite volume discretizations for deformable porous media

The development of cell-centered finite volume discretizations for deformation is motivated by the desire for a compatible approach with the discretization of fluid flow in deformable porous media. We express the conservation of momentum in the finite volume sense, and introduce three approximations methods for the cell-face stresses. The discretization method is developed for general grids in one to three spatial dimensions, and leads to a global discrete system of equations for the displacement vector in each cell, after which the stresses are calculated based on a local expression. The method allows for anisotropic, heterogeneous and discontinuous coefficients. The novel finite volume discretization is justified through numerical validation tests, designed to investigate classical challenges in discretization of mechanical equations. In particular our examples explore the stability with respect to the Poisson ratio and spatial discontinuities in the material parameters. For applications, logically Cartesian grids are prevailing, and we also explore the performance on perturbations on such grids, as well as on unstructured grids. For reference, comparison is made in all cases with the lowest-order Lagrangian finite elements, and the finite volume methods proposed herein is comparable for approximating displacement, and is superior for approximating stresses. © 2014 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.


INTRODUCTION
Coupling of fluid flow and mechanical deformation within a porous media is an important physical aspect of several surface and subsurface applications. Of particular interest for this presentation are engineered subsurface systems, such as CO 2 storage, oil-and gas production, and geothermal energy recovery. As concrete examples, the uplift at the In Salah storage site has potential to be an important monitoring mechanism [1], producing oil fields have experienced costly subsidence [2], and the long term poro-mechanical evolution, particularly with respect to fracture aperture changes, is critical for assessment of marginal geothermal resources [3].
Common simulation strategies for coupled fluid-mechanical systems rely to a large extent on specialized numerical discretization methods for the two individual problems of Darcy flow and mechanical deformation. The flow problem is most commonly treated by cell-centered finite volume methods (see e.g. [4] or [5]), and is implemented as such in the industry standard simulation software Eclipse by Schlumberger and Stars by CMG. In contrast the mechanical system is most commonly treated by finite element discretizations or with vertex-centered finite 401 The remainder of the paper is structured as follows. In the next section, we give the governing equations, the general cell-centered finite-volume framework, and discuss in detail the properties of various stress approximations. In section 3 we detail the construction of three different local stress approximations. The quantitative and computational aspects of the discretization are summarized in section 4, whereas comprehensive numerical examples are reported in section 5. The paper is then concluded in the final section.

FINITE VOLUME FRAMEWORK
In the most fundamental form, general equations of momentum balance are stated as conservation of momentum over an (arbitrary) volume. In Lagrangian coordinates this volume ! is considered fixed in space, and momentum conservation reads (see e.g. [24] for a clear introduction) Here D is the displacement from the initial state, T .n/ is the forces on the surface of !, identified by the outward normal vector n, and f are body forces acting on the material. The forces T .n/ may be referred to as surface traction or simply (surface) stress. In general the stress will be a non-linear function of the displacement, but for small deformations a linear approximation can be introduced. This linearization coincides with the symmetric Cauchy stress tensor , and we are therefore interested in the problem where The linearized strain is defined as the symmetric part of the deformation gradient, " D rD C .rD/ T =2 For a general material, the Cauchy stress tensor can be related to the strain through a fourth order tensor C through Hooke's law While potentially spatially discontinuous, the constitutive tensor functions C cannot be considered arbitrary, as they must satisfy certain symmetry and positive definiteness properties for the problem to be well posed. We will henceforth assume that the constitutive tensor satisfies the necessary conditions. Of special importance is the isotropic case, where the constitutive tensor has only two free parameters, and can be expressed as D C W" D 2 " C .tr "/ I Here the material coefficients and are known as the Lamé coefficients, and represent resistance to shear and compression, respectively. Most computational approaches to solving equations (1)(2)(3) or (1,2,4) involve first writing the integral balance equations (1-2) on a differential form by applying Gauss' theorem to the surface integral and subsequently noting that the integrands must satisfy (given sufficient regularity) The equation satisfies causality in time, and the temporal dimension can be discretized using a variety of standard methods, including an implicit Euler approach or Newmark time integration. Furthermore, in porous media applications the material inertia is usually negligible on the time-scale of fluid flow. In the remainder, we will therefore focus on the discretization of the spatial terms in equations (1) and (5). From the differential form given in equation (5), finite element approximations can be obtained by writing the equations on a weak form. Standard finite element approximations for the spatial terms in equations (4)(5) are well established, and are the method of choice for most mechanical problems [25]. We will therefore use finite element approximations as benchmarks in our numerical comparisons.
Our goal is in contrast to construct a finite volume discretization. By discretizing the integral equation (1) directly, we construct a method that exactly preserves the momentum balance. Such methods are known as finite volume, conservative finite difference, or control volume methods. Furthermore, by introducing a local (in space) approximation of equation (3), we obtain a method for which the stresses can be explicitly written in terms of displacements, such that the global system only has displacement as primary variable. This implies that we have only one degree of freedom per dimension per cell. We will term our local stress expression a multi-point stress approximation (MPSA), and construct the MPSA such that it is second-order consistent for general material properties and grids. The resulting method is directly applicable to non-isotropic and heterogeneous stress-strain relationships C , and is for range of problems second-order accurate for both displacement and stress. The method is applicable to both two and three spatial dimensions, structured and unstructured grids.
Finite volume approximations thus mimic the two components in the governing equations, in that they consider separately the momentum conservation and the constitutive law. Consequently, we will describe these two aspects in separate subsections.

Discrete momentum conservation
We partition the domain into a finite number of non-overlapping cells ! i . In a practice, this partition will in geological applications resolve the material discontinuities. For two cells i and j sharing a boundary, we denote the shared boundary e i;j . We are only interested in shared boundaries of dimension N 1, e.g. surfaces in 3D and lines in 2D, and consider e i;j as empty otherwise. Furthermore, keeping the 3D case in mind, we will refer to the shared boundaries in the following as cell faces. We are now equipped with the notation to write Equation (2) for each cell as By identifying D i and f i as the volume averaged displacement and force over the cell ! i and T i;j as the surface average stress over face e i;j , equation (6) can equivalently be written as We note that imposing a finite number of conservation volumes is the only approximation thus far, and that equation (7) is still algebraically equivalent to equation (1). In particular, if the right hand side is computed exactly for average surface and body forces, equation (7) exactly expresses the evolution of the average displacement of the grid cells.

Discrete stress
The finite volume method is completed through the choice of a discrete expression for the stress on cell faces. This mimics the physical modeling where the governing equations are closed through a constitutive rheology such as Hooke's law. Finite volume methods are commonly employed in subsurface simulation of flow, where the underlying constitutive rheology is Darcy's law, which gives the fluid flux. We will briefly review the qualitative aspects of various discrete flux expressions, and motivate the suitable framework to generalize to obtain an expression of discrete stress.
(1) Piece-wise polynomial interpolation. This approach considers a piece-wise polynomial interpolation (most commonly with a finite element basis) of the cell-centered displacement values. Consequently, it is known in various versions as the finite-element finite volume 403 method, box method [26], and control-volume finite elements [11]. Using this interpolation, the continuous constitutive law, equations (2)(3) are applied to obtain the stresses at cell boundaries. However, this approach requires that the material coefficients (e.g. C / are continuous over cell boundaries for the stress to be uniquely defined. However, in applications, the conservation grid is typically aligned such that discontinuities coincide with the cell boundaries. (3) Two-point approximations. On sufficiently regular grids, with isotropic material coefficients, it is sufficient to use the two neighboring cells to approximate the flux [4]. This leads to finite-volume/finite difference methods, see e.g. [14]. While this approach overcomes the problem of discontinuous coefficients, this approach fails to converge even in the scalar case in the presence of general grids or anisotropy [15]. (3) Multi-point approximations. More complex approaches to calculating local flux approximations are commonly referred to as multi-point approximations. In general, a separate expression is calculated for flux over each cell face, typically by considering subdivisions of the face [16]. While more expensive and complicated to implement, these flux approximations are known to be convergent for a wide range of grids, also when considering materials with large jumps in coefficients.
The two first approaches described above have previously been applied to mechanical deformation [11,14]. However, the multi-point approximations have to our knowledge not previously been extended beyond scalar equations, where they are known as multi-point flux approximations (MPFA). Given their superior qualitative properties, we seek to make such an extension herein termed multi-point stress approximation (MPSA).
In the scalar setting, there are three main types of MPFA methods. They are known by letters which mimic the connectivity of their approximation stencils O-, U-and L-methods. The three variants possess different qualities, and we will extend all three to mechanical problems. As such, we will define the MPSA O-, U-and L-methods below. Key references in the scalar case are [15] and [27] for the O-method, [15] for the U-method, and [28] for the L-method.

MULTI-POINT STRESS APPROXIMATIONS
Equation (7) is common for all finite volume methods, and the key aspect of a successful finite volume discretization lies in constructing an equivalent discrete constitutive relationship to equations (2)(3).
It follows from the continuous constitutive laws that we seek to express the surface stress T i;j as a linear function of the displacements D i . For the approximation to preserve the conservation property (7), we furthermore require that T i;j D T j;i , thus each face must have a unique, linear expression for the average stress, which we write explicitly as where t i;j;k D t j;i;k are referred to as stress weight tensors. A discrete stress-strain law on the form (8) is referred to as local if t i;j;k is non-zero only for cells k which are in some sense near neighbors to the face e i;j . As discussed in section 2.2, we will give the construction of the stress weights for what we term the MPSA. Here, multi-point emphasizes that the stress weights tensors depend on more than the two neighbor cells of the face (e.g. t i;j;k is non-zero for more than just k 2 i; j /. This is necessary in order to consistently approximate tangential derivatives to the face.

Geometry
The calculation of the stress weight tensors t i;j;k for face e i;j is split into several smaller calculations by partitioning the face into an equal number of subfaces to the number of corners on the face. Figure 1. The figures illustrate conservation cells with black lines, with cell centers marked by dots and grey lines delineating the subcells. The grey lines also partition the faces into subfaces. The three figures illustrate the grid structure for logically Cartesian (a), triangle (b) and unstructured grids (c). We will refer to unstructured grids such as the one illustrated, where no corner of the grid has coordination number (e.g. number of connected edges) higher than three, as dual-triangle grids, despite the dual grid not being strictly triangular. In the figures typical interaction regions are shaded (L-method, O-and U-methods, and all methods, respectively).
In 2D, this implies that every face is split into two subfaces, while in 3D the number of subfaces depends on the geometry and may in general change from face to face. The volume associated with corner l and the cell i will be termed a subcell, and denoted Q ! i;l see figures 1 for 2D examples. For concreteness, we say that the face is split such that the subface corners align with the centers of the face and edges. We will refer to the subfaces as Q e i;j;l , where l is the cell corner associated with the subface, and hence a common corner of cells i and j The stress weights for subface Q e i;j;l are denoted Q t i;j;l;k , and for consistency Considering a corner l and an adjacent subface Q e i;j;l we construct three local volumes, which will be of interest to the three different stress approximations, respectively. The local volumes will in general be referred to as the interaction regions.
Definitions, interaction regions: The interaction regions for the methods are defined as follows For the MPSA O-method, let the interaction region consist of all subfaces and all subcells adjacent to the corner l (see upper shaded area in figure 1b). For the MPSA U-method, let the interaction region for subface Q e i;j;l , consist of the subset of the MPSA O-interaction region which is neighbors or neighbors-neighbors to subface Q e i;j;l (see lower shaded area in middle figure 1b). For the MPSA L-method, let the interaction region for subface Q e i;j;l , consist of the subset of the MPSA U-interaction region retaining only the N C 1 closest subcells, where N is the physical dimension of the problem (see figure 1a).
The interaction regions are illustrated in figures 1. Notice that for dual-triangle grids, the interaction regions for the three methods become identical, as seen in figure 1c). For general 2D grids, the interaction region for the O-method will be the dual cells, while the interaction regions of the U-method and L-method form stylized "U" and "L" shapes straddling the subface, hence the names.
For the MPSA L-method, the notion of "closest subcells" must be made more precise, herein we use the simple Eucledian distance weighted by the square root of the second Lamé coefficient, which is the natural extension of the scalar selection criteria [28].
We will construct the subface stress weights Q t i;j;l;k within each interaction region, which together with equations (7-9) will define the discretization method. We can already see that the interaction region thus defines the size of the support of the discretization stencil, thus in terms of Equation (8), we see that each full stress weight t i;j;k will be zero unless cell k shares at least a corner with face e i;j . Our discretization thus has local stress weight tensors.

Calculation of stress weights
The stress weights for each subface Q e i;j;l are in principle calculated separately. For the calculation of the stress weights for subface Q e i;j;l , consider a linear approximation to displacement within each subcell. Then within each subcell Q ! i;l we approximate each component of displacement by a multilinear function of the spatial coordinates, such that Here we recall that D i is the cell-average displacement, which we associate in the linear expansion with the cell center x i . We constrain the degrees of freedom in the interaction region by two forms of continuity.

Continuity of surface stresses:
Since G i;l is constant within each subcell, the stress tensor calculated from Equation (3) will also be constant within each sub-cell. For planar subfaces, we thus impose pointwise continuity of surface stress, while for non-planar subfaces, as can result when 3D faces have more than three corners, we impose only integral continuity of surface stress. In either case, the continuity of stress can be expressed as

Continuity of displacement:
Full continuity of both surface stress and displacement will in general lead to an overdetermined system for linear approximations such as given in Equation (10). Instead, we will require continuity of displacement only for specific points Q x i ;m;l;n on Q e i;m;l : If a single continuity point Q x i;m;l;1 on the subface Q e i;j;l is chosen, we refer to this constraint as weak continuity of displacement. If N points Q x i;m;l;n , for n D 1 : : : N , are chosen on the subface Q e i;j;l , we refer to this constraint as strong continuity of displacement. For subfaces with a weak constraint, the location of continuity point may impact the method, (see [27,29] and [30] for discussions of this issue in the scalar case).
The above definition of interaction regions and constraints allow us to define the various MPSA methods.
Definition, MPSA L-method: The MPSA L-method is obtained by constraining the degrees of freedom of Equation (10) by requiring stress continuity and full displacement continuity over all subfaces in the interaction region of subface Q e i;j;l .
Definition, MPSA U-method: The MPSA U-method is obtained by constraining the degrees of freedom of Equation (10) by requiring stress continuity over all subfaces in the interaction region of subface Q e i;j;l . Furthermore, weak continuity of displacement is imposed over subface Q e i;j;l itself, while strong continuity of displacement is imposed over the remaining subfaces in the interaction region. In total, the above constructions yield linear systems for the unknown gradients G j;l within the interaction region for edge Q e i;j;l The local stress weights are defined from the gradients G j;l using either term in equation (11).
These local linear systems have equal number of constraints as degrees of freedom, and in general they are thus solvable. However, in certain special cases of practical importance, some constraints may become redundant, rendering the local systems without a unique solution. We note the following such counter-examples: A) The L-and U-method leads to underconstrained systems if the cell-centers align on a hyperplane (line in 2D or surface in 3D). B) The O-method leads to an under-constrained system if the primal grid has orthogonal grid lines relative to the anisotropy imposed by the constitutive tensor C . These counter-examples have important practical consequences. In particular, counter-example A) renders the U-and L-methods largely unsuitable for triangular grids, where it is common that cell centers are close to aligning. In contrast, the O-method is applicable for (strict) Delauney triangulations, but requires stabilization for logically Cartesian grids which frequently contain 90 ı angles. We provide this stabilization here.

Definition, MPSA O-method (general):
The MPSA O-method for general grids is obtained by constraining the degrees of freedom of Equation (10) by requiring stress continuity over all subfaces in the interaction region for subface Q e i;j;l . The remaining degrees of freedom are defined by the weighted least squares deviation from strong displacement continuity over all subfaces in the interaction region, where the weights are chosen as the harmonic mean of the largest eigenvalue of C of the adjacent cells.
We close this section by commenting on the choice of continuity points where weak continuity is imposed. In particular, it is known from the scalar case that for simplex grids in 2D and 3D, the MPFA O-method is symmetric if Q x i;m;l;1 is located 1/3 of the distance from the face center to the corner, as the effective part of the subcells then becomes parallelograms and parallelepipeds [30]. The same property holds for MPSA O-method, and for simplex grids we therefore prefer this version of the MPSA O-method over the general form. For the MPSA U-method, the choice of continuity point is less critical, but our numerical investigations have indicated that locating Q x i;m;l;1 at the face center provides the more robust results.

Boundary conditions
Implementation of periodic, Neumann, or Dirichlet boundary conditions is straight forward. If Q e i;j;l is a on the domain boundary, the corresponding lines in equations (11) and (12) are replaced by the given boundary conditions. If any other subface Q e i;m;l is on the boundary .m ¤ j /, then again the corresponding lines in equations (11) and (12) are replaced by the given boundary conditions, noting only that for Dirichlet boundary conditions, only weak displacement continuity can be enforced to avoid an over-determined system.

PROPERTIES OF THE MPSA METHODS
Herein we discuss a priori properties of the MPSA methods that can be deduced from their construction. In particular, we consider the computational cost in terms of both determining the discretization stencils as well as for the global system. Furthermore, we discuss the relative merits of node-centered vs. cell-centered variables.

Computational cost
The computational cost of the MPSA methods is in general dependent on two factors. The local systems in the interaction regions are of finite size and parallelize trivially, while the global system is relatively large with less possibilities for parallelization. We discuss these two issues below.
The MPSA U-and L-methods require one local calculation for each subface. The MPSA Omethod will have an identical local system for each subface in an interaction region, and thus requires only one local calculation for each cell in the dual grid. Let K be the number of subcells in the interaction region, and L the average coordination number for each corner of the grid. In 2D, K D 3; mi n.L; 4/; L for the L-, U-, and O-methods. The computational cost of the local solves depends on a linear system proportional to K, which needs to be solved for N K right-hand side for each of the L subfaces (for the L-and U-methods). Assuming an efficient solver with linear scalability, and recalling from equation (10) that there are N 2 unknowns in each subcell, we can calculate the relative computational cost of the local calculations per grid cell. In 2D, these are of order N 3 K 2 L for the L-and U-method and N 3 L 2 for the O-method.
The computational cost of the MPSA stress weights should be compared to the inner products for finite elements. For finite elements, no local systems need to be solved, and instead the inner products involving transformation matrices typically need to be evaluated at a small number of quadrature points within each element.
The second aspect of the computational cost relates to the global linear system. The MPSA methods lead to the smallest global linear system possible for a given grid, with one vector variable per cell. This is typically important in geological porous media where a large grid is specified to account for the domain size and material heterogeneity. Furthermore, the bandwidth of the system matrix is restricted to cells that share an interaction region. However, it is well known that finite volume methods typically do not lead to symmetric linear systems, since the displacement changes do not a priori need to yield reciprocal effects on the stresses on the (in some sense arbitrary) faces of the grid. Thus for large linear systems recursive iterative solvers such as CG have to be replaced by more expensive and memory consuming solvers such as GMRES. This loss is to a large extent mitigated for coupled problems, where a choice of finite volume discretization of the flow equations will lead to a non-symmetric component.
The exception is the MPSA O-method for triangular grids, which leads to a symmetric system matrix when the continuity point of displacement is chosen as described in section 3.2. Finite element discretizations also lead to symmetric linear systems, and the lowest-order finite elements have a comparable number of degrees of freedom to the finite volume methods.

Cell-centered variables
Common node-centered finite volume methods and finite element methods compute node-centered values of displacement. In contrast, the finite-volume approach advocated herein has cell-centered variables that are appropriately interpreted as the average displacement of the grid cell. It is worthwhile to briefly address the relative merits.
The node-centered variables, together with the function-space framework of finite elements, provide a natural point-wise interpolation of the solution vector from the nodal points to the full physical domain. In particular, this interpolation leads to continuous displacement, and a stress tensor that is uniquely defined almost everywhere, with the exception of cell faces.
In contrast the cell-centered variables, do not have a native interpolation, and represent the average displacement of the cell. This is sensible from a conservation-law point of view, but provides a weaker notion of the point-wise solution. Similarly, only the surface stress is recovered, defined on cell faces. This sparser set of information does come with advantages, in that it is typically exactly at the cell faces, aligned with material contrasts, where the aligned surface stress may be of greatest interest.
For both the node-centered and cell-centered variables, the deficiencies mentioned above can largely be rectified through post-processing. The main advantage of the cell-centered variables thus remains that the data structure coincides with the data structure typically used for the fluid flow equations in geological media. In many applications, the flow is of primary importance (e.g. oil recovery, CO 2 storage, or geothermal heat extraction), and the development of integrated codes for flow and deformation will therefore be facilitated by allowing for the deformation to be described on the flow grid.

Divergence of displacement
We note that it is straight forward to obtain a local stencil for the divergence of displacement based on the local calculations in Section 3.2. Indeed, the local calculations are conducted for each subface, and give approximations to the gradient of displacement in each of the two adjacent subcells. For the MPSA O-method, this is sufficient to obtain the approximation to the divergence of displacement.
For the MPSA U-and L-methods, the each subcell will have contributions from N subfaces that arise from different local calculations. Thus each subface calculation will only contribute to the divergence of displacement proportional to the triangle (pyramid in 3D) connecting the cell center and the subface.

NUMERICAL RESULTS
Our numerical examples aim at providing an assessment of the performance of the three MPSA methods in terms of heterogeneous coefficients and varying Poisson ratios on challenging structured and unstructured grids. As a benchmark, we will use the lowest-order Lagrange finite elements, chosen based on comparable computational cost and complexity.
We structure the presentation of the space of test problems around the following four axes. First, we consider the sensitivity to jump in material properties. Secondly, we consider the sensitivity of the methods to the Poisson ratio. Thirdly, we explore the robustness with respect to different classes of irregular grids. The first three sections are all conducted on logically Cartesian grids, and we finally make a comparison also on unstructured grids.
For all numerical tests, we consider the unit square or unit cube with zero Dirichlet boundary conditions (no displacement, or "clamped" boundary conditions) and choose the parameters and load functions such that the solution is a known analytically. We mention three error measures of relevance for applications. First, we consider the error in the displacement. Secondly, the divergence of displacement is the key coupling term between flow and deformation in the Biot equations for deformable porous media, and we therefore also monitor this measure. Finally, the stress acting on internal surfaces is important from the perspective of material fracturing, which is common in geological porous media, and we therefore also consider the error in the stress on all internal cell faces. For all examples considered herein, the stress and divergence of displacement are qualitatively and quantitatively similar, and in the interest of space only the stress will be reported in the figures.
For all MPSA cases, the displacement error of the discrete solution D h is taken as the relative weighted error measured in the discrete L 2 norm as compared against the cell centered displacement D i Á D .x i /: The error in the discrete divergence is calculated analogously. Similarly, the error in the discrete stress T h is again measured as the relative weighted error in the discrete L 2 norm as compared to the stress at the face center T i;j Á n i;jˇx i;j : For the finite element calculations, the displacement variable is weighted by the volume of the dual grid, while the full stress tensor at the cell center, weighted by the cell volume, is used in lieu of a natural surface stress variable.
All numerical experiments in section 5.1-5.3 are based on grids ranging from 3x3 to 96x96 cells. A 3x3 grid was chosen as basis for refinement to align with the material coefficients in Section 5.1 and to resolve the structure of the Kershaw grid in Section 5.3. For the comparison to the finite element method (FEM), a subdivision of the square cells into triangles was introduced according to the shortest diagonal. For the numerical examples herein, the generalization of the MPFA O-method calculates the largest eigenvalue of C based on the Lamé coefficients. The implementation is in Matlab, and uses the standard sparse libraries and the direct solvers provided therein. Furthermore, the FEM discretization used for comparison purposes is implemented using the Matlab PDE toolbox. All calculations are performed on a standard laptop.

Jump in material properties
As discussed in section 2.2, the MPSA approach is expected to accurately honor discontinuities in material properties. We investigate the properties of the proposed method in this respect.
Let the function .
x/ indicate a heterogeneity in the "middle" block of a three-by-three partitioning of the unit square, e.g.
Then for a material with discontinuity characterized by the constant Ä, we choose to consider Lamé coefficients with the structure For this material, we consider the solution which honors the zero Dirichlet boundary conditions as illustrated in Figure 2. This solution is continuous, with continuous stresses. In particular, we note that the stresses are independent of Ä. We find the body force by taking the divergence of the stress, and solve the problem with the various discretizations. The accuracy of the method is assessed while varying Ä from 10 6 to 10 6 , including the homogeneous case where Ä D 1. For this example, we have considered grid with uniform squares for the MPSA methods.  second line is the homogeneous case and can be considered the benchmark for the first and third lines of figures. We see that all methods converge as second order for displacement independent of the discontinuity in displacement. The MPSA L-method appears to converge asymptotically only as first order for all cases, including the homogeneous case, this is due to the indeterminacy in the local stencil for perfectly symmetrical grids such as the square lattice chosen here (see also next section). The first data point for the FEM is missing for the displacement error as the analytical solution is zero in the quadrature points for the coarsest grid, and the relative error is therefore undefined. In terms of absolute values of relative error in displacement, the finite volume methods appear to marginally outperform the finite element method for this example.
For stresses (and also divergence of displacement), the MPSA U-and O-methods are second order convergent, while the MPSA L-method and FEM are first order convergent. Again, all results are independent of material contrast. We note that the second-order convergence rate in stress enjoyed by the MPSA U-and O-methods is only possible when the continuous solution has sufficient regularity, such as in this example.

Poisson ratio
Numerical methods for deformation on quadrilateral grids may suffer from a locking of the degrees of freedom in the limit of incompressible media (see e.g. [31] for a recent study). This is commonly expressed in terms of the Poisson ratio This ratio is physically bounded above by 0.5, although most geophysical materials have a Poisson ratio between 0.06 (marble and granite) and 0.48 (soft clay). Nevertheless, it is also of importance to consider the limit as the Poisson ratio approaches 0.5, as non-linear models for deformation are incompressible in the failure limit [32].
In order to avoid the blow-up of the right-hand side as ! 1, we will in this section consider a divergence-free displacement field These solution vectors and associated forcing functions are shown in figure 4. We consider a material with unit shear modulus of D 1, and explore the accuracy and convergence of the methods for varying Poisson ratio by varying the second Lamé coefficient For this example, we have used a curvilinear grid as prototype of deformed reservoir models away from faults For all methods we use the smooth grid sequence given as sequence B in section 5.3. Figure 5 shows a representative sample of results for this test case. The first row represents the homogeneous case with unit Lamé coefficients, while the second and third rows have the second Lamé coefficient equal to D 10 and D 10 3 , respectively. This is equivalent to Poisson ratios of D ¹0:25; 0:45; 0:4995º, respectively. We see that the performance of all methods is more sensitive to the Poisson ratio than to material contrasts. The methods were also tested for lower (and negative) Poisson ratios, but the performance was similar to that for the Poisson ratio equal to 0.25.
For the error in the displacement, all methods ultimately approach an asymptotic convergence rate of second order. However, only the MPSA O-method remains unaffected by the Poisson ratio. We note that already for the homogeneous case, the MPSA U-method suffers from the relatively skewed grids of this example (see figure 6a), however the MPSA L-method performs better on these grids than on the square lattice used in the previous section. As the Poisson ratio increases, the MPSA O-method retains second order convergence for all grids, while the remaining MPSA methods and the finite element method all experience a plateau of reduced convergence rates. This   plateau can be inferred from the second row where D 0:45, and becomes evident for D 0:4995. For even higher Poisson ratios, the behavior of methods becomes further emphasized. In terms of absolute value of the displacement error, the methods are comparable for low Poisson ratios, and in the asymptotically convergent regime, while for high Poisson ratios the plateau in convergence appears at a lower accuracy for the finite element methods. The convergence results for stresses follow the same pattern as displacement, with the stresses from the MPSA O-method being the only robust stresses for high Poisson ratio.

Robustness on challenging grids
The previous two sections have considered regular square lattice grids and smooth grids. In this section we explore the suitability of the method on more general grids, obtained from transforming the original grids. Thus we only consider grids that retain the logically Cartesian structure, but we allow for arbitrary quadrilateral grids.
We consider four grid sequences, as illustrated in figure 6 at refinement level 3 (12x12 grid cells). Grid sequence A is chosen for reference as the standard square lattice. Grid sequence B is a smooth transformation of the square lattice grid, using the transformation .x; y/ ! .x C 0:5x .1 x/ sin 2 y;y C 0:5y.1 y/ sin 2 x/ While highly idealized, this grid sequence mimics the grids associated with geological porous media that historically have undergone continuous deformation. Grid sequence C is non-smooth transformation (Kershaw-type grid [33]), and is defined by dividing the unit square by a piece-wise linear function through the points .0; 0:2/; .1=6; 0:2/; .1=3; 0:8/; .2=3; 0:2/; .5=6; 0:8/, and .1; 0:8/. The area above and below this function are sub-divided equidistantly. While the Kershaw grid was originally motivated by unstable flow fields, we consider it as an example of highly skewed, contrasting, grid cells.
Grid sequence D is obtained by randomly perturbing the square lattice grid. The perturbation of each coordinate of each corner is a uniformly distributed number between˙x=2. This guarantees no non-convex grid cells. The grid sequence is chosen to highlight the robustness of the methods, emphasizing the convergence properties when the shape factors of the grids to not improve with grid refinement (e.g. when the refined grids do not converge to parallelograms). This notion of rough grids is of relevance in applications where grid refinement cannot be afforded [34].
For all grids, we use the intermediate Poisson ratio of D 0:45. To ensure that the results are not artificially impacted by using pure sine and cosine solutions, we modify equation (14) by a polynomial coefficient: The results are shown in Figure 7, where we have shown the results for grids A, C, and D (grid B is identified as the second row of Figure 5). All methods are convergent in displacement for all grids. Note that the decrease in accuracy on the Kershaw grid is due to that the coarsest grid, with only 3x3 cells, does not fully resolve the sawtooth pattern, and thus the refined grid (6x6) is more challenging. In terms of convergence rates, only the MPSA O-method retains full second order convergence for all grids. The MPSA L-method again suffers on the square lattice, and also for the highly skewed grid cells of the Kershaw grid, and only converges as first order for these grids. The MPSA U-method and the FEM perform very similarly on grids A, C, and D, with both struggling to obtain full second order convergence on the Kershaw grid, but otherwise showing robust qualities.
In terms of stresses and the divergence of displacement, the story is similar to displacement. On grids A and B, the methods perform as expected from sections 5.1 and 5.2. The MPSA O-method is the only method to retain second order convergence in stress on the Kershaw grid, while the MPSA L-method is not convergent. Finally, the randomly perturbed grid shows the challenge it poses with persistent lack of shape regularity, in that no method attains more than first order convergence rates in stress on this grid, although absolute value of the error appears smaller for the FV methods than the FEM method on this grid.

Extension to general grids
The FV-MPSA methodology presented herein is applicable to arbitrary grids. As pointed out in Sections 3 and 4, the three MPSA methods have identical interaction regions, and largely coincide for unstructured grids where the dual grid contains three cells. In contrast, for triangular grids the L-and U-methods are unsuited due to the degeneracy of the local systems. For these triangular grids, the O-method is still applicable as it leads to a symmetric system matrix which has favorable stability properties. We will verify these assertions here. For ease of comparison to the results on square grids, we return to grid sequence B from sections 5.2 and 5.3, and keep the same analytical solution as above.
To obtain a triangular grid, we use the same approach as for the finite element method comparisons of the previous sections, and subdivide each grid cell along its shortest diagonal. To obtain unstructured grids with triangular dual-grids, we invert the dual and primal configurations of the triangular grids. The resulting grid structures are illustrated for grid refinement 1 in figures 1. For consistency we use the MPSA O-method on all grids.
The results are shown in Figures 8. The grid types represent grids varying number of subcells in the cells of the dual grid, where the triangular grid has the most and the dual-triangular grid the least. The complexity of the method (and also the size of the discretization stencil) increases with the number of subcells in the interaction regions.
The MPSA O-method is convergent in both displacement and stress for all the grid types considered. However, the method performs better for grids with simpler interaction regions, i.e. the dual-triangular and quadrilateral grids. For the triangular grids, the convergence rate is reduced and only first order convergence is observed for displacements. We infer from these studies that the MPSA methods in general are more robust and favorable on grids with few subcells in the in the cells of the dual grid. Notably, this implies that it is advantageous to consider grids with low coordination numbers.

Computation of a solution with singularity
As a final validation, we consider a grid with a crack entering the domain. In particular, we are interested in the geological setting with fluid-filled fractures, thus internal to the crack we will assign Neumann boundary conditions mimicking a fluid of fixed pressure acting on the rock. For simplicity, we take all the Lamé parameters D D 1 This is implemented by considering the rectangle . 3; 3/ . 1; 1/, where the boundary also extends into the domain on the line segment which forms the negative x-axis, e.g. . 3; 0/ to .0; 0/. Thus the domain has essentially a U-shape. The solution on such domains will have a singularity at the origin [35]. For this case we assign zero Dirichlet (clamped) boundary conditions on the three sides without a crack, and zero Neumann (no stress) boundary conditions on the boundary of the rectangle where boundary enters the domain. On the boundary that enters the domain, we have chosen Neumann boundary conditions with T n D 1, where n is the outer normal vector to the domain. We thus expect the crack to open. The computed solution is shown in Figure 9, based on both uniform Cartesian grids and a Cartesian grid with cubic increasing grid spacing around the origin. We identify three characteristics of the solution: As emphasized by the logarithmic plot, the crack tip is consistent with the asymptotic . x/ 1=2 scaling which is known for this problem [36]. Away from the tip the crack width reaches a plateau value of 1=3, which is consistent with the analytical solution to 1D compression with the prescribed boundary conditions. Finally, towards the boundary of the domain there is an influence of the no-stress boundary condition.

Limitations
The results presented in the preceding four sections document that the proposed FV-MPSA methodology is applicable for general grids on in the presence of both material contrasts and moderately high Poisson ratios. For the experiments considered, the finite volume methods were competitive, or better, than the finite element method of comparable complexity except on the triangular grids. Furthermore, the results show that while the convergence rates deteriorate in some cases, the new methods are in general robust for a wide variety of grids, Poisson ratios and heterogeneity contrast. To complement these studies, we document in this section limitations we expect from the new method, which were not seen in the above studies. It is known in the scalar case that the FV methods suffer from oscillations highly skewed parallelogram grids [29]. In the scalar case, this is a fundamental limitation of the methodology, and we expect this limitation to also be present for the extension to elasticity. This may explain the poor results for the MPSA L-method on the Kershaw grids. For problems with either a very strong anisotropy, or significantly more skewed grids than those investigated in section 5.3, we expect also the other MPSA methods to suffer.
It is also known from a degree-of-freedom counting argument that the finite element method locks on triangles, as was observed in section 5.2 [32]. Conversely, the MPSA methods (including the O-method) are expected to lock for high Poisson ratios when the dual grids are logically triangles (grid type C in Figure 1). This was not revealed by the preceding numerical examples, but has been verified in separate test cases.
While we have strived to put the novel FV-MPSA methods for material calculations through a rigorous set of tests, numerical experiments will never be exhaustive. Due to the difficulty in establishing general frameworks, convergence proofs for finite volume methods in the scalar case are typically restricted in validity. For the scalar MPFA methods convergence proofs typically rely on analogies to either mixed finite element methods [34], or mimetic finite difference methods [37]. However, convergence proofs have also been established in the finite volume setting directly [38]. With the understanding gained from the numerical investigations herein, selecting the appropriate setting in which to develop a convergence theory for the MPSA methods will be an ongoing topic of research.
While we have shown no results for 3D, the method has been implemented and tested for Cartesian grids also in this setting. Preliminary testing indicates that the 2D results are representative for the performance of the method in 3D.

CONCLUSIONS
A novel set of cell-centered finite volume discretizations for the equations for deformation in the limit of linear elasticity has been presented. By construction, the MPSA methods are locally conservative in terms of the momentum conservation law. The methods are distinguished by their approximation of the stress on cell faces. The approach herein generalizes the approach for scalar conservation equations, and adopts the terminology of O-, U-and L-method approximations.
The methods contain the minimum of one degree of freedom per component of displacement per grid cell. This is similar to lowest-order finite elements, and numerical results show comparable accuracy for the displacement variable and superior accuracy in stress for a wide range of problems.
The numerical examples emphasize the robustness of the method with respect to both jumps in material coefficients, and with respect to variations in the Lamé parameters. Structured and unstructured grids have been investigated, including cases with challenging grids.

417
The three variants of the method presented herein have different characteristics. From our numerical examples, the O-method performs best overall, and also relies on only a single calculation per interaction region. The U-method appears to be also generally applicable, but locks for high Poisson ratios. While the L-method has the advantage of the simplest discrete stencil, it is both less accurate and less robust than the two other variants.
The MPSA methods are all formulated for general grids, both in 2 and 3 dimensions, and as such are directly applicable to a wide range of problems. In particular, due to its cell-centered finite volume structure, we expect the method to be particularly attractive for problems involving coupled flow and deformation in deformable porous media.