Nonlinear finite volume discretization of geomechanical problem

Elliptic differential operators describe a wide range of processes in mechanics relevant to geo‐energy applications. Extensively used in reservoir modeling, the Finite Volume Method with TPFA can be consistently applied to discretize only a specific type of application under severe assumptions. In this paper, we introduce a positivity preserving Nonlinear Two Point Stress Approximation (NTPSA) based on the recently developed collocated Finite Volume scheme for linear elastic mechanics. The gradient reconstruction is different from the one used in Nonlinear TPFA, but a similar form of weighting scheme is employed to reconstruct the traction vector at each interface. The convergence of the scheme is tested with a homogeneous anisotropic stiffness tensor. The motivation behind the implementation of a new discretization framework in mechanics is to develop a uniform discretization technique preserving monotonicity for generic poromechanics applications.

for decades. 1 Different numerical techniques have been developed to integrate the balance equations. However, the Finite Volume Method (FVM) with Two Point Flux Approximation (TPFA) has become an industry standard because of its simplicity and robustness. Nevertheless, the method requires a so-called K-orthogonal grid 2 which is barely feasible to construct for real-field models. 3 Energy transition extends the range of geological settings and physical processes to be taken into account in subsurface reservoir modeling. Land subsidence, induced seismicity and the stability assessment of hydrogen or carbon dioxide storages require full-field geomechanical and poromechanical modelling. The economical feasibility of geothermal resources strongly relies on the accurate prediction of energy redistribution in the subsurface. Many of these applications consider essentially anisotropic reservoirs or require advanced gridding that can not be resolved consistently with TPFA. Especially when we consider stress approximation or look into realistic fluid flow scenarios, the grids are often highly distorted because of the presence of impermeable or highly permeable zones, for example, fractures or faults. These features complicate the geological model and introduce nonphysical errors to the dynamic solution. So, the conventional TPFA does not provide a representative approximation of fluxes.
Numerous methodologies are developed and implemented to tackle this problem for both flow and mechanics. Some of the first were Multi-Point flux approximation (MPFA) schemes 2,4-7 where we utilize more than two cells for approximating fluid flux across the interface between control volumes. A similar form of discretization technique, namely multi-point stress approximation (MPSA), is employed to solve the linear elasticity problem. [8][9][10][11] Although the research in the field of simulation of momentum equations by Finite Elements (FE) is more mature, the advantage of using an FVM for mechanics is the support of a wide range of polyhedral meshes including Corner Point Geometry 12 that is ubiquitously used in reservoir modelling. Another advantage is that the use of collocated grids simplifies the formulation, solution and implementation of coupled fluid mass, energy and momentum balances within a unified FVM framework. [13][14][15] While MPFA gives a representative solution to the flow problem, it is known to be conditionally monotone 7,16 and may violate the discrete maximum principle. This non-monotone behavior often takes the form of spurious oscillations in numerical solution across the grid. In particular, it is more pronounced when solving the problem of flow in highly anisotropic porous media. 6 A nonlinear FV scheme for fluid flow was introduced almost two decades ago 17,18 and developed further by various researchers. [19][20][21][22] In this formulation, the linear elliptic equation is transformed to a nonlinear form such that the scheme becomes monotone. The idea of nonlinear Finite Volume (FV) approaches is that the flux approximations should have non-negative coefficients in front of elliptic unknowns which are pressure. This ensures non-positive off-diagonal and non-negative diagonal elements in Jacobian which are inverted to obtain the solution of a discrete system of equations. Also when there was a fully heterogeneous domain that is, all cells surrounding a control volume having different permeability tensors, the methods in literature 23,24 proved to be unsuccessful to approximate the flux accurately.
Another kind of approach involves using vertex interpolation of unknowns to construct a positive basis for gradient reconstruction. 22,25 But considering vertices will introduce additional unknowns to the discretization equation which becomes computationally expensive. Along the lines of the above research, the concept of Harmonic averaging is developed 26 and successfully extended to multi-phase flow and multi-dimensional problems 27,28 in the past few years.
In this study, we introduce the Nonlinear Two Point Stress Approximation (NTPSA), a new type of discretization method in the FVM domain. In analogy to the formulation of nonlinear schemes for fluid flow, we propose the derivation of nonlinear schemes for elasticity problems. The approximation is based on the weighting scheme for two single-side approximations of traction vectors on an interface. The weights are calculated from a two-point constraint. To guarantee the convergence of the nonlinear scheme, we do a positivity-preserving approximation based on the connection-based gradient reconstruction. The proposed scheme results in a nonlinear discrete system of equations even for a linear elastic response that can make the computational cost of purely geomechanical simulation unreasonably high. In reservoir modeling, the flow and transport already introduce significant nonlinearity that may alleviate the higher costs of the geomechanical part in a fully coupled simulation. 29,30 Homogenization function approach 27 valid for heterogeneous anisotropic elasticity is derived and utilized in the conjunction with MPSA. We show that the homogenization function approach has the potential to improve the preciseness of the solution which is crucial for the representative description of faults in the subsurface.
Both NTPFA and NTPSA were implemented in the open-source Delft Advanced Research Terra Simulator (DARTS) framework. It is a general-purpose reservoir simulator that supports applications ranging from advanced compositional [31][32][33] to geothermal 34,35 and induced seismicity. 15 This provides a good basis for a coupled NTPFA and NTPSA discretization within a unified simulation framework for reservoir modeling.

Conservation laws
Let us consider the system of partial differential equations (PDEs) of the following form where is a vector, is a matrix. In the case of Equation (1) represents the conservation of fluid mass balance in porous medium and is porosity, is the fluid density, is fluid viscosity, is permeability tensor, is fluid pressure. In the case of Equation (1) represents the conservation of energy in porous medium and is rock density, , are fluid and rock internal energies respectively, ℎ is fluid enthalpy, is heat conductivity tensor, is temperature. In the case of Equation (1) represents the conservation of linear momentum in elastic body and is density, is a vector of displacements, ∇ and (∇ ) stand for the Jacobian of the displacements and its transpose respectively, ℂ is the rank-four stiffness tensor. It is worth to be mentioned that the assumption of = in Equations (2), (4) leads to elliptic PDEs with respect to unknown pressure and the vector of displacements correspondingly. Equation (3) includes both convective and conductive heat transport so that in the case = it is close to elliptic PDE for Peclet number (ratio between advective and conductive transport rates) close to zero.
Equations (2)-(4) and their extension to multiphase multicomponent fluid flow, transport and coupled poromechanics are extensively used in reservoir modeling. 1,36,37 Robust solution of these equations requires reliable numerical schemes capable to handle severe heterogeneity, arbitrary material anisotropy and advanced gridding techniques. Moreover, the consistent numerical scheme must preserve local conservation of fluid mass balance to guarantee the proper solution of transport problem. 38 FV schemes can fit these requirements. They are both locally and globally conservative, and can be used to resolve highly heterogeneous anisotropic medium with arbitrary star-shaped polyhedral meshes. 24 Among the different kinds of FV schemes, NTPFA represents a good balance between robustness and complexity. 28 Although NTPFA is well-developed, to the best of our knowledge it has never been extended to elasticity problems. Below we introduce the NTPSA for momentum balance in a solid elastic body.

Static elasticity
We consider a static elasticity problem in Ω governed by the momentum balance equation as where is a vector of volumetric forces, is the rank-two stress tensor represented by a 3 × 3 symmetric matrix. In the case of elastic infinitesimal deformation, Hooke's law defines the stress tensor as where ℂ is the rank-four stiffness tensor, denotes the infinitesimal strain tensor, = { , , } is a vector of displacements. We define traction vector to a surface with a unit normal as The other representation of the traction vector is also possible 11 where ∇ , ∇ , ∇ are gradients of the components of displacement vector stacked into 9 × 1 vector ∇ ⊗ , is 3 × 9 matrix defined in the right-hand side of Equation (8) where are components of the 6 × 6 stiffness matrix written in Voigt notations. Elasticity problem in Equations (5), (6) is subjected to the boundary conditions written in the following general form 11 where , are traction and displacements at the boundary, respectively, is an identity matrix, denotes the outer product of column to row . In addition, , , , are coefficients that determine the magnitude of their corresponding boundary conditions, while , are the corresponding condition values. Below we consider Dirichlet ( = = 1, = = 0), Neumann ( = = 0, = = 1) and roller ( = = 0, = = 1) boundary conditions.

DISCRETIZATION FRAMEWORK
The conservation laws represented by Equation (1) can be written in integral form and discretized according to FVM. It implies the representation of integrals of divergence term ∇ ⋅ as the fluxes = at cell interfaces with the help of Gauss formula as where subscript denotes the values defined at cell and subscript means that the value is evaluated at the interface between cells and , superscripts and + 1 define variables taken at a previous and current time step, Δ is a time step, is a cell volume, is an area of the interface between neighboring cells. The backward Euler scheme is used here for time integration.
Below we introduce a nonlinear FV scheme for the elasticity problem. Its formulation consists of three parts. First, we present two single-side approximations of the traction vector. Second, we introduce the method for the reconstruction of gradients of the displacements which are required to accomplish the approximations. Third, we develop the linear combination of two single-side approximations. An extra condition is required to define the coefficients in this combination. In a similar fashion to fluid flux approximation, in this paper, we restrict the stencil of the final approximation to two points. As a result, we obtain a scheme based on NTPSA.

Local problem
Let us consider the local problem for momentum balance in an elastic body. The continuity of displacements and tractions between cell 1 and cell 2 can be written as follows where 1 , 2 are single-side approximations of traction vector, 1 , 2 and are the centers of cells 1, 2 and interface respectively, ⊗ denotes the Kronecker product, ⊗ ( − 1 ) represents 3 × 9 matrix constructed as follows Equation (12) assumes piecewise-linear displacements within each cell while Equation (13) ensures the local momentum balance between cells. Below we use the approximations derived from Equations (12), (13).

Single-side stress approximations
The approximation of traction vector is essentially multi-point because the displacement gradients in Equation (7) or in Equation (8) have to be approximated along multiple different directions even in the case of isotropic stiffness. The MPSA was first proposed in literature. 9 The recent work, 11,39 facilitate the different formulations of MPSA that do not use the dual grid. Following this approach, we represent the traction vector as where the following decompositions are used where, , are 3 × 3 and 3 × 9 matrices, = [ ⊗ ( − )][∇ ⊗ ] denotes transversal projection of displacement gradients onto the interface, ⊗ denotes the Kronecker product, ⊗ ( 1 − 2 ) , ⊗ and ⊗ ( − ) represent 3 × 9, 9 × 3 and 9 × 9 matrices constructed in a similar way to Equation (14), 1 , 2 are projections of cell centers 1 , 2 onto the interface, the direction of normal vector such that ( 2 − 1 ) > 0, 1 , 2 are the distances from the center of cells 1 or 2 to the center of the interface . The projections are shown in Figure 1. Co-normal decomposition in Equations (15)-(18) was initially proposed for diffusion operator. 40 It splits the flux into harmonic and transversal parts that improve the behavior of the scheme with respect to the locking issue. 27,40 Moreover, this representation allows the transversal term to be eliminated providing the interpolation technique with positive coefficients, called harmonic averaging point. 26 Recently, the harmonic averaging point was extended to mechanics 11 and poromechanics. 14 However, the harmonic averaging points may lay outside the interface between cells and strictly positive interpolation may be not possible. 41 In this paper, we use connection-based gradient reconstruction. 11 We also propose the homogenization function approach for momentum balance to find out positive interpolation across interfaces.

Gradient reconstruction
The approximation of the traction vector requires the approximation of displacement gradients. The following equation can be derived from the continuity of displacements in Equation (12) and tractions in Equation (13) on the interface between cell 1 and 2 11 ( At the boundary interface with boundary conditions defined in Equation (10) the following equation can be written 39 where = ( (20), (21) represent an equation with respect to the gradients of displacements ∇ ⊗ 1 . These equations can be written in every cell for every interface either internal or boundary. At least three of such equations are required to reconstruct gradients in three-dimensional space.
The positivity-preserving approximation is required for the convergence of nonlinear FV schemes written for diffusion operator. 42 It was achieved by choosing the right way of gradient reconstruction among many alternatives. 27 However, in contrast to a scalar diffusion operator, the traction vector includes nine scalar diffusion terms as shown in Equation (8).
To preserve the positivity of all coefficients in the approximation of traction, we perform the reconstruction of these nine terms individually which means that we use a specific triplet of equations for every scalar diffusion term.
When the stiffness is constant over the grid, the reconstruction of gradients ∇ , ∇ , ∇ can be done separately as Equation (20) reduces to The interpolation over four cells in 2D.
Even though the gradient reconstruction of the vector of displacements can be reduced to the independent reconstructions of its scalar components, it is not enough to guarantee the positivity of coefficients in the final approximation in Equation (8) or in Equation (17). For a better explanation, see the gradient reconstruction for fluid fluxes presented in Appendix B. Nevertheless, to find out this reconstruction one may iterate over different triplets of connections. The homogenization function approach was used to find out the admissible triplet of connections in nonlinear flux approximation in literature. 27 The approach proposes the methodology how constructing an interpolation over multiple connections which significantly expands the range of variants.
The search for an admissible triplet of connections that provides the reconstruction of gradients with positive coefficients may become expensive and unsuccessful, especially in the reconstruction of displacement gradients. In the case of three completely independent gradients, we need to make sure that 3 coefficients per gradient are non-negative. In heterogeneous case, the gradients of components are not independent. It results in a 9 × 9 matrix produced by the triplet of Equations (20), (21) to be inverted. In the end, we need to have a 3 × 9 matrix of coefficients to be non-negative.
It is worth to be noted that to reduce these costs the positivity constraint can be weakened and the search for admissible stencil for the reconstruction of gradients can be formulated as an optimization problem. 28

Homogenization function for mechanics
The homogenization function refers to the development of auxiliary conditions that relate the displacement gradients in cell 1 to any arbitrary cell. For example, in a three-dimensional hexahedral mesh, we are no longer limited only by six neighboring cells to reconstruct the gradients. Instead, this technology can provide an interpolation over neighboring cells, neighbors of neighbors and etc. For arbitrary point ∈ 2 Equation (20) can be rewritten as follows 27 where denotes the normal vector to the interface between cell 1 and cell 2, 1 = ( − 1 ) is a distance between the interface and the center of cell 1, 2 and 1 , 2 are matrices calculated for the interface according to their definitions in Equations (18), (8) respectively. Figure 2 illustrates the interpolation.
Taking the derivative of Equation (23) with respect to provides the expression for the gradient of displacement in cell 2 where diag( ) denotes 3 × 3 diagonal matrix with components of on diagonal.
Let us consider cell 3 which shares a common interface with cell 2. The continuity conditions from Equations (12), (13) can be written for this interface as Repeating the same derivations we can come up with an equation similar to Equation (23) written for any point ∈ 3 where , denote the points at the plane of interfaces and respectively, ∇ ⊗ 2 can be found from Equation (24). This interpolation can be extended to an arbitrary chain of cells as shown in Figure 2.
The homogenization function for displacement gradients can be defined in a similar way to how it was done for pressure gradient in Equation (B7) in Appendix B. It provides a relationship between displacement gradients in cell 1 to any other cell which are connected over a set of interfaces Σ where = ( ). Now we can derive the homogenization function for neighboring cells from Equation (23)  and for the neighbors of neighboring cells from Equation (27) One can see how the homogenization function changes while interpolating over a larger number of interfaces. We can generalize the expression of the homogenization function for the interpolation over an arbitrary chain of cells 1, … , + 1 connected through the set of interfaces Σ = { , , … }: where index corresponds to the final interface in Σ and + 1 corresponds to the final cell. Note that is normal to the interface between cells and + 1.
This auxiliary condition can be incorporated into gradient reconstruction in multi-point or nonlinear approximation schemes for momentum balance. In the MPSA, we can use these additional auxiliary conditions in the gradients calculated as a least square solution 11 that is, instead of interface neighbors, we will now consider additional terms which result in a larger multi-point stencil and can provide higher accuracy of flux approximation.

Weighting scheme
We have considered two approximations of the traction vector above. Both the straight definition of traction in Equation (8) and the transversal part in Equation (15) depend on the gradient. In flow approximation the weighting scheme in Equations (B12), (B13) is used for semi-fluxes included single-side approximation of gradients. The same methodology can be applied to the approximation of the traction vectors. The weighting scheme applied to Equations (8), (15) can be written as follows = ( 1 − 2 ) + 1 1 + 2 2 , where represents the approximation of traction vector, 1 , 2 are the single-side approximations of traction vector, 1 , 2 are single-side approximations of the transversal term defined in Equation (16), 1 , are 3 × 3 positive semidefinite matrices of the coefficients of convex combination such that 1 + 2 = .Taking 1 = 2 = ∕2 results in the Average MPSA scheme (in analogy to Average MPFA 43 ) presented in literature. 11 For derivations below it does not matter which form of weighting either Equation (32) or Equation (33) is used. Below we use notation meaning that the same logic can be applied for .
The following representation of single-side tractions 1 , 2 can be used once the gradient reconstruction described in the previous section has been performed where 1 , 2 denote the sets of inner cells participated in the approximation of 1 , 2 respectively,̃1,̃2 denote the sets of boundary conditions defined in Equation (10) and contributed to 1 , 2 , , are 3 × 3 matrices that represent the single-side approximations of 1 , 2 respectively, they are obtained from the reconstruction of displacement gradients in cell 1 and 2 being substituted to either Equation (8) or to Equation (15), = , + , represents the right-hand side of boundary condition in Equation (10) contributed to the single-side approximations, 1 , 2 denotes the sum of Neumann boundary conditions contributed to the approximation of 1 , 2 respectively. Substituting Equations (34), (35) in Equation (32) and rearranging the terms we can come up with the following equation where the following notations are used We can construct a two-point nonlinear scheme in the same way as it was done for flux approximation, by eliminating all the terms in Equation (36) except those that are proportional to 1 , 2 . Assuming that the matrices of weights 1 , 2 are diagonal, the following linear system of equations completely defines them which solution is where ( ) 1 , ( ) 2 represent the components of 1 , 2 vectors, stands for a small regularization parameter. 27 The weights 1 , 2 defined in Equation (41) introduce nonlinearity to the scheme as they depend on unknown displacements through 1 , 2 . It leads to nonlinear discrete equations for the linear continuous problem. This artificial complexity allows the numerical properties of the scheme to be adjusted. Here we reduce the stencil of the scheme to two-point as it leads to a monotone convergent scheme for a scalar diffusion operator. 27 However, the constraint in Equation (39) can be chosen in a different way. For example, keeping only the terms with a single displacement component in each row in Equation (36) it is possible to split nine coupled scalar diffusion operators in Equation (8) into three decoupled operators. If a two-point stencil is also preserved within this constraint then it can be seen that the full matrix will be M-matrix. 44 Here we impose the constraint in Equation (39) that produces a two-point approximation, but each semi-flux in Equations (34), (35) may include contributions from all displacements components.
The positivity of the solution is required for the purpose of nonlinear convergence. 27 For those setups where negative values of any displacement component may be exhibited, we add a large constant to the corresponding displacement component as an initial guess. Dirichlet boundary conditions were modified accordingly.

Sparcity pattern
The matrix which is assembled by utilizing the coefficient matrices from Equation (37), and the connections gathered from the traction flux approximation Equation (36) form a two-point stencil. However, the resulting Jacobian matrix has a sparsity pattern similar to the MPSA scheme since all nonlinear entries are consistent with the MPSA stencil. This discretization is equivalent to a lower-order FE approximation (e.g., Brezzi-Douglas-Marini order 1 scheme. Equation (36) can be used to provide an approximation to the traction vector in a discrete system from Equation (11). The system of nonlinear equations needs to be solved with respect to displacements serving as nonlinear unknowns. Here we adopt the Newton-Raphson method to linearize the equations. The resulting linear system of equations on each nonlinear iteration can be expressed as where stands for a Jacobian matrix, represents the update in the vector of displacements over nonlinear iteration, denotes the vector of residuals. We repeat Newton-Raphson iterations until the method converges.

Convergence study
In the first part, we prove the convergence of NTPSA and Average MPSA discretization techniques by using the homogeneous anisotropic tensor test suggested in literature. 11 Let us consider the following reference solution of the elasticity problem in the cubic domain Ω = [0, 1] 3 ( , , ) = We substitute the reference solution and the stiffness matrix to the left-hand side of Equation (5) and use the result of the calculation as the vector of volumetric forces in the right-hand side of Equation (5). We also subject the domain to Dirichlet boundary conditions according to the reference solution. The solution of all three components of displacement can be seen in Figure 3.
We reconstruct stress tensor in the centers of cells. This is done by the procedure described in literature. 11 As we can see in Figure 4, both NTPSA and AvgMPSA schemes converge with the second order with respect to displacements and we observe superlinear convergence with respect to stresses. We use cubic mesh in these calculations. In Figure 4, is the inverse of the cube root of a number of cells in the domain in the log scale. (43) for displacement components (left) (center) (right).

F I G U R E 4
The convergence of numerical solution made with nonlinear and linear schemes to the reference with grid refinement. The second-order convergence for displacement vector (A) and superlinear convergence for stress tensor (B) is obtained.
The residual over NTPSA iterations is shown in Figure 5 in a semi-log scale. The structured cubic and unstructured extruded wedge grids of various resolutions are considered. The residual drops more than by two orders of magnitude every iteration. Once a cut-off tolerance 10 −10 (orange line) is reached nonlinear iterations are stopped. Figures 4 and 5 demonstrate good behavior of scheme in a homogeneous domain.

Compression and shear
Geological formations in the subsurface contain faults and fractures that may highly affect the hydromechanical response of porous media. Contact mechanics governs the behavior of fault and fractures and it requires special treatment. 15,[45][46][47] Therefore, the preciseness of stresses calculated at faults presents a high interest for applications. Let us consider the two-dimensional unit square domain ( = = 1) shown in Figure 6. The domain is subjected to shearing ( = 0.01) and smaller compressive ( = −0.001) displacements at the top boundary, the bottom boundary is kept fixed whereas the left and right sides are free. The plane strain setup is considered. Vertical and horizontal displacements are presented in Figure 7.
In this test case, we analyze stress profiles over the horizontal line in the center of the continuous domain shown in Figure 6. Normal and tangential components of the traction vector over the line are presented in Figure 8 for two types  Table 1. We compare the reference solution against the profiles obtained with three different schemes: Average MPSA with positivity-preserving gradient reconstruction, Average MPSA with least square gradient reconstruction and NTPSA. The magnitude of tangential traction is higher than the magnitude of the normal one which is consistent with displacements applied at the top boundary. Thus, the relative accuracy for tangential traction is noticeably higher than for normal traction. All considered schemes exhibit up to 3% deviation in tractions in this example. Even though regular square mesh produces much smoother traction profiles than the adaptive wedge mesh, it can not be used for more complex geometry.
In order to preserve a positive solution in the NTPSA scheme, we shift the reference state for displacements to {10, 10} to make sure that they remain positive over nonlinear iterations. NTPSA scheme takes 3 nonlinear iterations to converge to the residual of 10 −10 .

Homogenization for linear MPSA
The above results (e.g., Figure 8) clearly show that the NTPSA scheme does not provide a smooth traction profile in the solution of elasticity problem. An alternative option for obtaining a smooth traction profile is explored here. To achieve this, we apply the homogenization approach discussed in section 3.3.1 directly in the least-squares reconstruction. When we are performing gradient reconstruction using least-squares, we take an extended stencil that not only include immediate neighbors but also several adjusted cells which are in the vicinity of these neighbors. To make this study grid independent, we proceed to try the same three types of grids (exclude hexahedrons) as discussed in section 4.2.
The parameters and , listed in Table 1 represent the characteristic cell size of the whole mesh and in the vicinity of fault respectively. A lower value of , relative to indicates finer mesh near the fault compared to the rest parts of mesh. We also refer to levels of homogenization which indicate how many neighboring cells are being considered in gradient reconstruction. For example, the 2nd level of homogenization used in cell means that all cells belonging to the level of homogenization ≤ 2 including cell can participate in the gradient reconstruction in cell . For each level we choose among multiple variants of reconstruction to achieve the approximation with positive coefficients. The levels of homogenization are shown in Figure 9. The 0th level of homogenization represents the stencil for the least squares gradient reconstruction. Figure 10 demonstrate the resolution study of two numerical strategies with respect to a normal traction profile conducted for the adaptive wedge meshes. As resolution increases the calculated normal traction matches the reference solution at the tips of the meshed line first ( Figure 10A) and only then in the middle of the line ( Figure 10C). This can be explained by the distance to boundaries and grid refinement nearby the line. We observe homogenization method gives a slightly different solution than the least squares counterpart. Specifically, in Figure 10B and Figure 10C we can see that the homogenized solution is closer to the reference solution than the least squares method. Figure 11 illustrates a similar resolution study conducted with two square meshes. As expected, regular meshes produce smoother and more accurate solutions than adaptive meshes. In this case, the profiles obtained with homogenization coincide with the solution obtained with least squares, slightly deviating from the reference solution.  Next, we analyze the extent of homogenization and observe how increasing the level of homogenization will impact the solution profile. We consider four levels of homogenization after the least-squares stencil. Figure 12 shows the response we get as we increase the level of homogenization. We can see for the mesh which has 3398 cells, the initial least-square solution is noisy, but increasing the number of cells in the stencil gives us a better profile (especially at the endpoints), although it never actually comes as close as the refined solution in Figure 10C even after four levels of homogenization.
In the case of square meshes, the higher number of homogenization levels produces the solutions which practically coincide with the results obtained with zero level homogenization scheme shown in Figure 11. A deeper homogenization does not bring any accuracy benefit on these meshes as the solution is already precise enough.
Overall, homogenization behaves as expected that is, considering several other connections in our vicinity apart from neighbors is bound to give a better understanding of displacement gradient on specific cell-face pairs. This could also be extended to a heterogeneous case, which will be discussed in future work.

CONCLUSION AND DISCUSSION
Spatial discretization is a significant piece of the numerical scheme that determines the spatial accuracy of the simulation. Unconditionally stable linear TPFA is inconsistent to describe fluid flux on non-K-orthogonal grids. It also can not be used to resolve momentum fluxes. Instead, the linear multi-point approximations were developed and successfully applied to fluid flow, elasticity and poroelasticity. However, it is well-known that multi-point approximations are not monotone and may introduce some non-physical features to the solution. Nonlinear approximation allows these features to be avoided.
The main advantage of nonlinear schemes is that they provide flexibility in spatial discretization by taking a convex combination of different approximations instead of having a single approximation. The weights in such combination can be chosen from various premises like the preserving of monotonicity or discrete maximum principle. It is known that the solution to the elasticity problem does not satisfy the discrete maximum (minimum) principle. Instead, Korn's inequality 48 used to prove the existence, uniqueness and well-posedness of elasticity problem should be maintained. 10 In this study, we proposed a FV scheme based on positivity-preserving NTPSA for the anisotropic elasticity problem. To preserve the positivity of approximation, we developed a homogenization approach for momentum fluxes over an arbitrary chain of connections (cells). Moreover, the homogenization approach allows for improving the accuracy of approximation which can be extremely helpful in the presence of faults.
We demonstrated the convergence of the scheme in a homogeneous anisotropic domain. The extension of the scheme to heterogeneous cases makes the preserving of positivity challenging. In this case, the equations in gradient reconstruction are coupled which implies having much more coefficients to be non-negative for a given stencil. The potential solution can include relaxed requirements for the positiveness and reformulation of the search as an optimization problem. 28 Besides, the use of harmonic averaging points in gradient reconstruction can help. 11,14 We tested the accuracy of approximation with respect to the normal and tangential components of the traction vector. We found that the proposed nonlinear scheme does not provide better, precision than linear multi-point schemes. However, the use of higher-order spatial discretization techniques, for example, homogenization function approach, can increase the accuracy of the solution on adaptive meshes.
The proposed scheme can also form a basis to build a robust comprehensive FV-based reservoir simulator that integrates momentum balance with multi-phase fluid flow, compositional transport and energy balance in porous media.

A C K N O W L E D G M E N T S
Aleksei Novikov was sponsored by the project "Science4Steer: a scientific basis for production and reinjection strategies to minimize induced seismicity in Dutch gas fields" (with project number DEEP.NL.2018.046) of the research program "DeepNL" which is financed by the Dutch Research Council (NWO). We thank Matthias Moller, Hadi Hajibeygi and Martin Schneider for useful suggestions and discussions.

D ATA AVA I L A B I L I T Y S TAT E M E N T
All input and simulation data will be provided by request.
where we have two phases = , as water and oil respectively, stands for porosity, is phase density, denotes phase saturation, is the source (sink) of phase mass, is the phase Darcy's velocity, which is defined as where -relative phase permeability, -permeability tensor, -phase viscosity. For simplicity, we do not consider gravity and capillarity contributions in this paper, although they could be taken into account in nonlinear approximation. 24 We use the following definition of the phase flux and the flux where = ∕ denotes phase mobility. The equations above are supplemented with the boundary condition to be satisfied on the boundaries of domain Ω, where , , are prescribed constants that define boundary condition, stands for unit normal vector to surface. Two different types of boundary conditions, Dirichlet ( = 1, = 0) and Neumann ( = 0, = 1) are considered in Equation (A4). Moreover, we prescribe in the case of influx at the boundary.

APPENDIX B: FLUID FLUX APPROXIMATION
For flow problems, splitting the fluid flux defined in Equation (A3) into harmonic and transversal terms may help to avoid locking issues. 40 We use this decomposition, so-called co-normal decomposition, in the approximation of the fluid flux 27 where denotes the unit normal, , are the centers of cells = 1, 2 and the interface lying between cells 1 and 2, the direction of normal vector such that ( 2 − 1 ) > 0, = | ( − )| , = 1, 2 is a distance from the center of cell to the center of the interface , 1 = 1 + 1 , 2 = 2 − 2 are projections of cell centers 1 , 2 onto the interface, = , = ( − ) are co-normal and transversal components of permeability tensor , ∇ is a pressure gradient in cell 1. Note that the first term in Equation (B1) is called harmonic and the second is transversal.

Local problem
We repeat the derivation of NTPFA carried out in literature 21,24 but with the use of co-normal decomposition for flow representation. 27 We consider an interface between cells 1 and 2 and assume pressure remains piecewise linear. The continuity of pressure and flux across this interface can be written as where 1 , 2 are pressures in the cells, 1 , 2 are cell centers, 1 , 2 are permeability tensors defined at the cells = 1, 2, 1 , 2 are single-side approximations of flux. The point on interface , where requirements in Equations (B3), (B4) are imposed, can be the harmonic averaging point (HA) on interface, a concept introduced by literature. 17 In this work, we follow harmonic averaging equation specified in literature 24 to find the position and pressure related to the harmonic point.

Gradient reconstruction
The pair of Equations (B3), (B4) allows one of two pressure gradients to be excluded. However, another pressure gradient remains in the expression for fluxes. In order to derive an expression for pressure gradient, we formulate the auxiliary conditions which will help us to reconstruct it specifically for every cell-face pair. In three-dimensional space, three conditions are required to reconstruct the gradient. In this step, we make sure to choose the auxiliary conditions such that a positive basis is obtained for the scheme. The auxiliary conditions can be written as a system of equations with respect to ∇ 1 where is a 3 × 3 matrix of rows ( ℎ1 − 1 ) , = 1, 2, 3; ℎ , = 1, 2, 3 represents the HA points located anywhere in the plane of the interface, 24,26 ℎ , = 1, 2, 3 are the pressures in these points; 1 is the pressure in cell 1. Substitution of Equation (B5) in the left part of Equation (B4) gives where it is clearly seen that all components of ( 1 ) −1 have to be non-negative to have a flux approximation with positive coefficients. 24 In the case of heterogeneous anisotropic permeability on a nonorthogonal grid, the HA point may lay outside the particular interface and the gradient reconstruction with positive coefficients may run into difficulty. 41 Some formulations of NTPFA will require only a sum of coefficients in the flux approximation ( 1 ) −1 to be nonnegative which is a less restrictive condition. 22 Even in this case, the triplet of interfaces providing a positive sum of coefficients may not exist. In, 28 the searching algorithm was proposed to overcome this restriction at the cost of a wider stencil used in the approximation. The homogenization function approach provides a pressure interpolation over the chain of cells and may be used to find a triplet of interfaces that guarantees only positive coefficients. 27 Homogenization function  Σ 1, ( − 1 ) can be defined as a generalization of equations written in the system in Equation (B5) cabable to relate gradient in cell 1 to any other cell which are connected over a set of interfaces Σ  Σ 1, ( − 1 ) ∇ 1 = − 1 , where represents pressure defined at . Using pressure and flux continuity requirements from Equations (B3), (B4) over interfaces Σ homogenization function may be represented as follows 27 where +1 = +1 , 1 = ( − 1 ) is the distance between the center of cell 1 and the plane of interface defined by its center .
The homogenization function approach allows us to interpolate pressure over a heterogeneous domain and the gradient reconstruction is no longer limited by immediate neighbors to cell 1, but it can involve any cell. The stencil selection procedure will give priority to the cells which are encountered first that is, stencils with a lower level of homogenization and positive basis will be considered in the formulation. For example, = 0 level corresponds to immediate neighbors of the current cell