Virtual Element Approximations for Two Species Model of the Advection-Diffusion-Reaction in Poroelastic Media

. This paper proposes virtual element methods for approximating the mathematical model consisting of coupled poroelastic and Advection-Diffusion-Reac-tion (ADR) equations. The space discretization relies on virtual element spaces containing piecewise linear polynomials as well as non-polynomials for displacement, pressure and concentrations, and piecewise constant for total pressure; a backward-Euler scheme is employed for the approximation of time derivative. Using standard techniques of explicit schemes, the well-posedness of the resultant fully discrete scheme is derived. Moreover, under certain regularity assumptions on the mesh, optimal a-priori error estimates are established by introducing suitable projection operators. Several numerical experiments are presented to validate the theoretical convergence rate and exhibit the proposed scheme’s performance.


Introduction
Considering the computational efficiency of polygonal meshes, in the last decade, there have been many developments in virtual element methods (VEMs)numerical techniques evolved from mimetic finite difference methods, see [1,5] and references therein.Contrary to finite element methods, in VEM, the mass and stiffness matrices are computed (with the help of degrees of freedom) without explicit knowledge of the basis functions.This is desirable while dealing with polygonal meshes and demanding more accurate solutions, i.e., higherorder approximations.Virtual element spaces involved in the discrete formulations contain the polynomial functions and non-polynomial functions (solution of suitable PDEs).However, neither the implementation nor the theory demands the expression of these unknown non-polynomial functions.The presence of the polynomial functions in the virtual element spaces would help in demonstrating the convergence rates of the proposed VEM.We remark that convergence analysis of VEM can be carried out similarly as in finite element methods by introducing certain projection operators onto the space of polynomial functions.Recently, VEM has been also developed in [12,14,24,25,26] for the Biot's equation and for nonlinear problems in [9,13].
In this paper, we propose a virtual element method for a three-field formulation of the time-dependent poromechanics equations coupled with advectiondiffusion-reaction equations for two species.We base the development of current work as the formulation proposed in [20] and [22] for the stationary Biot system and extend the discrete analysis to include the quasi-steady case, and the development in time-dependent problems [9,27].Several other works on the Biot's equation with various numerical techniques such as finite element method, hybrid high order method, discontinuous Galerkin method are seen in [10,11,19,23,30] and references within.

Governing equations: Coupled poroelastic and a system of advectiondiffusion-reaction equations
Let us consider the model describing the fluid flow of two chemical species in poroelastic media on a domain Ω ⊂ R 2 .The coupled problem is stated as follows [15,28]: For given body load b and a mass source ℓ, one seeks for each time t ∈ (0, t final ], the displacements of the porous skeleton u s , the pore pressure of the fluid p f , the volumetric part of the total stress, or total pressure ψ, concentration of first species w 1 and concentration of second species w 2 such that − div(2µε(u s ) − ψI) = ρb + div σ act in Ω × (0, t final ], ψ − αp f + λ div u s = 0 in Ω × (0, t final ], here κ(x) is the hydraulic conductivity of the porous medium (possibly anisotropic), ρ is the density of the solid material, η is the constant viscosity of the interstitial fluid, c 0 is the constrained specific storage coefficient, α is the Biot-Willis consolidation parameter, and µ, λ are the shear and dilation moduli associated with the constitutive law of the solid structure, D 1 , D 2 are positive definite diffusion matrices.Nevertheless, for sake of fixing ideas and to specify the coupling effects also through a stability analysis, that is conducted in [15], they chose the coupling terms f and g as 1 w 2 ) + γ w 2 div u s , where β 1 , β 2 , β 3 , γ are positive model constants.The active component of total Cauchy stress, or active stress is given as σ act = −τ r(w 1 , w 2 )k⊗k, that is, σ act operates primarily on a given constant direction k, and its intensity depends on a scalar field r(w 1 , w 2 ) and on a positive constant τ .
We endow the problem (1.1) with appropriate initial data for i = 1, 2 and boundary conditions in the following manner [2µε(u s ) − ψ I + σ act ]n = 0 and p f = 0 on Σ × (0, t final ], where the boundary ∂Ω = Γ ∪ Σ is disjointly split into Γ and Σ where we prescribe clamped boundaries and zero fluid normal fluxes; and zero (total) traction together with constant fluid pressure, respectively.Moreover, zero concentrations normal fluxes are prescribed on ∂Ω.We point out that, if we would like to start with a model in terms of the divergence (div(w i u s ) instead of u s • ∇w i in (1.1), i ∈ {1, 2}), we need to assume zero total flux (including the advective term, see, e.g., [2]).Homogeneity of the boundary conditions is only assumed to simplify the exposition of the subsequent analysis.This paper is structured as follows.In Section 2, we propose the weak formulation by introducing suitable Sobolev spaces.The fully discretized scheme is presented and proved that the scheme is well-posed in Section 3. We established a-priori error estimates in Section 4 with the help of Stokes and elliptic projection operators.In Section 5, we implement the numerical examples to verify the convergence results.Finally, we conclude the work by providing the future directions in Section 6.

Weak formulation
We will use the following notations for the Sobolev spaces in this article.
Let us multiply (1.1) by adequate test functions and integrate by parts (in space) whenever appropriate.Incorporating the boundary conditions (1.3), the variational problem formulated as: For a given t > 0 and initial conditions (1.2), find where the formulation with the bilinear forms a 1 : and linear functionals F b,r : V → R (for r known), G ℓ : Q → R, J f , J g : W → R (for known f and known g), are defined as where Remark 1.Throughout this paper, z stands for the coupling terms f and g .
We will consider that the initial data (1.2) are non-negative and regular enough.Moreover, throughout the text we will assume that the anisotropic permeability κ(x) and the diffusion matrices D 1 (x), D 2 (x) are uniformly bounded and positive definite in Ω.The latter means that, there exist positive constants κ 1 , κ 2 , and Also, for a fixed u s , the reaction kinetics f (w 1 , w 2 , •), g(w 1 , w 2 , •) satisfy the growth conditions, that is, for given w 1 , w 2 ∈ R, the scalar field r(w 1 , w 2 ) and reaction kinetics f Math.Model.Anal., 27(4):668-690, 2022.
In addition, the terms in (2.1) fulfill the continuity bounds, coercivity and positivity bounds, and C k,1 and C k,2 are the positive constants satisfying and c p is the Poincaré constant.Moreover, the bilinear form b 1 satisfies the inf-sup condition (see, e.g., [17]).

Discrete formulations and wellposedness
In this section, by following virtual element methods for space and Euler backward scheme for time discretizations, we present a fully discrete scheme corresponding to (2.1).We also address the stability, existence and uniqueness of solution for the discrete problem.

Virtual element discretizations
Let the domain Ω be discretised into the family of the polygonal meshes T h with mesh size h and elements K, vertices on each element K as K number of vertices in K, and any edge in the polygonal mesh is denoted by e.For any natural number k, let P k (S) and M k (S) represent the space of polynomials and monomials, or scaled polynomials of degree less than or equal to k for any S ⊂ R 2 .We will write a ≲ b for the expressions a ≤ Cb, where C is a generic constant independent of parameters h and ∆t.
We also suppose that the polygonal mesh satisfy the assumptions (refer [5,8]).
Before proceeding to define the virtual element spaces, we will take few notations as, for any smooth enough u, v, Note that the above definition on the H 1 scalar product is defined up to a constant, and in order to determine the constant, we define another projection P 0 K as follows The vectorial energy projection from the vector space [H 1 (K)] 2 to [P 1 (K)] 2 is denoted as Π ∇ K .A variant of the vectorial projection Π ∇ K and supported by the bilinear form a and the function p 1 ∈ ker(a K 1 (•, •)) are again taken from Similar to the energy projections, the projections Π 0,0 2 are identified as the vectorial L 2 projection onto constants and linear polynomials, respectively.We stress that these operators not only help us in deriving the optimal error estimates but are also useful in the computation of discrete bilinear forms (to be defined later).
Then the local VE spaces are introduced as follows [1,4]: where the boundary space is given by The approximation space for displacement is V h , for total pressure, Z h , and for pressure and concentration is Q h (K).We have the following degrees of freedom depending on the corresponding spaces (refer [1,4,29] for details on unisolvance) as, for displacement: the values at all vertices of the element K, and value of v s • n e K at mid point on each edge e ∈ ∂K; for total volumetric stress: value of ψ at any point in K; and for pressure, or concentration: the values of q f at vertices of the element K.
Then any global space X h corresponding to space X (choices are X=V, Q, W ) is given as Note from above that the local discrete spaces and its degrees of freedom for pressure and concentration are the same, and thus have the same local projection operators on each element K ∈ T h .
We define the local discrete bilinear forms by considering the computability, consistency and stability, as: where κ and Di 's denoting the average values of the respective parameters, and the stabilization terms on each K, with Ndof denoting the dimension of the associated space to the variables (for instance, we have Ndof as dimension of Q h (K) for the variable p h , q h ), as Then the global discrete bilinear forms and discrete functional are defined naturally for any discrete form The discrete functionals are defined in terms of L 2 projections as The stabilization terms satisfies the stability condition with respect to the norm associated with the respective bilinear forms, that is, for where c * , c * , c * , c * , ĉ * , ĉ * are constants independent of h K .With the help of the stability of the projection operators and stabilizations terms, we obtain the following lemma ( for details, see [6,9,12] and references within).
Moreover, the coercivity properties can be obtained through similar arguments, which implies that for all Also, the discrete inf-sup condition hold on (V h , Z h ) (refer [7,8,29]).By introducing the adequate discrete spaces associated with velocity, pressure, total volumetric stress, and concentrations, the fully discrete formulation is described in the next subsection.

Fully discrete scheme
Let us discretise the time interval (0, t final ] into N equispaced points and time step ∆t with denoting n th time step as t n = n∆t and n = 1, . . ., N , and use the following general notation for the first order backward difference ∆t δ t X n := X n − X n−1 .In this way, we can write a discrete form of (2.1): From the given initial data u s,0 h , p f,0 h , ψ 0 h , w 0 1,h , w 0 2,h (which will be projections of the continuous initial conditions of each field) and starting with n = 1, we first solve for Then we seek the concentrations w n 1,h , w n 2,h ∈ W h for given displacement u s,n h ∈ V h (solution of (3.2)) and known initial data w 0 1,h , w 0 2,h such that ) The above process of solving the problem continues iteratively for n = 2, . . ., N .We note that the above systems of equations are linear for each n since we have considered the explicit scheme in time discretization.

Existence and Uniqueness
We will show the well-posedness of the discrete scheme through stability and then uniqueness of the linear system of equations.We start it by introducing the discrete-in-time l 2 − norm as follows, Now, we collect the following important results for the further analysis: ) The next lemma is referred from our previous work [12] and it is needed for the analysis.
Lemma 2. We have the following bound, for all q f,n h , where Next, we recall a well known lemma utilized to handle analysis of the nonlinear and time dependent problems.Also, the proof of the following lemma can be consulted from [18,Lemma 5.1].
Lemma 3 [Generalised Discrete Gronwall lemma].Let k, B, and a j , b j , c j , γ j for integer j ≥ 0 be non-negative numbers such that, for n > 0, we have Proof.The linear problem (3.2) in the form of an uncoupled fully discrete scheme for poroelasticity problem is well-posed for a given data and can be referred from [12].We will ensure the existence of a unique solution of linear uncoupled ADR equation (3.3) by virtue of the Lax-Milgram lemma.In order to proceed, we define the bilinear forms for each i = 1, 2 and given u s,n h as solution of the problem (3.2a)-(3.2c)as, We can rewrite the uncoupled ADR problem (3.3a)-(3.3b)for all s h ∈ W h as and i = 1, 2. The continuity of right hand side obtained using the bounds of the linear functionals m h (w n−1 1,h , •) and J h,n z (w n−1 1,h , w n−1 2,h , u s,n h ; •).Now, we will prove the coercivity of bilinear form C h i (•, •).For all s h ∈ W h , the usage of Inverse inequality for polynomials leads to Now, for any s h ∈ W h , the use of above bound for c K h (u s,n h ; •, •), coercivity of bilinear forms m h (•, •) and a h 3+i (•, •), and Young's inequality, we get Math.Model.Anal., 27(4):668-690, 2022.
Choosing ∆t > 0 small enough and require D a i so that α * ≥ we achieve the coercivity of C h i (•, •).Also, the continuity of the bilinear form is followed from the boundedness of the discrete bilinear forms (3.1).Thus, the Lax-Milgram lemma deduces the existence of a unique solution.⊓ ⊔ Theorem 2 [Stability of the fully discrete scheme].Assume that the solution 3) satisfies, for constant C > 0 (independent of h and ∆t), ), and for time step n and n − 1 in Equation (3.2c) with ϕ h = ψ n h , and adding these equations gives Summing over n, and a use of (3.4) and bound of L n h from Lemma (2), we get Here, the term N is defined as, Use of Young's inequality gives Using inf-sup condition and Equation (3.2a), we obtain Summing over n gives Again use of these arguments for Equation (3.3b) with s h = w n 2,h .Then adding the resultant bounds and assuming M := max < ∞, then choosing ϵ, ∆t so that, we have ∆t(ϵ + C 1 M ) < 1.Thus, the use of Lemma 3 completes the stability of the discrete solution.⊓ ⊔

Error analysis
We will assume the following regularity assumptions to generate the convergence results: For all t > 0, the displacement of porous medium u s (t) ∈ [H 2 (Ω)] 2 , the fluid pressure p f (t) ∈ H 2 (Ω), the total pressure ψ(t) ∈ H 1 (Ω) and concentrations w 1 , w 2 ∈ H 2 (Ω).We will assume the regularity in time as, Lemma 4. For each u s ∈ V ∩ [H 1+r (Ω)] 2 with 0 ≤ r ≤ 1 and under the regularity assumption on the polygonal mesh (mentioned in Section 3.1), there exist an interpolant u s I ∈ V h satisfying where (A u h , A ψ h ) and A p h , A w1 h , A w2 h are standard Stokes and elliptic projection operators respectively.

Now, we define the following projection operator
For these operators satisfy the following estimates (see, for instance, [27,29]): ) ∈ V h be the unique solution to the system of equations (2.1) and projection operators, respectively.Then, To derive the theoretical error estimates for fully discrete scheme, we decompose the error using the projection operator A h as follows : and for each n = 1, . . ., N .
From the weak formulation (2.1) and fully discrete scheme (3.2)-(3.3),and an appeal to projection operator A h , the error equations for the fully discrete scheme in terms of η n ξ , where ξ = u, p, ψ, w 1 , w 2 , are for i = 1, z = f and i = 2, z = g.Now, we will divide the derivation of the error estimates of the fully discrete scheme (3.2)-(3.3)into the two lemmas: one containing the error bounds from the poroelasticity equations and other-regarding ADR equations.We start here by recalling/mentioning the bounds to be used in succeeding lemmas.
The use of Taylor's expansion shows the following corollary holds, Corollary 1.For any ξ, we have Now, we derive the result below by using the Lipshitz condition of function r(w n 1 , w n 2 ) for each n, Corollary 1 and standard inequalities.
the solution to the system (3.2) for each n = 1, . . ., N .Then the following estimate holds, with constant C independent of h and ∆t, where ϵ 3 > 0 chosen in subsequent analysis.
Proof.Set v s h = δ t η n u in (4.2a), q f h = η n p in (4.2b), and ϕ h = η n ψ in (4.2c) then summing over each n the resultant equation with the use of Lemma 2 gives The error term R 1 can be written as The Young's inequality with some constant ϵ 1 > 0 gives The bound of R b 1 is through the property of function r(w 1 , w 2 ), triangle's inequality, and bound (4.3a) yields Use of polynomial approximation and consistency of ãh,K 2 (•, •), using Cauchy Schwarz, Poincaré and Young's inequalities, estimate (4.3b) and by Taylor's expansion, we get An application of Lemma 1 and Young's inequality implies Use of the estimates (4.1a) and Lemma 1 as seen in bound of R 2 with constant ϵ 3 > 0 gives Combining the bounds of R i 's, we get Using inf-sup condition of b 1 (•, •) and error equation (4.2a), we obtain Therefore, the adequate choices of ϵ's and choosing the initial conditions using projections as u s,0 h := u s I (0) and p f,0 h := p f I (0) conclude the error bounds (4.4).⊓ ⊔ Note that the initial conditions are chosen so that the estimates of η n u , η n p and η n ψ are known, and can consider another such choices for analysis and computations.Next, we approach the remaining error equations to avail the following lemma.
Use of estimates for projection Π 0 K , triangle's inequality and projection (4.1) gives The estimate (4.3a), and Young's inequality with ϵ 1 ϵ 4 = C 2 /4 and ϵ 1 , ϵ 4 > 0 gives The use of consistency for bilinear form m h (•, •) and bounds (4.3a)-(4.3b)with Young's inequality for ϵ 4 > 0 gives Assume that ∇w j 1 and u s,j h are bounded for each j then the values of ∇w j (by use of inverse estimate) are finite respectively.
Thus, by applying the Cauchy-Schwarz and Young's inequalities implies Thus, the bounds of A i 's gives Similar to the above bounds, taking s h = η n w2 in (4.2d) for i = 2, we get the error bounds in the terms of η w2 .Thus, the addition of the bounds for w 1 and w 2 concludes the proof.⊓ ⊔ Application of the discrete Gronwall's inequality (Lemma 3.5) with the combination of results from Lemmas 6-7 yield the final result: ,h ∈ W h be the solution to fully discrete problem (3.2)-(3.3)for each n = 1, . . ., N .Then the following estimate holds, with constant C independent of h and ∆t,

Numerical experiments
We define the L 2 and H 1 errors for the approximation space of order k as , where v can be u s , p f , ψ and w 1 , w 2 .The convergence rates of the errors E k (v) and E ′ k (v) with k = 0, 1 for the corresponding mesh sizes h and h ′ respectively, are calculated as r k

Space and time convergence
We initiate the tests in the domain Ω 1 := (0, 1) × (0, 1) and verify the spatial convergence rate of the VEM for given exact solutions by discretizing the domain into elements containing non-convex polygons seen in Figures 1-3., and scalar function r(w 1 , w 2 ) := w 1 + w 2 .Also, the exact concentration solutions are given as w 1 (x, y, t) = w 2 (x, y, t) = t sin(πx) sin(πy), and the reaction kinematics with unit value of D 1 , D 2 , β 1 , β 2 , β 3 and γ = 0.1 making the concentration equations with different load functions.However, the load functions ℓ and b, and the exact global pressure ψ is obtained from the problem (1.1) in domain Ω with τ = 1, ρ = 1.We see that the trend of the computed errors with h = ∆t in the Figure 4 is the one expected from the theoretical results in Section 4. We can note that we have considered the parametric values to check the extend of the method with small c 0 and D 1 , and large values of λ.We show the computed rate of convergence with ∆t = 0.005 and varying mesh sizes h on a distorted triangular meshes (shown in Figure 2) in the Table 1.We have computed the relative errors Ēk for any variable v as Ēk (v) = E k (v)/ ∥v∥ k,Ω and r k (v) as their rate of convergence.

Conclusions
By extending the analysis of [12,29], we have discussed and analyzed the lowest order conforming virtual element methods for the approximation of coupled poroelasticity and ADR equations.The major contributions of this article are: showing the well-posedness of fully discrete schemes and establishing the optimal a-priori error estimates for all the variables that naturally appeared in the weak formulation.The possible extensions of this work include the study of general flow-transport problems, and the coupling with other phenomena such as diffusion of solutes in poroelastic structures [28], interface elasticityporoelasticity problems [3,16], multilayer poromechanics, or multiple-network consolidation models [21].

For
given u s,n h as solution of the problem (3.2a)-(3.2c),and taking s h = w n 1,h in (3.3a) then the use of Young's inequality with appropriate choice of ϵ and the bound of trilinear form c h (•; •, •) implies Taking s h = η n w1 in (4.2d) for i = 1, then multiplying with ∆t and summing over n enable us to get (

Figure 4 .
Figure 4. Computed errors for meshes N h , H h and D h in first three figures, and last one for variable κ = x + y + 10 on mesh H h .
Lemma 7 [Coupled CDR error bounds].Let (w n 1 , w n 2 ) ∈ [W ] 2 be the solution to the continuous problem (2.2)-(2.3)and (w n 1,h , w n 2,h ) ∈ [W h ] 2 be the solution to the discrete problem (3.3) for each n.Then the following estimate holds, with constant C independent of h and ∆t,