Hybrid mimetic finite-difference and virtual element formulation for coupled poromechanics

We present a hybrid mimetic finite-difference and virtual element formulation for coupled single-phase poromechanics on unstructured meshes. The key advantage of the scheme is that it is convergent on complex meshes containing highly distorted cells with arbitrary shapes. We use a local pressure-jump stabilization method based on unstructured macro-elements to prevent the development of spurious pressure modes in incompressible problems approaching undrained conditions. A scalable linear solution strategy is obtained using a block-triangular preconditioner designed specifically for the saddle-point systems arising from the proposed discretization. The accuracy and efficiency of our approach are demonstrated numerically on two-dimensional benchmark problems.


Introduction
Modeling hydro-mechanical coupling is essential to accurately simulate a wide range of subsurface processes involving fluid flow and mechanical deformation, such as oil and gas recovery [1], geological CO 2 storage [2], and geothermal energy production [3]. In life sciences, this coupling also plays a central role in the modeling of bone deformation [4] and blood-vessel interaction in hemodynamics [5]. Studies involving the numerical solution of Biot's equations of poroelasticity on a computational mesh are routinely used to investigate these processes. In most subsurface applications, generating a mesh that faithfully represents the structure of the porous medium is a difficult task. Geological formations often exhibit a high heterogeneity characterized by stratigraphic layering and the presence of faults and fractures. As a result, there is strong interest in using unstructured polyhedral meshes that conform to the complex geological features of the porous medium [6]. To preserve accuracy and reduce computational cost, it is appealing to solve both the porous flow problem and the mechanical problem on the same mesh. In this work, we address this issue by developing a robust numerical scheme and fully coupled solution strategy for the displacementvelocity-pressure formulation of Biot's equations on arbitrary polyhedral meshes.
In recent years, considerable efforts have been invested in the development of stable and convergent numerical schemes for general second-order elliptic problems on polyhedral meshes [7][8][9][10]. In poromechanical simulations, the mechanical problem has traditionally been solved with the Finite-Element Method (FEM). To overcome the mesh restrictions imposed by the standard FEM and handle arbitrary cell shapes, generalized finite-element discretizations have been proposed for polyhedral meshes, such as the polyhedral finite-element method (see [11] and references therein). For completeness, we mention here that discontinuous Galerkin methods, hybrid high-order methods [12], and finite-volume approaches [13,14] have also been designed to solve the mechanical component of Biot's equations on polyhedral meshes. Recently, using concepts inspired by Mimetic Finite Difference (MFD) methods, the introduction of the Virtual Element Method (VEM) provided a variational framework to construct a consistent and stable scheme on arbitrary polyhedral meshes [8,15]. An attractive feature that distinguishes VEM from other generalized FEMs is that the assembly of the VEM discrete equations does not require knowledge of the analytical expression of the (non-polynomial) basis functions, which allows a generic and robust treatment of complex cell geometries. For the same reason, hanging nodes can be dealt with in a simple and convenient fashion. The scheme's properties have been studied extensively [16][17][18][19][20], and have been assessed numerically in a wide range of physical simulations, including geomechanical applications [21][22][23][24]. In this work, we take advantage of the appealing features discussed above, using a low-order VEM relying on vertex-based degrees of freedom (dofs) to approximate displacements in the discrete mechanical problem.
For the flow problem, a large number of numerical schemes for polyhedral meshes have been proposed [25], and a complete review is out of the scope of the present work. We focus here on locally conservative, low-order, linear schemes. Pioneering linear finite-volume approaches have relied on a cell-centered discretization, such as Multi-Point Flux Approximation (MPFA) schemes [26,27]. However, linear schemes combining cell-centered dofs with face-centered dofs-such as the family of Hybrid Mimetic Mixed (HMM) methods [28][29][30]-or with vertexbased dofs-such as the Vertex Approximate Gradient (VAG) scheme [31]-have been shown to be convergent on polyhedral meshes. A recent numerical comparison of some of these schemes with nonlinear methods can be found in [32], and a general analysis framework is provided by [10,33]. Motivated by the fact that VEM and mimetic schemes share a common theoretical basis, we make a natural choice-not studied before in a poromechanical context to the best of our knowledge-to couple VEM with a low-order hybrid MFD discretization that approximates pressure with cell-centered dofs and velocity with face-centered dofs. In the hybrid formulation, the MFD scheme also involves facebased Lagrange multipliers, which allows a cell-wise assembly of the discrete flow equations. Static condensation is used to eliminate the velocities during assembly and reduce the size of the linear systems.
The proposed hybrid MFD-VEM numerical scheme uses the lowest order Virtual Element space for the displacement field and a piecewise-constant interpolation for pressure. This combination of displacement and pressure approximation spaces does not intrinsically satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) inf-sup stability condition [34,35]. For incompressible problems approaching undrained conditions, this results in the development of spurious modes in the pressure field for specific mesh topologies [36]. To circumvent this issue, the scheme is stabilized using a local pressure-jump method inspired by [37,38] and adapted here to unstructured meshes. We design a fully implicit, fully coupled solution strategy for the stabilized discretization. The workhorse of the algorithm is a block-triangular preconditioner constructed specifically for the linear systems arising from the hybrid MFD-VEM scheme using a methodology presented in [38]. Using this, we demonstrate three key features of our numerical framework applied to arbitrary polygonal meshes with highly distorted cells. First, the scheme is convergent upon space-time refinement and matches the analytical solution of the well-known benchmark problems [39]. Second, the scheme exhibits an accurate behavior for incompressible problems with small time steps. We note in particular that no stabilization is needed on a large class of arbitrary polygonal meshes [40], and that local stabilization effectively damps spurious pressure modes on those meshes triggering the instability. Third, the linear solution strategy is robust and scalable for all the polygonal meshes considered here. These encouraging results will be extended to three-dimensional cases in future work. This paper is organized as follows. We review the strong and weak forms of Biot's poroelasticity equations in Section 2. In Section 3, we introduce the hybrid MFD-VEM scheme to solve the initial boundary value problem. We also introduce the local pressure-jump method used to stabilize the scheme. Section 3.6 is dedicated to the presentation of the linear solution algorithm for the coupled systems. In Section 4, we demonstrate the accuracy of the proposed numerical scheme on polygonal meshes and we illustrate the scalability of the solution algorithm. The paper ends with a few concluding remarks regarding future work.

Initial-boundary value problem in strong form
We consider a displacement-velocity-pressure formulation of Biot's poroelasticity equations [41,42] in a twodimensional domain Ω ∈ R 2 . Let I = (0, T ) denote the time interval. The three-field strong form of the initialboundary value problem (IBVP) consists of a linear momentum balance equation, a mass balance equation, and Darcy's law. Using an excess pressure formulation, the displacement u : Ω × I → R 2 , the Darcy velocity q : Ω × I → R 2 , and the excess pore pressure p : Ω × I → R satisfy: in Ω × I (momentum balance), (1a) in Ω × I (mass balance).
In the momentum balance equation, the total Cauchy stress tensor is σ(u, p) = σ (u)−αpI, where σ (u) is the effective stress tensor, α is Biot's coefficient, and I is the second-order unit tensor. In this work, we consider isotropic linear elastic materials. Hence, the effective stress can be defined as σ (u) = 2Gε(u)+λ trace(ε(u)), with ε(u) = 1 2 (∇u+∇u T ) the linearized strain tensor, and λ and G the Lamé parameters of the porous medium. The vector b(x, t) denotes body forces. In Darcy's law, κ denotes the intrinsic permeability tensor of the porous medium divided by the fluid viscosity, which is assumed constant. Using the superposed dot,(), to denote a derivative with respect to time, the fluid increment [41] is defined asζ(u, p) = α divu + S εṗ , where S ε is the constrained specific storage coefficient, i.e. the inverse of Biot's modulus.
To define the boundary conditions in the mechanical and flow problems, the domain boundary Γ is decomposed as Γ = Γ u ∪ Γ σ and Γ = Γ p ∪ Γ q such that Γ u ∩ Γ σ = ∅ and Γ p ∩ Γ q = ∅. We prescribe the following boundary conditions to, respectively, the displacement, total Cauchy stress tensor, Darcy velocity, and excess pore pressure fields: where the (space-and time-dependent) boundary values on the right-hand sides are denoted by the notation (·). In (2b), t is the prescribed traction, and n is the outward normal vector on Γ.
To complete the definition of the problem, we impose the following initial condition on the excess pore pressure field: with p 0 its initial value. Frequently, though not always, p 0 = 0. The initial displacement u 0 and velocity q 0 are then computed such that Eqs. (1a) and (1b), respectively, are satisfied.

Weak statement of the model problem
In the weak statement of the model problem, we retain a mixed structure in which mass balance and Darcy's law are kept as separate equations. This mixed formulation will be used in Section 3.2 to construct a consistent discretization of the flow problem based on a hybrid MFD scheme. The weak form of the IBVP (1) involves the following spaces: where H 1 (Ω) and H(div; Ω) are the Sobolev spaces containing, respectively, the vector functions whose first derivatives belong to L 2 (Ω), and the vector functions whose divergence is in L 2 (Ω). Using the notation (·, ·) to denote the inner product on the functional space specified in the subscript, the weak form of the mixed IBVP (1) consists in finding the triplet (u(t), q(t), p(t)) ∈ U × Q × P that, for all t ∈ I, satisfies: (α divu, χ) L 2 (Ω) + (div q, χ) L 2 (Ω) + (S εṗ , χ) L 2 (Ω) = 0 ∀χ ∈ L 2 (Ω).
A detailed analysis of this problem can be found in [43]. In the next section, we discretize the weak form (5) on arbitrary polygonal meshes by applying a low-order VEM to the momentum balance equation and a hybrid MFD scheme to the mass balance and Darcy's law.

Polygonal meshes
To describe the arbitrary polygonal meshes considered in this work, we rely on a mesh structure consisting of the triplet {T , F , V} (see also [13,44]). In this triplet, T denotes the set of non-overlapping open polygonal cells in the mesh such that Ω = ∪ K∈T K. Let ∂K = K \ K be the boundary of cell K ∈ T . The second entry in the triplet, F , is the set of mesh faces corresponding, in the two-dimensional case considered here, to one-dimensional hyperplanes of R 2 . The set of mesh faces is decomposed into three non-intersecting subsets F = F int ∪ F p ∪ F q , containing respectively the interior faces, the boundary faces located on Γ p , and the boundary faces located on Γ q . For all K ∈ T , F K is a subset that contains the faces of cell K such that ∂K = ∪ f ∈F K f . Denoting by T f = {K ∈ T : f ∈ F K } the subset consisting of the cell(s) adjacent to interface f , we require that the mesh is conforming in the sense that card(T f ) = 2 if f ∈ F int and card(T f ) = 1 otherwise. The third entry of the triplet, V, is the set of mesh vertices. The subset of the vertices of cell K is V K , and the subset of the vertices of face f is V f .
The notation | · | refers to the d-measure of a d-dimensional quantity (cell or face). We denote the positions of the center of cell K, the center of face f , and the vertex v by x K , x f , and x v , respectively (see Fig. 1). The geometrical computations in the next sections involve the outward normal to face f ∈ T K with respect to cell K, denoted by n K, f , and the vector connecting x K to x f ( f ∈ F K ), denoted by c K, f (see Fig. 1). Using these notations, a mesh satisfies the κ-orthogonality condition [25] if: We make the regularity assumption that there exists γ > 0 such that for all K ∈ T , K is star-shaped with respect to a ball of radius larger than γh K , where h K is the cell diameter. Figure 2 illustrates the different types of meshes considered in our numerical examples. They include skewed meshes (Skewed) obtained by perturbing the vertices of a uniform Cartesian mesh (Cartesian) using the function: x + 0.07 sin (4πx) cos 4πy + π 2 y + 0.07 sin (4πx) cos 4πy + π 2 The test suite also involves hybrid meshes (Hybrid) composed of triangles and quadrilaterals generated with GMSH [45], and arbitrary polygonal meshes (Polymesher1 and Polymesher20) generated with respectively one and 20 smoothing steps of Polymesher [46]. We stress the fact that the meshes of types Skewed, Hybrid, Polymesher1, and Polymesher20 do not satisfy, in general, the κ-orthogonality condition (6).

Fully discrete coupled scheme
In this section, we present the system of algebraic equations arising from a hybrid MFD-VEM discretization of the weak form (5). To simplify the presentation, we will momentarily delay the description of the construction of certain VEM and MFD operators, to first describe the overall form of the discrete equations. A subsequent section will then revisit these important but subsidiary elements in detail.
To discretize the momentum balance equation (5a), we consider a low-order VEM in which the displacement field is approximated in the following functional spaces: The local space U h,K is "virtual" in the sense that the analytical expression of its basis functions is not known and not needed for the construction of the scheme. For the low-order VEM, a function of U h,K is uniquely defined by its values at the vertices of K. These vertex-based values are, for each component of the displacement, the dofs of the VEM. The VEM methodology provides two key operators acting on the displacement dofs: a coercive bilinear form a h (·, ·) : U h × U h → R approximating the [L 2 (Ω)] 2×2 -inner product of Eq. (5a), and a discrete divergence operator discretizing the divergence term appearing in the L 2 (Ω)-inner product of Eq. (5a). We denote this VEM divergence operator by div vem h : U h → P h , where P h is: The definition of the local virtual space U h,K as well as the construction of the VEM operators a h (·, ·) and div vem h are reviewed in Section 3.3. The right-hand side of (5a) involving body forces is computed by defining, for each K ∈ T , b K = 1 |K| K b (where the integral of a vector is intended to be performed component-wise) and approximating the local right-hand side as follows: In the low-order MFD scheme, we adopt a purely discrete representation of the solution fields for the velocity, pressure, and Lagrange multiplier variables. To represent the Darcy velocity field, we denote by W h the set of one-sided face-based discrete fields. In the hybrid formulation, a discrete field w h = (w K, f ) K∈T , f ∈F K ∈ W h contains one dof per boundary face, and two dofs per interior face that are not necessarily equal. Each dof approximates the average Darcy velocity over a face such that To discretize the weak form of the Darcy equation (5b), the MFD scheme involves a discrete weighted inner product [·, ·] W h : W h × W h → R used to approximate the [L 2 (Ω)] 2 -inner product of Eq. (5b). For the discretization of the divergence term present in the L 2 (Ω)-inner product of Eq. (5b), we also define a discrete divergence operator, div mfd h : W h → P h , where the set of cell-based discrete fields, P h , is isomorphic to the space of piecewise constant functions, P h . The construction of these two MFD operators is detailed in Section 3.4.
The pressure field is represented as a cell-based discrete field of P h . A discrete pressure solution p h = (p K ) K∈T ∈ P h is a collection of dofs approximating the cell-based pressure averages, i.e., Equivalently, the pressure solution can be viewed as a piecewise-constant function of P h using the decomposition classically defined as the sum of local inner products: In the hybrid formulation, the computation of the Darcy velocity involves Lagrange multipliers that are represented as a face-based discrete field. We remark that, considering the set of face-based discrete fields, L h , the field π h = (π f ) f ∈F ∈ L h contains one dof per face that can be viewed as an approximation of the face pressure average, i.e., The subset L h,0 contains the discrete fields We consider a fully implicit (backward-Euler) temporal discretization of the coupled system. The superscript n denotes the time level at which the degrees of freedom are evaluated. Using the notations introduced above, the coupled problem in discrete weak form reads: given two functions {u 0 , p 0 } defining the initial state, find the function u n h ∈ U h and the discrete fields {w n h , p n h , π n h } ∈ W h × P h × L h such that for n ∈ {1, . . . , N}: where the right-hand side of (15c) is: In the hybridized system (15), the use of Lagrange multipliers allows a local, cell-wise computation of the one-sided face velocities in the discrete Darcy equation (15b). To ensure that the hybrid scheme remains mass conservative, the set of algebraic constraints (15d) imposes the continuity of the velocities at the mesh faces. We stress the fact that, despite a relatively large number of dofs, the hybridized system is amenable to static condensation, which is used to locally eliminate the one-sided face velocities during the assembly. The resulting algebraic system solved by the linear solver is discussed in Section 3.6. In the following sections, we focus on the terms of (15) that have not been fully defined yet. In Section 3.3, we define the virtual space U h and review the construction of the VEM operators a h (·, ·) and div vem h . In Section 3.4, we show that a similar methodology is used in MFD to form the operators [·, ·] W h and div mfd h . For simplicity, we drop the superscript n denoting the time level.

Local virtual space and VEM operators
We start the section with the local virtual space introduced in the low-order VEM [8,15] to approximate the displacement variable. Let K ∈ T be a polygon of the tessellation of Ω. Following [15], we define the following set of scaled monomials on K: where we recall that x K = (x K , y K ) is the center of K and h K is the diameter of K. We define the projection operator Π ∇,K 1 : Then, we define the following functional space: where P 1 (K) is the set of polynomials of degree ≤ 1 defined on K, and P 1 ( f ) the set of polynomials of degree ≤ 1 defined on the face f . A function η h ∈ U h,K is completely defined by its values at the vertices of K [8]. Functions in U h,K are called "virtual" because their analytical expressions is not known. Instead, we know that they are polynomials of degree ≤ 1 on each edge. For each η h ∈ U h,K , its projection Π ∇,K 1 η h is computable from the degrees of freedom. Indeed, the right-hand sides of (18) can be computed knowing only the analytical expression of η h on ∂K. In particular, the first right-hand side is computed by integrating by parts and applying Green's theorem. Exploiting the fact that ∆m = 0 ∀m ∈ P 1 (K), we get The last condition on the functions in U h,K (see Eq. (19)) allows us to compute the integral mean of η h on K. Indeed, Finally, the integral mean of the gradient of η h on K is also computable knowing only the analytical expression of η h on ∂K: . We use the same notation in the following for vectorial and tensorial functions, where integral means are performed on each component. With reference to the notation introduced in Section 3.2, we discretize the displacement u n h defining the space Furthermore, the discrete bilinear form of the problem is obtained by first defining the following discretizations of the strains and effective stresses, for each K ∈ T and each η h = η h,x η h,y ∈ U h : where that is used in (15a) and (16).
To complete the definition of VEM discrete bilinear forms, following [15], we define, for each K ∈ T , the continuous and coercive bilinear form The bilinear form S K h is defined in order to be computable from the degrees of freedom of the two functions involved and to ensure the coercivity of a K h . The most common choice, and the one we make here, is where we recall that V K denotes the set of the vertices of K. With the above definitions, we can build the discrete VEM bilinear form, that is used in (15a), as Moreover, Eq. (21) applied to each component of a test function η h ∈ U h , allows us to compute the approximation of the right-hand side term involving the body force as described in (10). Finally, we note that, since the VEM basis functions are known to be polynomials on faces, we can compute the right-hand side of (15a) involving t as we would do for a classical finite element discretization.
Remark 1. The bilinear form a K h defined by (27) is composed of two terms: the first one (first addend in (27)) accounts for the consistency of the scheme on polynomials of degree 1, the second one (second addend in (27), defined by (28)) accounts for coercivity. Indeed, whenever either η h ∈ [P 1 (K)] 2 or u h ∈ [P 1 (K)] 2 is the identity for polynomials of degree 1. Moreover, since the operators σ K and ε K are consistent on polynomials of degree 1 (see (24) and (25)), we have that a K h is exact on [P 1 (K)] 2 × [P 1 (K)] 2 . The same structure can be observed for the MFD bilinear form [· , ·] W h , that is described in Section 3.4 (see (33)).

MFD operators
In this section, we focus on the two key discrete MFD operators-namely, the divergence operator and the weighted inner product-acting on the set of one-sided face-based fields W h that appear in the discrete Darcy and mass balance equations (15b)-(15c). Let w h , q h ∈ W h be two discrete fields representing one-sided face velocities. The discrete divergence operator, div h,MFD : The weighted inner product [·, ·] W h : W h × W h → R is written as the sum of local inner products: The mimetic discretization framework provides a methodology to construct a symmetric positive definite matrix of size card(F K ) × card(F K ), denoted by M K , that only depends on the permeability and geometric properties of cell K.
As in Section 3.3, an expression for M K is obtained by imposing consecutively a consistency condition and a stability condition [7], which can then be used to demonstrate that the scheme is convergent for elliptic problems in mixed form [47]. For a low-order mimetic finite-difference scheme, the consistency condition states that the local inner product [·, ·] W h|K must be exact when one of the two arguments is the projection on W h|K of a constant function κ |K · ∇q with q ∈ P 1 (K). Setting q to q 1 (x) = x − x K and q 2 (x) = y − y K in the consistency condition results in the following algebraic constraint on the entries of M K [7]: Let f i denote the i-th face in F K . By construction, N K and R K are two full-rank matrices of size card(F K ) × 2 such that the i-th row of N K is n K, f i κ, and the i-th row of R K is | f i |c K, f i . Condition (32) defines a family of consistent loworder mimetic finite-difference schemes, but does not provide a unique expression for M K . It is worth noting that the symmetric matrix 1 |K| R K κ −1 R K satisfies (32), but is only positive semidefinite. To enforce the SPD structure of M K and obtain a coercive bilinear form, a stabilization term is introduced. Let C K be a matrix of size card(F K ) × (card(F K ) − 2) whose columns form a basis of ker(N K ), and let U K be a SPD matrix of size (card(F K ) − 2) × (card(F K ) − 2). As shown in [48], writing M K in the generic form: yields a coercive local bilinear form that still satisfies the algebraic consistency condition (32) since C K N K = 0. In (33), it is clear that the structure of the mimetic inner product is analogous to that of the VEM operator a K h as both bilinear forms can be split into a term ensuring consistency and a term enforcing stability (see Remark 1). Multiple definitions of the stabilization term in (33) have been proposed in previous work [25]. Here, we obtain an inner product satisfying the MFD stability condition by setting the scaling coefficient as: and by writing U K = (C K C K ) −1 , which yields after some simplifications [7]: Remark 2. We note that setting: yields a matrix containing the inverses of the standard TPFA half-transmissibilities divided by | f | 2 . This can be used to recover the cell-centered linear TPFA scheme after algebraic elimination of the face variables [49]. However, it is well known that the linear TPFA scheme is not consistent when the mesh does not satisfy the κ-orthogonality condition (6).

Local pressure-jump stabilization
The stability of the coupled scheme described in (15) is subject to the well-known Ladyzhenskaya-Babuška-Brezzi (LBB) inf-sup stability condition [34,35]: The proposed scheme relies on approximation spaces for the displacement and pressure variables that do not, in general, satisfy the discrete inf-sup condition given in (37). Specifically, for incompressible solid and fluid constituents, i.e. S ε = 0, in the presence of undrained conditions-resulting for instance from small time steps and/or low permeability-the instability may manifest itself by the presence of spurious modes in the pressure field (checkerboarding).
It is worth noting that, for specific mesh topologies, the spurious pressure modes do not appear and no stabilization is required. In particular, the numerical examples of Section 4.3 indicate that the Polymesher1 and Polymesher20 mesh families do not exhibit any checkerboarding. These results are in agreement with the findings of [36,40,50], in which equal-order interpolation schemes applied to the Stokes problem are shown to be stable whenever no vertex in the mesh is connected to more than three faces.
To eliminate checkerboarding in the unstable mesh configurations-in this work, the Cartesian, Skewed, Hybrid mesh types-we rely on the local pressure-jump stabilization technique introduced originally for the Stokes problem [51,52], and more specifically, we adapt the methodology used in [37,38] to unstructured meshes. To achieve this, we partition the mesh into non-overlapping cell aggregates referred to as macro-elements using the algorithm presented in Appendix A. We emphasize that these macro-elements are unstructured, in contrast to the traditional concept of structured macro-elements on a logically-nested grid. Considering the set of macro-elements E, a macro-element E ∈ E is constructed such that E = ∪ K∈E K and ∪ E∈E E = T . We denote by F E,int the set of mesh faces that are in the interior of E. Following [37,38], the local pressure-jump stabilization consists in introducing artificial fluxes at the internal faces f ∈ F E,int of each macro-element E to prevent the development of spurious pressure modes. This is done by adding the following term in the mass balance equation (15c): where the pressure jump for an internal face f ∈ F E,int is [[p n h ]] f := p n K − p n L (K, L ∈ T f with K < L). We introduce a geometric coefficient, Υ f , representing a characteristic area associated with an internal face f ∈ F E,int and defined as follows. Considering cell K ∈ E, we associate each node v ∈ V K with the measure, denoted by m K,v , of the quadrilateral whose vertices (appropriately ordered) are v, the centers of the two faces adjacent to v in K, and the centroid of K (see Fig. 3). Using that, we define Υ f as: The two-dimensional stabilization coefficient [38,53], denoted by β, is computed using Biot's coefficient and the Lamé parameters as: Since the artificial fluxes are only added at the internal faces of the macro-elements, the stabilized scheme is not cellwise mass conservative but remains mass conservative at the level of the macro-elements-in the sense that the sum of the fluxes at the external faces of a macro-element is equal to zero for the incompressible setting considered here. The stability properties of the scheme and their impact on the performance of the iterative linear solver are assessed in Section 4.3.

Solution strategy
The discrete weak form (15) produces a sequence of 4 × 4 block systems of algebraic equations of the type where • A uu is a symmetric positive definite (SPD) matrix corresponding to the elasticity block; • A ww is a block-diagonal SPD matrix, having as many blocks as the number of cells, each one with a size equal to the number of faces of the corresponding cell; • A pp consists of two contributions: (i) a diagonal matrix that depends cell-wise on the specific storage coefficient S ε , and (ii) a local stabilization contribution (Eq. (38)), if needed, whose sparsity pattern is a subset of that of a discrete (cell-centered) Laplace operator; • A up and A wp are matrices generated by the inner products (α I P h χ n h , div vem h η h ) L 2 (Ω) and [χ n h , div mfd h ϕ h ] P h , respectively; • A wπ is a matrix whose columns correspond to a unique mesh interface, having two (respectively one) negative unit entries for interfaces belonging to F int (respectively F p ∪ F q ).
Following standard practice [54], we take advantage of the block-diagonal structure of A ww by reducing (41) via static condensation to a 3 × 3 linear system with A pp = (A pp + ∆tA wp A −1 ww A wp ), A pπ = (A wp A −1 ww A wπ ), and A ππ = (A wπ A −1 ww A wπ ). We note that A pp shares the same sparsity pattern as A pp since the matrix arising from the static condensation, i.e. ∆tA wp A −1 ww A wp , is diagonal. The non-symmetric linear system (42) is solved by a preconditioned Krylov subspace method. It shares the same properties of the Biot system addressed in [38], which is obtained based on a mixed hybrid finite element formulation of the same three-field formulation considered in this work. Thus, we adopt the same block-triangular preconditioning strategy, used in conjunction with a right-preconditioned generalized minimal residual (GMRES) method [55]. Denoting by P −1 the preconditioning operator, we work with the modified system: The preconditioner in factorized form reads: with A −1 uu a suitable approximation of A −1 uu , B −1 pp a suitable approximation of the inverse of the first-level SPD Schur complement B pp = (A pp + A up A −1 uu A up ), and C −1 ππ a suitable approximation of the inverse of the second-level SPD Schur complement C ππ = (A ππ − ∆tA pπ B −1 pp A pπ ). We define B −1 pp considering the following expression: where the so-called fixed-stress assumption [56][57][58] is used to introduce a sparse approximation of the triple product A up A −1 uu A up . In (45), diagm(·) is an operator constructing a diagonal matrix by extracting the diagonal entries of the input matrix, and D uu = diagm(A uu ). Based on (45), the action of B −1 pp on a vector is always expressed by a single 1 -Jacobi iteration [59,60]. Therefore, B −1 pp is available explicitly and allows us to replace, for preconditioning purposes, the exact second-level Schur complement with The approximations considered for A −1 uu and C −1 ππ are provided in Section 4. For a detailed analysis of the features of preconditioner (44) the reader is referred to [38].

Numerical examples
We now consider three sets of numerical experiments. In the first set, Mandel's problem (Fig. 4a), a classical benchmark of linear poroelasticity for which an analytical solution is available [39], is used to validate our MFD-VEM formulation. In the second set, we numerically investigate the scheme's convergence properties based on a manufactured regular solution [61], solving the IBVP (1) on the unit square domain Ω = [0, 1] 2 . In the third set, a cantilevered square block (Fig. 4b) is considered to demonstrate the robustness of the proposed stabilization with respect to pressure oscillations in the incompressible limit. In particular, we emphasize its beneficial effects on the iterative solver convergence. Details of the material and parameter values used in the three test cases are summarized in Table 1.
In all tests, we assume the zero vector as initial guess for (non-restarted) GMRES and terminate the iterations when the initial residual has been reduced by a factor of 10 6 . As to the preconditioner (44), the operatorÃ −1 uu is defined based on the so-called separate displacement component part of the stiffness matrix [62,63], i.e. a sparse approximation to A uu in which x-and y-displacement dofs are decoupled. We consider either a sparse direct solver or algebraic multigrid (AMG) preconditioning for bothÃ −1 uu andC −1 ππ . Specifically, we use a classic AMG method [64] as provided by the HSL MI20 package [65] with default parameters (symmetric Gauss-Seidel smoother, single V-cycle, and a direct coarse solver) except for the coarsening failure criterion control (c fail), which is set to 2. An 1 -Jacobi [59,60] smoother is always employed for B −1 pp .

Mandel's problem
Mandel's problem [39] consists of a sample of saturated, isotropic, poroelastic material that is loaded under planestrain conditions by a constant compressive force of magnitude 2F applied at time t = 0 on rigid, frictionless, impermeable plates. The sample has dimensions 2a × 2b (Fig. 4a). The left and right sides (x = ±a) are stress free, drained and kept at constant ambient pressure. The exact analytical solutions of this problem are known (see [66,Appendix A.2]).
Given the symmetry of the problem, we solve it on a quarter of the domain, represented by the square Ω = (0, a) × (0, b) (Fig. 4a). The simulation setup is the same as that used in [66]. We use the unit square as the domain, i.e. a = b = 1 m. The rigid plate constraint is accounted for by prescribing the vertical displacement at the loaded boundaries using the closed-form solution of the problem.
In Fig. 5 we compare analytical and MFD-VEM solutions for the pressure along a cross-section at y = 0.5 at different times. The results are relative to 400-cell meshes from the mesh families in Fig. 2. A constant time step size ∆t = 10 −4 · T c is used in both cases, with T c = a 2 /(κ(λ + 2G)) the characteristic time of the consolidation process. We can see that the computed values match the expected profile very well, with relative error values of the order of 10 −4 . Note that the parameters used are representative of the limit case of incompressible solid (b = 1.0) and fluid (S ε = 0) constituents, which maximizes the hydromechanical coupling, i.e. the so-called Mandel-Cryer effect. These results confirm the robustness of the method with respect to badly shaped elements. Because of the lack of regularity of the pressure field for Mandel's problem [67,68], we do not perform a mesh refinement study to evaluate numerically the convergence rate of quantities of interest here. This is addressed in the next section.
The following error measures are considered: The above quantities are computed on the families of meshes illustrated in Fig. 2. These meshes are refined together with the time discretization parameter ∆t, in such a way that, if (N cells k , ∆t k ) is the pair of discretization parameters at the k-th refinement (number of cells and time discretization parameter, respectively), then (N cells k+1 , ∆t k+1 ) = (4N cells k , ∆t k /2). In Fig. 6 we display the behavior of the above quantities with respect to the maximum diameter of the discretization h, at each refinement step, reporting the experimental convergence rates in the legend, computed using the last two points of each line. We can see that the rates of convergence are in very good accordance with the theory on the more regular meshes (Polymesher20 and Cartesian), while on less regular meshes we still see a preasymptotic behavior, due to the fact that h does not scale exactly as 1 √ N cells . To check the correct convergence rates in time, in Figs. 7 and 8 we display the behavior of e p and e u for fixed h, as functions of the time discretization parameter ∆t. We can see that, once the time component of the error becomes dominant, we obtain the expected rates of convergence. Finally, Figs. 13 to 10 display the errors computed for each choice of ∆t and h and confirm the behaviors expected from theoretical results. Indeed, it can be checked that all quantities are asymptotically decreasing linearly with respect to h; e p and e u decrease linearly with respect to ∆t, and e σ is asymptotically constant with respect to ∆t.

Cantilevered square block problem
We consider the cantilever problem discussed in [69] and later used in [38] to demonstrate the robustness and computational efficiency of coupled schemes relying on the unstructured macro-element stabilization. A sketch of the setup of the problem is shown in Fig. 4b. We assume the problem domain Ω to be the unit square. The mechanical boundary conditions fix the displacement (u = 0) on the left, and impose a unit downward traction at the top and a zero traction on the right and at the bottom. We impose no-flow boundary conditions on the four sides. The problem parameters are identical to those used in [38] and are provided in Table 1. Sequences of refined meshes for each family are employed for a total of 6 levels ( Table 2), with a problem size ranging three orders of magnitude in terms of global number of unknowns.
First, we focus on the effectiveness of the proposed stabilization technique in preventing the occurrence of oscillations in the discrete pressure field. A contour plot of the pressure MFD-VEM solution on the level 0 mesh of each                 Table 2). In panels a, b, and c, the star symbol ( ) superscript indicates that the local pressure-jump stabilization was introduced, with the unstructured macro-element mesh highlighted using thicker edges. The remaining panels (d to h) were obtained with the unstabilized formulation. family after a single timestep ∆t = 1 × 10 −5 s is given in Fig. 14. The use of an unstabilized formulation produces checkerboard oscillations for the Cartesian, Skewed, and Hybrid meshes, which are evident in the bottom panels ( Fig. 14f, g, h). Such oscillations are successfully removed by the pressure-jump stabilization (Fig. 14a, b, c). As expected, stable results are obtained with meshes Polymesher1 and Polymesher20 without stabilization (Fig. 14d, e) since all the vertices are connected to at most three faces [36]. Conversely, for larger timestep sizes the MFD-VEM formulation becomes intrinsically stable with no stabilization required (Fig. 15). Mesh level
Finally, the algorithmic weak scalability of the iterative linear solver strategy described in Section 3.6 is numerically investigated. We solve (42) only once for two different timestep sizes. Figs. 16-17 display the number of GMRES iterations needed in the sequence of refined problems detailed in Table 2 When nested direct solvers are used to approximate the action of A −1 uu and C −1 ππ (Fig. 16), two behaviors may be observed in case of a small timestep relative to the characteristic consolidation time, ∆t = 1 × 10 −5 s. For a stable formulation-stabilized (Cartesian ( ) , Skewed ( ) , Hybrid ( ) ) or intrinsically stable ((Polymesher1, Polymesher20)-the preconditioner P −1 performs very well, with only a modest increase in the iteration count between levels 0 and 5. Such almost-optimal behavior of P −1 with respect to mesh size is lost in case of an unstable formulation. Indeed, GMRES shows erratic iteration counts that reflect the sensitivity of Krylov-based solvers to the presence of near-singular pressure modes. For larger timestep size, i.e. ∆t = 1 × 10 −1 s, the system is far from the incompressibility limit (Fig. 17). Hence, the linear solver always exhibits robust convergence and the stabilization effects become irrelevant. Replacing nested direct solvers by AMG preconditioning for A −1 uu and C −1 ππ preserves the desired behavior at a substantially smaller cost.

Closure
In this work, we have presented a numerical scheme coupling hybrid MFD and VEM to discretize Biot's equations of poroelasticity on arbitrary polygonal meshes. A key feature of the discretization is that it remains convergent in the presence of highly distorted cells with arbitrary shapes, as demonstrated numerically in Section 4. For incompressible problems approaching undrained conditions, the discretization is stabilized with a local pressure-jump technique applicable to unstructured meshes to prevent the development of spurious pressure modes (checkerboarding). The proposed simulation framework also includes a fully coupled linear solution strategy, based on a block-triangular preconditioner, with excellent scalability.
Future work will focus on the extension of these results to three-dimensional polyhedral meshes to further demonstrate the potential of the numerical framework and confirm scalability in large-scale problems. A key goal is to address both fully unstructured meshes and stratigraphic corner-point grids [25]. The latter remain the industry-standard approach in reservoir engineering studies, but contain deformed (and sometimes degenerate) hexahedral cells with non-matching faces that cannot be handled by the standard FEM.
Algorithm 1 Unstructured macro-element construction algorithm 1: Mark all vertices as not visited yet. 2: for v in V int do 3: if v is marked not visited yet then 4: Form a macro-element with the cells adjacent to v.

5:
Mark all the vertices of the cells forming the new macro-element as visited. 6: end if 7: end for 8: Collect cells not yet assigned to a macro-element in list C. 9: while C is not empty do 10: Pick cell K at the top of list C.

11:
if K is adjacent (through faces) to at least one macro-element then 12: Assign K to the neighboring macro-element with the smallest number of cells. 13: Remove K from list C.
14: else 15: Put K at the back of list C. 16: end if 17: end while Appendix A. Unstructured macro-element construction for local pressure-jump stabilization Here, we review the partitioning of the mesh involved in the local pressure-jump stabilization. This is done using the mesh connectivity in two steps described in Algorithm 1. We obtain non-overlapping macro-elements containing at least three cells.
To initialize the algorithm, we mark all the vertices as not visited yet. In the first step, we loop over the internal vertices of the mesh. If a vertex is marked as visited, it is skipped and we proceed to the next vertex in the list. If not, we form a macro-element made of the cells adjacent to this vertex, and we mark this vertex as visited, as well as all the vertices of the cells forming the macro-element. At the end of this first step, some cells are, in general, still unassigned to a macro-element.
In the second step, we collect these unassigned cells in a list, and for each cell in the list, we proceed as follows. If an unassigned cell is adjacent (through faces) to at least one macro-element, this cell is assigned to its neighboring macro-element with the smallest number of cells, and is then removed from the list. If an unassigned cell is surrounded by unassigned cells, it is placed at the back of the list, and will be processed later. This procedure is applied until the list is empty, at which point all cells are assigned to a macro-element.