A priori and a posteriori error bounds for the fully mixed FEM formulation of poroelasticity with stress-dependent permeability
Arbaz Khan Bishnu P. Lamichhane Ricardo Ruiz-Baier Segundo Villa-Fuentes

Tensor network square root Kalman filter 
 for online Gaussian process regression
    [
    Received 16 July 2024; accepted 04 September 2024
=================================================================================



§ ABSTRACT 
We develop a family of mixed finite element methods for a model of nonlinear poroelasticity where, thanks to a rewriting of the constitutive equations, the permeability depends on the total poroelastic stress and on the fluid pressure and therefore we can use the Hellinger–Reissner principle with weakly imposed stress symmetry for Biot's equations. The problem is adequately structured into a coupled system consisting of one saddle-point formulation, one linearised perturbed saddle-point formulation, and two off-diagonal perturbations. This system's unique solvability requires assumptions on regularity and Lipschitz continuity of the inverse permeability, and the analysis follows  fixed-point arguments and the Babuška–Brezzi theory. The discrete problem is shown uniquely solvable by applying similar fixed-point and saddle-point techniques as for the continuous case. The  method is based on the classical PEERS_k elements, it is exactly momentum and mass conservative, and it is robust with respect to the nearly incompressible as well as vanishing storativity limits. We derive a priori error estimates, we also propose fully computable residual-based a posteriori error indicators, and show that they are reliable and efficient with respect to the natural norms, and robust in the limit of near incompressibility. These a posteriori error estimates are used to drive adaptive mesh refinement. The theoretical analysis is supported and  illustrated by several numerical examples in 2D and 3D. 


Mixed finite elements,  
stress-based formulation, nonlinear poroelasticity, fixed-point operators, error estimates.
65N30, 65N15, 65J15, 76S05, 35Q74.





§ INTRODUCTION


Nonlinear interaction between flow and the mechanical response of saturated porous media is of a great importance in many applications in biophysics, geomechanics, and tissue engineering, for example. One of such models is the equations of nonlinear poroelasticity, whose mathematical properties were studied in great detail, for example, in the references <cit.>. In these works, it becomes clear that a distinctive property of nonlinear poroelasticity models targeted for, e.g., soft tissue (cartilage, trabecular meshwork, brain matter, etc.), is that the nonlinear permeability (the hydraulic conductivity,  defined as how easily pore fluid escapes from the compacted pore spaces) that depends on the evolving total amount of fluid, does not entail a monotone operator, and therefore one cannot readily apply typical tools from monotone saddle-point problems. 

Our interest is in deriving mixed finite element (FE) formulations (solving also for other variables of interest), and for this we can cite in particular <cit.>, where  fully mixed formulations based on the Hu–Washizu principle are studied. Writing the poroelasticity equations in terms of the strain tensor was motivated in particular in <cit.> because the permeability – at least in the regime we focus here – depends nonlinearly on the total amount of fluid, which is a function of strain. 

The upshot here compared to <cit.> is that we are able to rewrite the constitutive equation for permeability to depend on the total poroelastic stress and on the fluid pressure (similarly as in, e.g., <cit.>). This allows us to revert to the more popular Hellinger–Reissner type of mixed formulations for poroelasticity  <cit.> (without solving explicitly for the strain). Consequently, another appealing advantage with respect to the formulation in <cit.> is that, as in the Hellinger–Reissner formulation, the model becomes robust with respect to the Lamé constants. Also in contrast to <cit.>, in this work we use a mixed form for the fluid flow (adding the discharge flux as additional unknown), which gives the additional advantage of mass conservativity. 

Regarding the well-posedness analysis, the aforementioned non-monotonicity of the permeability suggest, for example, to use a fixed-point argument. We opt for freezing the arguments  of permeability, turning the double saddle-point structure with three perturbations coming from the stress trace operator and from the L^2 pressure blocks, into two decoupled saddle-point problems whose separate solvability can be established from the classical literature for weakly symmetric elasticity and mixed reaction-diffusion equations. Banach fixed-point theorem is then used to show well-posedness of the overall problem. This analysis needs to verify conditions of ball-mapping and contraction of the fixed-point map, and this imposes a  small data assumption, which can be carried over to the external load, mass source, boundary displacement, and boundary fluid pressure. Compared to <cit.>, these conditions are less restrictive and imply also a less restrictive discrete analysis (which follows closely the continuous one), due to the analysis being performed using the inverse of the Hooke tensor, which allows us to achieve robustness with respect to the first Lamé parameter λ. 

Note that at the discrete level we can simply use conforming FE spaces. Discrete inf-sup conditions are already well-known for the chosen FE families of PEERS_k and Raviart–Thomas elements used for the solid and fluid sub-problems (but several other inf-sup stable spaces that satisfy a discrete kernel characterisation are also possible). We emphasise that, similarly to <cit.>, all estimates hold uniformly in the limit of nearly incompressibility (implying that the formulation is Poisson locking-free) as well as when the constrained storage coefficient vanishes (poroelastic locking-free), and therefore they are free of non-physical pressure oscillations. 

An additional goal of this work is to derive efficient and reliable residual a posteriori error estimators for the nonlinear poroelasticity equations. The approach follows a similar treatment as that of <cit.> (which focuses on mixed formulations of stress-assisted diffusion equations), with the difference that here we do not need to include augmentation terms for the mixed form of the mixed diffusion problem. The main ingredients in the analysis of these estimates are Helmholtz decompositions and a global inf-sup (together with boundedness and Lipschitz continuity of coupling terms), local inverse and trace estimates, bubble-based localisation arguments, and properties of Clément and Raviart–Thomas interpolators. See also <cit.> for estimators in a similar multiphysics context and, e.g.,  <cit.> for mixed linear elasticity.  Note that for the reliability of the estimator the aforementioned Helmholtz decompositions – for both tensor-vector and vector-scalar cases – should be valid for mixed boundary conditions. For this we follow <cit.> and <cit.>, from which we inherit an a convexity assumption on the Neumann sub-boundary (where we impose traction and flux boundary conditions). 

Outline. The rest of the paper is organised as follows. The remainder of this Section has a collection of preliminary definitions and notational convention, as well as the statement of the governing partial differential equations. The weak formulation and proofs of the uniform boundedness of the bilinear forms and suitable  inf-sup conditions are shown in Section <ref>. The fixed-point analysis of the coupled problem is carried out in Section <ref>. Section <ref> then focuses on the Galerkin discretisation, including its well-posedness analysis and definition of specific FE subspaces that provide momentum and mass conservativity. In Section <ref> we show a Céa estimate and using appropriate approximation properties we derive optimal a priori error bounds including also the higher order case. The definition of a residual a posteriori error estimator and the proofs of its reliability and efficiency are presented in Section <ref>. We conclude in Section <ref> with some numerical tests that both validate and underline the theoretical properties of the proposed discretisations. 

Notation and preliminaries.
Let ^2(Ω) be the set of all square-integrable functions in Ω⊂^d where d ∈{2,3} is the spatial dimension, and denote by ^2(Ω)=^2(Ω)^d its vector-valued counterpart and by ^2(Ω)=^2(Ω)^d× d its tensor-valued counterpart. We also write  

    :={∈^2(Ω):  = -^ t},

to represent the skew-symmetric tensors in Ω with each  component being square-integrable. Standard notation will be employed for Sobolev spaces ^m(Ω) with m≥ 0 (and we note that ^0(Ω)=^2(Ω)). Their norms and seminorms are denoted as ·_m,Ω and |·|_m,Ω, respectively (as well as for their vector and tensor-valued counterparts ^m(Ω), ^m(Ω)) see, e.g., <cit.>. 

As usual 𝕀 stands for the identity tensor in ^d× d, 
and |·| denotes the Euclidean norm  in ^d. Also, for any vector field =(v_i)_i=1,d  we set the gradient and divergence operators as
:= (∂ v_i/∂ x_j)_i,j=1,d and  ÷ := ∑_j=1^d ∂ v_j/∂ x_j. 
In addition, for any tensor fields =(τ_ij)_i,j=1,d
and = (ζ_ij)_i,j=1,d, we let  be the divergence operator ÷ acting along the rows of , and define the transpose, the trace, the tensor inner product, and the deviatoric tensor  as
^ := (τ_ji)_i,j=1,d, 
() := ∑_i=1^dτ_ii,  
: := ∑_i,j=1^nτ_ijζ_ij, and 
^ :=  - 1/n (), respectively. 
We also recall the Hilbert space

    (÷;Ω) := {∈^2(Ω): ÷ ∈L^2(Ω)},

with norm _÷;Ω^2:=_0,Ω^2+÷ _0,Ω^2, and introduce its tensor-valued version 

    (;Ω) := {∈^2(Ω):  ∈^2(Ω)}.








Governing equations. 
Let us consider a fully-saturated poroelastic medium (consisting of a mechanically isotropic and homogeneous fluid-solid mixture) occupying the open and bounded domain Ω in ^d, the Lipschitz boundary ∂Ω is partitioned into disjoint sub-boundaries ∂Ω:= Γ_∪Γ_, and it is assumed for the sake of simplicity that both sub-boundaries are non-empty |Γ_|·|Γ_|>0.
The symbol  will stand for the unit outward normal vector on the boundary. Let ∈^2(Ω) be a prescribed body force per unit of volume (acting on the fluid-structure mixture) and let g ∈ L^2(Ω) be a net volumetric fluid production rate.  

The balance of linear momentum for the solid-fluid mixture is written as 

    -  =     in Ω,

with  being the total Cauchy stress tensor of the mixture (sum of the effective solid and fluid stresses), whose dependence on strain and on fluid pressure  is given by the constitutive assumption (or effective stress principle) 

    =  () -α p 𝕀    in Ω.

Here the skeleton displacement vector  from the position ∈Ω is an unknown, the tensor () := 1/2 ( + []^ t ) is the infinitesimal strain, by  we denote the fourth-order elasticity tensor, also known as Hooke's tensor (symmetric and positive definite and  characterised by := λ()𝕀 + 2μ), 𝕀 is the identity second-order tensor, λ and μ are the Lamé parameters (assumed constant and positive), 0≤α≤ 1 is the Biot–Willis parameter, and p denotes the Darcy fluid pressure (positive in compression), which is an unknown in the system. 

We also consider the balance of angular momentum, which in this context states that the total poroelastic stress is a symmetric tensor 

= ^ t.

To weakly impose it, it is customary to use the rotation tensor

    = 1/2 ( - []^ t ) =  - ().


The fluid content (due to both fluid saturation and local volume dilation) is given by 

    ζ = c_0 p + α÷,

where c_0 ≥ 0 is the constrained specific storage 
coefficient. 

Using Darcy's law to describe the discharge velocity in terms of the fluid pressure gradient,  the balance of mass for the total amount of fluid is ∂_tζ - ÷ (κ∇ p)  = g in Ω× (0,t_end), where κ is the intrinsic permeability 
of 
the medium, a nonlinear function of the porosity.  In turn, in the small strains limit the porosity can be approximated by a linear function of the fluid content ζ (see for example <cit.>), and so, thanks to 
(<ref>), we can simply write κ( (),p).  Furthermore, after a backward Euler semi-discretisation in time with a constant time step and rescaling appropriately, we only consider the type of equations needed to solve at each time step and therefore we will concentrate on the form 

    c_0 p + α () - ÷ (κ( (),p) ∇ p)  = g     in Ω.

Typical constitutive relations for permeability are, e.g., exponential or Kozeny–Carman type (cf.  <cit.>) 

    κ( (),p) = k_0/μ_f𝕀 + k_1/μ_fexp(k_2 (c_0p + α ()))𝕀,   κ( (),p) =   k_0/μ_f𝕀 +  k_1(c_0p + α ())^3/μ_f(1-(c_0p + α ()))^2𝕀,

where μ_f denotes the viscosity of the interstitial fluid and k_0,k_1,k_2 are model constants. We note that in the case of incompressible constituents one has c_0 = 0 and α = 1, indicating that permeability depends only on the dilation () = ÷ (see, e.g., <cit.>). We also note that even in such a scenario (of incompressible phases) the overall mixture is not necessarily incompressible itself. 

More precise assumptions on the behaviour of the permeability are postponed to Section <ref>.

Next we note that from (<ref>) we can obtain


    = (dλ+2μ)÷ -dα p 
    
    ^-1 + αdλ+2μ p 𝕀  =  ()     in Ω.



Then, from the first equation in 

(<ref>) we get 

    () = 1dλ + 2μ + dαdλ+2μ p,

and therefore 
the dependence of κ on () and p (cf.  (<ref>)) can be written in terms of  and p as follows

    [                                   κ(,p) = k_0/μ_f𝕀 + k_1/μ_fexp(k_2dλ+2μ( (c_0(dλ+2μ)+dα^2)  p + α) )𝕀,; κ(,p) =   k_0/μ_f𝕀 +  k_1( (c_0(dλ+2μ)+dα^2)  p + α)^3/(dλ+2μ)μ_f(dλ+2μ-((c_0(dλ+2μ)+dα^2)  p + α))^2𝕀. ]

In addition, putting together the second equation in (<ref>) and (<ref>) we obtain:

    ^-1 + αdλ+2μ p 𝕀 =  -    in Ω.


Finally, we introduce the discharge flux   as an unknown defined by the constitutive relation 

    κ(,p)^-1 = ∇ p,

and combining (<ref>) and (<ref>), we are able to rewrite the mass balance equation as 

    c_0 p + αdλ+2μ + dα^2dλ+2μ p - ÷  = g     in Ω.


To close the system, we consider  mixed boundary conditions for a given _∈^1/2(Γ_) and p_∈^1/2(Γ_): 

    = _ p = p_Γ_,     · = 0  = 0Γ_.








§ WEAK FORMULATION AND PRELIMINARY PROPERTIES



 §.§ Derivation of weak forms

Let us define the following spaces

    _(;Ω) := {∈(;Ω): =0Γ_},  _(÷;Ω) := {∈(÷;Ω): ·=0 Γ_}.

We test equation (<ref>) against ∈^2(Ω),  equation (<ref>) against ∈_(;Ω), impose the symmetry of  weakly, test equation (<ref>) against q∈^2(Ω),   equation (<ref>) against ∈_(÷;Ω), integrate by parts and using the boundary conditions (<ref>) naturally, and then reorder the resulting equations. Then we arrive at 

    ∫_Ω^-1: + αdλ+2μ∫_Ω p   + ∫_Ω· + ∫_Ω:   = ⟨,_⟩_Γ_     ∀ ∈_(;Ω),
    ∫_Ω·    =  -∫_Ω·     ∀ ∈^2(Ω),
    ∫_Ω:   =0      ∀ ∈,
    ∫_Ωκ(,p)^-1· +  ∫_Ω  p ÷   = ⟨·,p_⟩_Γ_     ∀ ∈_(÷;Ω),
    (c_0 +dα^2dλ+2μ)∫_Ω p q +  αdλ+2μ∫_Ω q  - ∫_Ωq÷   = ∫_Ω  g  q      ∀  q∈^2(Ω),

where ⟨·,·⟩_Γ_ denotes the duality pairing between ^-1/2(Γ_) and its dual ^1/2(Γ_) with respect to the inner product in L^2(Γ_), and we use the same notation, ⟨·,·⟩_Γ_, in the vector-valued case. 

Introducing bilinear forms a:_(;Ω)×_(;Ω)→, b:_(;Ω)×^2(Ω) ×→, c:_(;Ω)×^2(Ω)→, the nonlinear form a_, p:_(÷;Ω)×_(÷;Ω)→, and the bilinear weak forms b:_(÷;Ω) ×^2(Ω) → and c:^2(Ω) ×^2(Ω) →, defined by

    [                         a(,) := ∫_Ω^-1:,   
    b(,(,)) := ∫_Ω·  + ∫_Ω:,   c(,q):=αdλ+2μ∫_Ω q ,; [3ex]
     a_, p(,) := ∫_Ωκ(, p)^-1· ,   b(,q) := ∫_Ω  q ÷,   c(p,q):=(c_0 +dα^2dλ+2μ)∫_Ω p q , ]

respectively, and  linear functionals H∈_(;Ω)', F∈ (^2(Ω)×)', H∈_(÷;Ω)',  G∈^2(Ω)'  

    H(): = ⟨,_⟩_Γ_,    
    F(,):= -∫_Ω·,    H(): = ⟨·,p_⟩_Γ_,    
    G(q):= -∫_Ω  g  q ,

we arrive at: 

find (,,,,p)∈_(;Ω)×^2(Ω)××_(÷;Ω)×^2(Ω), such that:


    a(,)      +   b(,(,))     +   c(,p)   =  H()   ∀ ∈_(;Ω),
     b(,(,))            = F(,)   ∀ ∈^2(Ω), ∀ ∈,
     a_,p(,)      +   b(,p)        =  H()   ∀ ∈_(÷;Ω),
     b(,q)     -   c(p,q)     -   c(,q)     = G(q)   ∀  q∈^2(Ω).





 §.§ Stability properties and suitable inf-sup conditions

For the sake of the analysis, we allow the permeability κ(,p) to be anisotropic but still require κ(,p)^-1 to be uniformly positive definite  in ^∞(Ω) and  Lipschitz continuous with respect to p∈^2(Ω). That is, there exist positive constants κ_1,κ_2 such that 

    κ_1||^2 ≤^ tκ(·,·)^-1, 
    
        κ(·,p_1)^-1 - κ(·,p_2)^-1_^∞(Ω)≤κ_2 p_1 -p_2_0,Ω,

for all ∈ℝ^d∖{}, and for all p_1,p_2∈^2(Ω).

We start by establishing the boundedness of  the bilinear forms a, b, c, b, c:


    |a(,)|≤1μ_;Ω_;Ω,    |b(,(,))|≤_;Ω(_0,Ω + _0,Ω),  
    |c(,q)|≤γ_;Ωq_0,Ω,
    
    | b(,q)|≤_÷;Ωq_0,Ω ,    | c(p,q)|≤γp_0,Ωq_0,Ω,


where

    γ:=α√(d)dλ+2μγ:=c_0 +dα^2dλ+2μ.


On the other hand, using Hölder's and trace inequalities we can readily observe that the right-hand side functionals are all bounded

    |H()|  ≤__1/2,Γ__;Ω, 
         |F(,)|  ≤_0,Ω_0,Ω≤_0,Ω(_0,Ω + _0,Ω), 
     
    | H()|  ≤p__1/2,Γ__÷;Ω,     |G(q)|  ≤g_0,Ωq_0,Ω.

Let us now denote by  and  the kernels of b and b, respectively. They are characterised, respectively, as


    = {∈_(;Ω) :    =  = ^ tΩ},
    
    
        = {∈_(÷;Ω) :  ÷  = 0   in  Ω}.

From <cit.> we easily deduce that there exists c_a>0 such that

    a(,) ≥ c_a _;Ω^2   ∀ ∈.

The following inf-sup conditions are well-known to hold (see, e.g., <cit.>):

 
    sup_≠∈_(;Ω) b(,(,))/_;Ω   ≥β(_0,Ω +_0,Ω)
      ∀ (,)∈^2(Ω)×,
    sup_≠∈_(÷;Ω) b(,q)/ψ_÷;Ω   ≥βq_0,Ω  ∀ q∈^2(Ω).
   

Finally, we observe that c is elliptic over ^2(Ω)

    c(q,q) ≥γ q_0,Ω^2.





§ ANALYSIS OF THE COUPLED PROBLEM

We now use a combination of the classical Babuška–Brezzi and Banach fixed-point theorems to establish the well-posedness of (<ref>) under appropriate assumptions on the data.



 §.§ A fixed‑point operator

We adopt a similar approach to, e.g., <cit.>. 


Firs, we define a closed ball of ^2(Ω) centred at the origin and of given radius r>0

    := { p ∈^2(Ω)  :   p_0,Ω≤ r }.

Then, for a given (, p)∈_(;Ω)×, thanks to the assumptions on the nonlinear permeability, we can infer that the form a_, p (cf. (<ref>)) is continuous, as well as coercive over 


    | a_, p(,)|    ≤ C_ a _÷;Ω_÷;Ω,
     a_, p(,)    ≥κ_1 _÷;Ω^2   ∀ , ∈.




Then, we define the auxiliary operators 
:⊆^2(Ω)→_(;Ω)×(^2(Ω)×) and : _(;Ω)×→_(÷;Ω)×^2(Ω), given by

    ( p):=(R_1( p),(R_2( p),R_3( p)))=(,(,))   ∀  p∈,

with (,(,))∈_(;Ω)×(^2(Ω)×) satisfying

    [            a(,)     +   b(,(,)) =  H() - c(, p)       ∀ ∈_(;Ω),;  [1ex]
     b(,(,))                          = F(,)   ∀ (,)∈^2(Ω)×, ]

and

    (, p):=(S_1(, p),S_2(, p) )=(,p)   ∀  (, p)∈_(;Ω)×,

where (,p) is such that

    [       a_, p(,)      +   b(,p)          = H()     ∀ ∈_(÷;Ω),;   [1ex]
     b(,q)     -   c(p,q) = G(q) + c(,q)   ∀  q∈^2(Ω) . ]


By virtue of the above, by defining the operator :⊆^2(Ω)→^2(Ω) as

    ( p):=S_2(R_1( p), p),

it is clear that (,,,,p) is a solution to (<ref>) if and only if p∈ solves the fixed-point problem 

    (p)=p.

Thus, in what follows, we focus on proving the unique solvability of (<ref>). 



 §.§ Well‑definedness of Lg

From the definition of  in  (<ref>) it is evident that its well-definedness requires the well-posedness of problems (<ref>) and (<ref>). We begin by analysing that of (<ref>).

Let	p ∈ (cf. (<ref>)). Then, there exists a unique (,(,)) ∈_(;Ω)×^2(Ω)× solution to (<ref>).
In addition, there exist C_1, C_2>0, such that

    [                   _;Ω≤ C_1  (__1/2,Γ_ + _0,Ω) + 1c_a γ  p_0,Ω,; _0,Ω + _0,Ω≤ C_2  (__1/2,Γ_ + _0,Ω) + 1β(1 + 1μ c_a) γ  p_0,Ω. ]



It is a direct consequence of the Babuška–Brezzi theory <cit.>, using  (<ref>) and (<ref>)  with

    C_1:=(1c_a + 1β)(1 + 1μ c_a)  C_2:= 1β(1 + 1μ c_a)(1 + 1μβ);

we omit further details.


Next, we provide the well-definedness of , or equivalently, the well-posedness of (<ref>).

Let	(, p)∈_(;Ω)×. Then, there exists a unique (,p) ∈_(÷;Ω)×^2(Ω) solution to (<ref>).
In addition, there exist C>0 such that

    _÷;Ω + p_0,Ω≤ C  (g_0,Ω + p__1/2,Γ_ + γ_;Ω ).



The existence of a unique solution (,p) to (<ref>) is straightforward given the properties of the forms a, b, and c. By examining (<ref>), (<ref>), and (<ref>), we can confirm that the assumptions of <cit.> are satisfied. In addition,  and p satisfy the following bounds 

    _÷;Ω   ≤( 1κ_1 +  C_1 +  C_1√(γ))(2max{ C_2,  C_3})^1/2 (g_0,Ω + p__1/2,Γ_ + γ _;Ω ),
    
    p_0,Ω   ≤1β(1 + κ_2 r (1κ_1 +  C_1 +  C_1√(γ)) (2max{ C_2,  C_3})^1/2) (g_0,Ω + p__1/2,Γ_ + γ _;Ω ),


where

    C_1:=1β(1 + C_ aκ_1),   C_2:= 1κ_1 +  C_1 + γ  C_1^2   C_3:=  C_1 ( 1 + C_ aβ +  C_1 C_ a^2 γβ^2);

and the above implies (<ref>). We leave out additional minor details.


Given r>0, let us assume that 

    C (1 + γ C_1 )( 
     g_0,Ω+ p__1/2,Γ_ + __1/2,Γ_ + _0,Ω) +  C c_a γ^2 r≤ r,

where C_1, C_1 and γ are defined in (<ref>), (<ref>) and (<ref>), respectively. Then, for a given p∈ (cf. (<ref>)), there exists a unique p ∈ such that (p) = p.


From Lemmas <ref> and <ref>, we ascertain that the operators  and , respectively, are well-defined, thereby ensuring the well-definition of . Furthermore, from  (<ref>) and (<ref>), for each p ∈, we deduce that

    ( p)_0,Ω    = S_2(R_1( p), p)_0,Ω
       ≤ C  (g_0,Ω + p__1/2,Γ_ ) +  C γ (R_1( p)_;Ω
       ≤ C  (g_0,Ω + p__1/2,Γ_ ) +  C γ C_1  (__1/2,Γ_ + _0,Ω) +  Cc_a γ^2 p_0,Ω,

this, combined with assumption (<ref>), implies ()⊆, which concludes the proof.



Another option for defining the operator  (see (<ref>)) is to introduce the perturbation c on the right-hand side of the system, given by 

    [                  a_, p(,)                 +   b(,p)                     = H()                ∀ ∈_(÷;Ω),;              [1ex]
     b(,q)                           = G(q) + c(,q) +  c( p,q)               ∀  q∈^2(Ω). ]

But in this case, the assumption of small data in (<ref>) (as well as in other instances, later on) would also involve the storativity parameter c_0, making the analysis slightly more restrictive.




 §.§ Existence and uniqueness of weak solution
 

We begin by establishing two lemmas deriving conditions under which the operator  is a contraction.

Given p_1,  p_2, ∈, the following estimate holds

    R_1( p_1) - R_1( p_2)_;Ω≤1c_a γ  p_1-  p_2_0,Ω.



Let (_1,(_1,_1)), (_2,(_2,_2)) ∈_(;Ω)×(^2(Ω)×), such that ( p_1)=(_1,(_1,_1)) ad ( p_2)=(_2,(_2,_2)). Then, from the definition of  (cf. (<ref>)), we have

    [             a(_1-_2,) +   b(,(_1-_2,_1-_2))   =  - c(, p_1-  p_2)             ∀ ∈_(;Ω),;   [1ex]
     b(_1-_2,(,))                                         = 0         ∀ (,)∈^2(Ω)×. ]

Since _1-_2 ∈ (cf. (<ref>)), taking = _1-_2 in (<ref>), and utilising the ellipticity of a on  (cf. (<ref>)) along with the bound of c (cf. (<ref>)), we obtain:

    c_a _1-_2_;Ω^2 ≤ a(_1-_2,_1-_2)  =  - c(_1-_2, p_1-  p_2) ≤ γ _1-_2_;Ω p_1-  p_2_0,Ω,

which concludes the proof.


Given (_1,  p_1), (_2,  p_2), ∈_(;Ω)×, the following estimate holds

    [                                                           S_2(_1,  p_1) - S_2(_2,  p_2)_0,Ω; ≤2κ_2  Cmin{γ, κ_1} (g_0,Ω + p__1/2,Γ_ + γ_2_;Ω)  p_1 -  p_2_0,Ω + 2min{γ, κ_1} γ _1-_2_;Ω. ]



Let (_1,p_1), (_2,p_2)∈_(÷;Ω)×^2(Ω), such that (_1, p_1)=(_1,p_1) and (_2, p_2)=(_2, p_2). Then, from the definition of  (cf. (<ref>)), and employing similar arguments to those in Lemma <ref>, we have

    a__1, p_1(_1,_1-_2) -  a__2, p_2(_2,_1-_2) +  c(p_1-p_2,p_1-p_2) = - c(_1-_2,p_1-p_2),

by adding ± a__1, p_1(_2,_1-_2) in the last equation, we obtain

    a__1, p_1(_1-_2,_1-_2) +  c(p_1-p_2,p_1-p_2)
    =   a__2, p_2(_2,_1-_2) -  a__1, p_1(_2,_1-_2) - c(_1-_2,p_1-p_2).

Then, using the the first assumption for κ (cf. (<ref>)), the ellipticity of c (see (<ref>)), the definition of a_, p (cf. (<ref>)) and the continuity of the form c (see (<ref>)), we deduce

    [ κ_1 _1-_2_0,Ω^2 + γp_1-p_2_0,Ω^2 ≤ a__1, p_1(_1-_2,_1-_2) +  c(p_1-p_2,p_1-p_2);             = ∫_Ω (κ(_2, p_2)^-1-κ(_1, p_1)^-1) _2 · (_1-_2) - c(_1-_2,p_1-p_2);   ≤κ(_2, p_2)^-1- κ(_1, p_1)^-1_^∞(Ω) _2_0,Ω_1-_2_0,Ω +  γ _1-_2_;Ωp_1-p_2_0,Ω. ]

From the last equation, by utilising the second assumption regarding κ (see (<ref>)), we obtain

    [ 1/2min{γ,κ_1}(_1-_2_0,Ω + p_1-p_2_0,Ω)^2 ≤κ_1 _1-_2_0,Ω^2 + γp_1-p_2_0,Ω^2;                  ≤κ_2 p_2- p_1_0,Ω_2_0,Ω_1-_2_0,Ω +  γ _1-_2_;Ωp_1-p_2_0,Ω;        ≤(κ_2 p_2- p_1_0,Ω_2_0,Ω  +  γ _1-_2_;Ω) (_1-_2_0,Ω + p_1-p_2_0,Ω), ]

the last, together with the fact that _2 satisfies (<ref>), leads to the following bound 

    [ 1/2min{γ,κ_1}(_1-_2_0,Ω + p_1-p_2_0,Ω) ≤κ_2 p_2- p_1_0,Ω_2_0,Ω  +  γ _1-_2_;Ω;            ≤κ_2 p_2- p_1_0,Ω C  (g_0,Ω + p__1/2,Γ_ + γ_2_;Ω )  +  γ _1-_2_;Ω, ]

and this yields (<ref>), concluding the proof.


The following theorem presents the main result of this section, establishing the existence and uniqueness of the solution to the fixed-point problem (<ref>), or equivalently, the well-posedness of problem (<ref>).

Given r>0, assume that ∈^2(Ω), g ∈ L^2(Ω), _∈^1/2(Γ_), p_∈^1/2(Γ_) and γ satisfies


    [ 2 max{1, κ_2}min{γ, κ_1 , r} { C (1 + C_1γ)(g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω)+γ^2c_a(1κ_2 + C   r)}   < 1. ]

Then,  (cf. (<ref>)) has a unique fixed point p∈. Equivalently,  (<ref>) has a unique solution (,,,,p)∈_(;Ω)×^2(Ω)××_(÷;Ω)×. In addition, 
there exists C>0, such that

    _;Ω + _0,Ω + _0,Ω + _÷;Ω + p_0,Ω≤ C (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω + γ r ).




Recall that (<ref>) ensures the well-definedness of . Let p_1,  p_2,  p_1,  p_2 ∈, such that ( p_1)=p_1 and ( p_2)=p_2. From the definition of  (see (<ref>)), and the estimates (<ref>) and (<ref>), we deduce

    p_1-p_2_0,Ω = ( p_1)-( p_2)_0,Ω = S_2(R_1( p_1), p_1)-S_2(R_1( p_2), p_2)_0,Ω
       ≤2κ_2  Cmin{γ, κ_1} (g_0,Ω + p__1/2,Γ_ + γR_1( p_2)_;Ω)  p_1 -  p_2_0,Ω+ 2min{γ, κ_1} γ R_1( p_1)-R_1( p_2)_;Ω
       ≤2κ_2  Cmin{γ, κ_1} (g_0,Ω + p__1/2,Γ_) p_1 -  p_2_0,Ω + 2κ_2  Cmin{γ, κ_1}γR_1( p_2)_;Ω p_1 -  p_2_0,Ω
           + 2c_a min{γ, κ_1} γ^2  p_1-  p_2_0,Ω,

the above, along with the fact that R_1( p_2) satisfies (<ref>) and p_2∈, implies

    p_1-p_2_0,Ω≤2κ_2  Cmin{γ, κ_1} (g_0,Ω + p__1/2,Γ_) p_1 -  p_2_0,Ω+ 2c_a min{γ, κ_1} γ^2  p_1-  p_2_0,Ω
                    +2κ_2  Cmin{γ, κ_1}γ( C_1  (__1/2,Γ_ + _0,Ω) + 1c_a γ r  )  p_1 -  p_2_0,Ω
         ≤2min{γ, κ_1} {κ_2  C (1 + C_1γ)(g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω)+γ^2c_a(1+κ_2  C   r)} p_1-  p_2_0,Ω,

which together with (<ref>) and the Banach fixed-point theorem yields that  has a unique fixed point in .
Finally,  (<ref>)  is derived analogously to the estimates in (<ref>) and (<ref>), which completes the proof.


The operator  (see (<ref>)) could be also defined, for example :→, with := { (, p) ∈_(;Ω)×^2(Ω)  :  _;Ω +  p_0,Ω≤ r } and (, p):=(R_1( p),S_2(, p))=(,p), with R_1 and S_2 defined as in (<ref>) and (<ref>), respectively.




§ FINITE ELEMENT DISCRETISATION

In this section, we present and analyse the Galerkin scheme for problem (<ref>). It is worth mentioning upfront that the well-posedness analysis can be straightforwardly extended from the continuous problem to the discrete case. Therefore, we omit many of the details.




 §.§ Finite element spaces and Galerkin scheme

Let us consider a  regular partition 𝒯_h of Ω̅ made up of triangles K (in ℝ^2) or tetrahedra K (in ℝ^3) of diameter h_K, and denote the mesh size by h := max{ h_K:  K ∈𝒯_h}. 


Given an integer ℓ≥ 0 and K ∈𝒯_h, we first let P_ℓ(K) be the space of polynomials of degree ≤ℓ defined on K, whose vector and tensor versions are denoted _ℓ(K)  :=  [P_ℓ(K)]^d and ℙ_ℓ(K)
 = [P_ℓ(K)]^d × d, respectively. Also, we let 𝐑𝐓_ℓ(K)  := _ℓ(K) ⊕P_ℓ(K) be the local Raviart–Thomas space of order ℓ defined on K, where  stands for a generic vector in ^d, and denote by ℝ𝕋_k(K) the tensor-valued counterpart of this space. 

For each K∈𝒯_h we consider the bubble space of order k, defined as 

    𝐁_k(K):=
     𝐜𝐮𝐫𝐥^t(b_KP_k(K))   in  ℝ^2,
    ∇× (b_K𝐏_k(K))   in  ℝ^3,

where b_K is a suitably normalised cubic polynomial on K, which vanishes on the boundary of K (see <cit.>).

We recall the classical PEERS_k elements (cf. <cit.>) to define the discrete subspaces for the stress tensor , the displacement , and the rotation tensor 

    ^_h   :={_h ∈_(;Ω):  _h|_K∈ℝ𝕋_k(K)⊕[𝐁_k(K)]^d  ∀  K∈𝒯_h},
    ^_h   := {_h ∈^2(Ω):   _h|_K∈P_k(K)  ∀  K∈𝒯_h},
    ^_h   :={_h∈∩ℂ(Ω)  and  _h|_K∈ℙ_k+1(K)  ∀  K∈𝒯_h},

and the following estimates are proven for the  PEERS_k elements (cf. <cit.>)


    sup_≠_h∈^_hb(_h,(_h,_h))/_h_;Ω   ≥β^*(_h_0,Ω +_h_0,Ω)
      ∀ (_h,_h)∈^_h×^_h,
     a(_h,_h)    ≥ c_a _h_;Ω^2   ∀ _h∈_h,

where _h denotes the discrete kernel of b, that is

    _h   :=  {_h∈^_h :   b(_h,(_h,_h)) = 0   ∀ (_h,_h)∈^_h×^_h}.


 
Additionally, for  and the pressure p, we define the FE subspaces

    ^_h   :={_h∈_(÷;Ω):  _h|_K∈𝐑𝐓_k(K)  ∀  K∈𝒯_h},
    ^p_h   := {q_h ∈^2(Ω):   q_h|_K∈P_k(K)  ∀  K∈𝒯_h},

and it is well known that b satisfies the inf-sup condition (see, e.g.,  <cit.>)

    sup_≠_h∈^_h b(_h,q_h)/_h_÷;Ω≥β^* q_h_0,Ω  ∀ q_h∈^p_h.


Note that it is of course possible to consider other conforming and inf-sup stable spaces such as Arnold–Falk–Winther and Brezzi–Douglas–Marini instead of (<ref>) and (<ref>), respectively.   
The Galerkin scheme for   (<ref>) reads: find (_h,_h,_h,_h,p_h)∈^_h×^_h×^_h×^_h×^p_h, such that:

    [               a(_h,_h)      +   b(_h,(_h,_h))          +   c(_h,p_h)               =  H(_h)              ∀ _h∈^_h,;   [1ex]
     b(_h,(_h,_h))                                                           = F(_h,_h)     ∀(_h,_h)∈^_h×^_h ,; [1ex]
     a__h,p_h(_h,_h)          +   b(_h,p_h)                                      =  H(_h)              ∀ _h∈^_h,;       [1ex]
     b(_h,q_h)         -   c(p_h,q_h)          -   c(_h,q_h)               = G(q_h)           ∀  q_h∈^p_h. ]






 §.§ Analysis of the discrete problem

In this section, we analyse the Galerkin scheme (<ref>). It's worth noting that establishing well-posedness can be readily achieved by extending the results derived for the continuous problem to the discrete setting.

Firstly, and similarly to the continuous case, we define the following set

    _h := { p_h ∈^p_h :   p_h_0,Ω≤ r }.

Next, for a fixed p_h in _h, we have that the bilinear form a__h,p_h satisfies

    a__h, p_h(_h,_h) ≥κ_1 _h_÷;Ω^2   ∀ _h∈_h,

where _h is the discrete kernel of b

    [ _h   :=  {_h∈^_h :   b(_h,q_h) = 0   ∀ q_h∈^p_h}. ]


Additionally, we define the discrete operators _h:_h⊆^p_h→^_h×(^_h×^_h) and _h: ^_h×_h →^_h×^p_h, respectively, by

    _h( p_h):=(R_1,h( p_h),(R_2,h( p_h),R_3,h( p_h)))=(_h,(_h,_h))   ∀  p_h∈_h,

where (_h,(_h,_h))∈^_h×(^_h×^_h) is the unique solution of 

    [              a(_h,_h)     +   b(_h,(_h,_h)) =  H(_h) - c(_h, p_h)             ∀ _h∈^_h,;  [1ex]
     b(_h,(_h,_h))                                  = F(_h,_h)    ∀ (_h,_h)∈^_h×^_h, ]

and

    _h(_h, p_h):=(S_1,h(_h, p_h),S_2,h(_h, p_h) )=(_h,p_h)   ∀  (_h, p_h)∈^_h×_h,

where (_h,p_h) is the unique tuple in ^_h××^p_h such that 

    [     a__h, p_h(_h,_h)        +   b(_h,p_h)             =  H(_h)            ∀ _h∈^_h,;     [1ex]
     b(_h,q_h)       -   c(p_h,q_h) = G(q_h) + c(_h,q_h)         ∀  q_h∈^p_h. ]


Employing properties (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>) and proceeding exactly as for the continuous case (Lemmas <ref> and <ref>), it can be
easily deduced that both operators are well-defined. Then, analogously to the continuous case, we define the following fixed-point operator 

    _h:_h⊆^p_h→^p_h,   p_h↦_h( p_h):=S_2,h(R_1,h( p_h), p_h),

which is clearly well-defined (since R_h and S_h are). Further, it can be easily deduced that _h(_h)⊆_h if 

    C^* (1 + γ C_1^* )( 
     g_0,Ω+ p__1/2,Γ_ + __1/2,Γ_ + _0,Ω) +  C^* c_a γ^2 r≤ r,

where C^* and C_1^* (depending on c_a, μ, κ_1, κ_2, C_ a, β^*, β^*) are the discrete versions of the constants C and C_1 (cf. (<ref>) and (<ref>)). Finally, it is clear that (_h,_h,_h,_h,p_h) is a solution to (<ref>) if and only if p_h satisfies 


    _h(p_h)=p_h.

The main outcome of this section is presented in the following theorem, establishing the existence and uniqueness of a solution to the fixed-point problem (<ref>),  equivalently proving the well-posedness of problem (<ref>).

Given r>0, assume that 

the data and γ satisfy

    2 max{1, κ_2}/min{γ, κ_1 , r} { C^* (1 + C_1^*γ)(g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω)+γ^2c_a(1/κ_2 + C^*   r)}   < 1.

Then,  _h (cf. (<ref>)) has a unique fixed point p_h∈_h. Equivalently, problem (<ref>) has a unique solution (_h,_h,_h,_h,p_h)∈^_h×^_h×^_h×^_h×_h.
In addition, 

there exists C^*>0, such that

    _h_;Ω + _h_0,Ω + _h_0,Ω + _h_÷;Ω + p_h_0,Ω
    
    ≤ C^* (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω + γ r ).



First, we observe that, similar to the continuous case (as seen in the proof of Theorem <ref>), assumption (<ref>) ensures the well-definedness of _h and that _h(_h)⊆_h. Now, by adapting the arguments used in Section <ref> (cf.  Lemmas <ref> and <ref>), one can derive the following estimates


    R_1,h( p_1) - R_1,h( p_2)_;Ω   ≤1c_a γ  p_1-  p_2_0,Ω,
    S_2,h(_1,  p_1) - S_2,h(_2,  p_2)_0,Ω   ≤2κ_2  C^*/min{γ, κ_1} (g_0,Ω + p__1/2,Γ_ + γ_2_;Ω)  p_1 -  p_2_0,Ω
          + 2/min{γ, κ_1} γ _1-_2_;Ω,

for all p_1,  p_2∈_h and _1, _2 ∈^_h, which together with the definition of _h (see (<ref>)), yield

    _h( p_1)-_h( p_2)_0,Ω   ≤2min{γ, κ_1} {κ_2  C^* (1 + C_1^*γ)(g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω)
                      + γ^2c_a(1+κ_2  C^*   r)} p_1-  p_2_0,Ω,

for all p_1,  p_2∈_h. In this way, using estimate (<ref>) we obtain that _h is a contraction
mapping on _h, thus problem (<ref>), or equivalently (<ref>) is well-posed.
Finally, analogously to the proof of Theorem  <ref> (see also Lemmas <ref> and <ref>) we can obtain (<ref>), which concludes the proof.






§ A PRIORI ERROR ESTIMATES


In this section, we aim to provide the convergence of the Galerkin scheme (<ref>) and derive the corresponding rate of convergence. From now on we assume that the hypotheses of Theorem <ref> and Theorem <ref> hold.



 §.§ Preliminaries


Let the tuples (,,,,p)∈_(;Ω)×^2(Ω)××_(÷;Ω)×^2(Ω) and (_h,_h,_h,_h,p_h)∈^_h×^_h×^_h×^_h×^p_h be the unique solutions of (<ref>) and (<ref>), respectively.


Let us write _ =  - _h, _ =  - _h, _ =  - _h, _ =  - _h and _p = p - p_h.
As usual, for a given (_h,(_h,_h) )∈^_h×(^_h×^_h) and (_h,q_h)∈^_h×^p_h, we shall then decompose these errors into

    _ = _ + _,    _ = _ + _,    _ = _ + _,    _ = _ + _,    _p = ξ_p + χ_p,

with
_ :=  - _h,
_ := _h - _h,  
_ :=  - _h, 
_ := _h - _h, 
_ :=  - _h,  
_ := _h - _h, 
_ :=  - _h,  
_ := _h - _h, 
ξ_p := p - q_h, and  
χ_p := q_h - p_h.


Considering the first two equations of  problems (<ref>) and (<ref>), the following identities hold

    [           a(,)    +   b(,(,)) =  H() - c(,p)      ∀ ∈_(;Ω),; [1ex]
     b(,(,))                        = F(,) ∀ ∈^2(Ω), ∀ ∈, ]

and

    [             a(_h,_h)    +   b(_h,(_h,_h)) =  H(_h) - c(_h,p_h)            ∀ _h∈^_h,; [1ex]
     b(_h,(_h,_h))                                = F(_h,_h)    ∀(_h,_h)∈^_h×^_h. ]

From these relations we can obtain that for all (_h,(_h,_h))∈^_h×(^_h×^_h), there holds

    [             a(_,_h)     +   b(_h,(_,_))      =  - c(_h,_p),; [1ex]
     b(_,(_h,_h))                                     = 0 ]

which together with the definition of the errors in (<ref>), implies that 

    [                          a(_,_h) + b(_h,(_,_)) + b(_,(_h,_h)); = -a(_,_h) - b(_h,(_,_)) - b(_,(_h,_h)) - c(_h,_p) - c(_h,_p) ]

for all (_h,(_h,_h))∈^_h×(^_h×^_h).

Next, considering the last two equations of both problems (<ref>) and (<ref>), we obtain

    [        a_,p(,)      +   b(,p)         =  H()     ∀ ∈_(÷;Ω),;   [1ex]
     b(,q)     -   c(p,q) = G(q) + c(,q)    ∀  q∈^2(Ω). ]

and

    [      a__h,p_h(_h,_h)        +   b(_h,p_h)             =  H(_h)            ∀ _h∈^_h,;     [1ex]
     b(_h,q_h)       -   c(p_h,q_h) = G(q_h) + c(_h,q_h)         ∀  q_h∈^p_h. ]

Then, using arguments similar to those in Lemma <ref>, by adding  ± a__h, p_h(,_h), we have

    [ a__h,p_h(__h,_h) +  b(_h,_p_h) +  b(__h,q_h) - c(_p_h,q_h)= - ∫_Ω(κ(,p)^-1-κ(_h,p_h)^-1) ·_h + c(_,q_h), ]

which together with (<ref>), implies that

    [             a__h,p_h(_,_h) +  b(_h,χ_p_h) +  b(_,q_h) - c(χ_p,q_h)+  a__h,p_h(_,_h); =-  b(_h,ξ_p) -  b(_,q_h) + c(ξ_p,q_h) - ∫_Ω(κ(,p)^-1-κ(_h,p_h)^-1) ·_h + c(_,q_h). ]




 §.§ Derivation of Céa estimates


There exist C_3^*,  C_4^*>0, independent of h, such that

    __;Ω + __0,Ω + __0,Ω≤ C_3^* (__;Ω + __0,Ω + __0,Ω  + _p_0,Ω  ) + C_4^* γχ_p_0,Ω.
    


From the properties of a and b (refer to (<ref>) and (<ref>)), and <cit.>, we derive the following discrete global inf-sup condition

    __;Ω + __0,Ω + __0,Ω
    
    ≤ (C_1^* +C_2^*) sup_≠(_h,_h,_h)∈^_h×^_h×^_ha(_,_h) + b(_h,(_,_)) + b(_,(_h,_h))/_h_;Ω + _h_0,Ω + _h_0,Ω,

with C_1^*, C_2^*>0 independent of h, are the discrete version of the constants C_1, C_2 defined in (<ref>). Then, combining the last inequality with (<ref>), and the continuity properties of a and b (see (<ref>)), 
we obtain

    [                                                       __;Ω + __0,Ω + __0,Ω; ≤ (C_1^* + C_2^* )(1μ__;Ω + __0,Ω + __0,Ω + __;Ω + γχ_p_0,Ω + γξ_p_0,Ω  ), ]

which implies (<ref>) with C_3^*:=(C_1^* + C_2^* )(1/μ+1+γ) and C_4^*:=C_1^* + C_2^* , and
concludes the proof.




There exist C_5^*,  C_6^*>0, independent of h, such that

    [                            __÷;Ω + _p_0,Ω≤ C^*_5 (__÷;Ω + ξ_p_0,Ω + __;Ω); +  C^*_6 ( (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω + γ r ) χ_p_0,Ω + γ__;Ω). ]
    


Similarly to Lemma <ref>, using the properties of a, b and c (refer to (<ref>), (<ref>) and (<ref>)), and <cit.>, we derive the following discrete global inf-sup condition

    __÷;Ω + χ_p_0,Ω≤ 2 C^* sup_≠(_h,q_h)∈^_h×^p_h a__h,p_h(_,_h) + b(_h,χ_p) + b(_,q_h) -  c(χ_p,q_h)/_h_÷;Ω + q_h_0,Ω,

with C^* defined as in (<ref>). Then, using the equation (<ref>), the bound properties of a, b and c (see (<ref>) and (<ref>)), and the second assumption for κ (cf. (<ref>)), we obtain

    __÷;Ω + χ_p_0,Ω
       ≤ 2 C^* (C_ a__÷;Ω + ξ_p_0,Ω + __÷;Ω + γξ_p_0,Ω + κ(,p)^-1 - κ(_h,p_h)^-1_^∞(Ω)_÷;Ω + γ__;Ω)
       ≤ 2 C^* (C_ a__÷;Ω + ξ_p_0,Ω + __÷;Ω + γξ_p_0,Ω + κ_2_p_0,Ω_÷;Ω + γ__;Ω),

hence, using the fact that  satisfies (<ref>) and the error decomposition (<ref>), we have 

    __÷;Ω + χ_p_0,Ω   ≤ 2 C^* (C_ a__÷;Ω + ξ_p_0,Ω + __÷;Ω + γξ_p_0,Ω + κ_2ξ_p_0,Ω_÷;Ω  + γ__;Ω) 
            + 2 C^* ( κ_2χ_p_0,Ω_÷;Ω + γ__;Ω)
       ≤ C^*_5 (__÷;Ω + ξ_p_0,Ω + __;Ω)  
            + 2 C^* ( κ_2 C (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω + γ r ) χ_p_0,Ω + γ__;Ω),

the last equation implies (<ref>), with C^*_5 := 2 C^* ( C_ a+1 +γ +γ + κ_2 C (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω + γ r )     ) and C^*_6 := 2 C^*(κ_2 C + 1), and concludes the proof.



Assume that

    (C_4^* +  C^*_6 +  C^*_6 r) γ +  C^*_6 (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω )  ≤12,

with C_4^* and C^*_6 being the constants in Lemmas <ref> and <ref>. Furthermore, assume that the hypotheses of Theorem <ref> and Theorem <ref> hold. Let (,,,,p)∈_(;Ω)×^2(Ω)××_(÷;Ω)×^2(Ω) and (_h,_h,_h,_h,p_h)∈^_h×^_h×^_h×^_h×^p_h be the unique solutions of (<ref>) and (<ref>), respectively.
Then, there exists C_Céa>0, independent of h, such that

    [     __;Ω + __0,Ω + __0,Ω + __÷;Ω + _p_0,Ω; ≤ C_Céa ( (,,,,p), ^_h×^_h×^_h×^_h×^p_h). ]



Combining (<ref>) and (<ref>), and using the assumption (<ref>), we deduce 

    [ __;Ω + __0,Ω + __0,Ω + __÷;Ω + χ_p_0,Ω≤ C_3^* (__;Ω + __0,Ω + __0,Ω +  ξ_p_0,Ω  );                           +  C^*_5 (__÷;Ω + ξ_p_0,Ω + __;Ω) + 12χ_p_0,Ω + 12__;Ω. ]
 
And from the latter inequality we obtain

    [                           __;Ω + __0,Ω + __0,Ω + __÷;Ω + χ_p_0,Ω; ≤ 2(C_3^* +  C^*_5) (__;Ω + __0,Ω + __0,Ω + __÷;Ω +  ξ_p_0,Ω  ). ]
 
In this way, from (<ref>), (<ref>) and the triangle inequality we obtain

    [     __;Ω + __0,Ω + __0,Ω + __÷;Ω + _p_0,Ω≤__;Ω + __;Ω + __0,Ω + __0,Ω;                   + __0,Ω + __0,Ω + __÷;Ω + __÷;Ω + χ_p_0,Ω+  ξ_p_0,Ω; ≤ (2C_3^* + 2 C^*_5 + 1) (__;Ω + __0,Ω + __0,Ω + __÷;Ω +  ξ_p_0,Ω  ), ]

which combined with the fact that (_h,(_h,_h) )∈^_h×(^_h×^_h) and (_h,q_h)∈^_h×^p_h are arbitrary (see (<ref>)), concludes the proof.




 §.§ Rates of convergence

In order to establish the rate of convergence of the Galerkin scheme (<ref>), we first recall the following approximation properties 
associated with the FE spaces  specified in Section <ref>.


For each 0 < m ≤ k+1 and for each ∈^m(Ω)∩_(;Ω) with ∈^m(Ω), 
∈^m(Ω),  ∈^m(Ω)∩^2_(Ω),  ∈^m(Ω)∩_(;Ω) with ÷ ∈^m(Ω), and q∈^m(Ω), 
there holds 



    (,^_h)    := inf__h ∈^_h - _h_;Ω≲ h^m {_m,Ω +  _m,Ω},
    
    
    
    (,^_h)    := inf__h ∈^_h - _h_0,Ω≲  h^m _m,Ω,
    
    
    
    (,^_h)    := inf__h ∈^_h - _h_0,Ω≲  h^m _m,Ω,
    
    
    
    
    (,^_h)    := inf__h ∈^_h - _h_÷;Ω≲ h^m {_m,Ω + ÷ _m,Ω},
    
    
    
    (q,^p_h)    := inf_q_h ∈^p_hq - q_h_0,Ω≲  h^m q_m,Ω.


For (<ref>), (<ref>) and (<ref>) we refer to <cit.>, whereas (<ref>) and (<ref>) are provided in <cit.> and <cit.>, respectively.




With these steps we are now in position to state the rates of convergence associated with the Galerkin scheme (<ref>).


Assume that the hypotheses of Theorem <ref> hold and let (,,,,p)∈_(;Ω)×^2(Ω)××_(÷;Ω)×^2(Ω) and (_h,_h,_h,_h,p_h)∈^_h×^_h×^_h×^_h×^p_h be the unique solutions of the continuous and discrete problems (<ref>) and (<ref>), respectively.
Assume further that  ∈^m(Ω), ∈^m(Ω), ∈^m(Ω), ∈^m(Ω), ∈^m(Ω), ÷ ∈^m(Ω) and p∈^m(Ω), for 1≤ m ≤ k+1.
Then there exists C_rate>0, independent of h, such that 

    [                              __;Ω + __0,Ω + __0,Ω + __÷;Ω + _p_0,Ω; ≤ 
    C_rate h^m{_m,Ω +  _m,Ω + _m,Ω + _m,Ω + _m,Ω + ÷ _m,Ω + p_m,Ω}. ]

 

The result follows from Céa estimate (<ref>)
and the approximation properties (<ref>). 
	


Similarly to <cit.>, the analysis  developed in Sections <ref>-<ref> can be adapted to a formulation without the variable  (_h in the discrete problem), imposing  symmetry of  by taking  ∈ := {∈:  ∈^2(Ω) } and :={∈^2(Ω):  = ^ t}, utilizing results from <cit.> (<cit.> for the discrete problem), and adapting the strategy used in, e.g., <cit.>.




§ A POSTERIORI ERROR ESTIMATES

In this section we derive residual-based a posteriori error estimates  and  demonstrate the reliability and efficiency of the proposed estimators. Mainly due to notational convenience, we confine our analysis to the two-dimensional case. The extension to three-dimensional case should be quite straightforward (see, e.g., <cit.>).

Similarly to <cit.>, we introduce additional notation. Let ℰ_h be the set of edges of 𝒯_h, whose corresponding diameters are denoted h_E, and define

    ℰ_h(Ω) := { E∈ℰ_h :  E⊆Ω },  ℰ_h(Γ_) := { E∈ℰ_h :  E⊆Γ_ }, ℰ_h(Γ_) := { E∈ℰ_h :  E⊆Γ_ } .

On each E∈ℰ_h, we also define the unit normal vector _E:=(n_1,n_2)^  and the tangential vector _E:=(-n_2,n_1)^. However, when no confusion arises, we will simply write  and  instead of _E and _E, respectively. Also, by d/d we denote the tangential derivative. The usual jump operator [[ · ]] across internal edges are defined for piecewise continuous matrix, vector, or scalar-valued functions. For sufficiently smooth scalar ψ, vector := (v_1,v_2)^, and tensor fields := (τ_ij)_1≤ i,j≤2, we let

    (ψ):=(∂ψ/∂ x_2 , -∂ψ/∂ x_1)^ t ,  ():=∂ v_2/∂ x_1-∂ v_1/∂ x_2 ,  ()=[ (v_1)^; (v_2)^ ]()=[ (_1); (_2) ].

In addition, we denote by Π_h the Raviart–Thomas interpolator and by I_h the Clément operator (see, e.g., <cit.> for their properties). In what follows, we denote by _h the tensor version of Π_h, which is defined row-wise by Π_h and by _h the corresponding vector version of I_h which is defined componentwise by I_h.

In what follows, we will assume that the hypotheses of Theorems <ref> and <ref> are satisfied. Let _h, _h, _h, _h, p_h denote the FE solutions of (<ref>). We define the residual-based and fully computable local contributions to the error estimator Ξ_K^2, defined as the sum of Ξ_s,K^2 and Ξ_f,K^2, where Ξ_s,K and Ξ_f,K pertain to the solid (mixed elasticity) and fluid (mixed Darcy) components, respectively: 


    Ξ_s,K^2    := +_h _0,K^2+_h-_h^ t^2_0,K+h_K^2𝒞^-1_h + αdλ+2μ p_h 𝕀 - _h +_h_0,K^2
          +h_K^2𝐜𝐮𝐫𝐥(𝒞^-1_h + αdλ+2μ p_h 𝕀  +_h)_0,K^2+∑_E∈∂ K ∩ℰ_h(Ω)h_E [[(𝒞^-1_h + αdλ+2μ p_h 𝕀 +_h)]]_0,E^2
          +∑_E∈∂ K ∩ℰ_h(Γ_)h_E (
    𝒞^-1_h + αdλ+2μ p_h 𝕀 +_h) -d_/d_0,E^2 
    +∑_E∈∂ K ∩ℰ_h(Γ_)h_E _ - _h _0,E^2,
    
    Ξ_f,K^2    := -g+c_0 p_h + αdλ+2μ_h + dα^2dλ+2μ p_h - ÷_h_0,K^2
    +h_K^2κ(_h,p_h)^-1_h- ∇ p_h_0,K^2
          +h_K^2(κ(_h,p_h)^-1_h)_0,K^2+∑_E∈∂ K ∩ℰ_h(Ω)h_E [[(κ(_h,p_h)^-1_h)·]]_0,E^2 
          + ∑_E∈∂ K ∩ℰ_h(Γ_)h_E κ(_h,p_h)^-1_h· -d p_/d_0,E^2 
     + ∑_E∈∂ K ∩ℰ_h(Γ_)h_E p_ - p_h_0,E^2.




Then, we define the global estimator 

    Ξ^2:= ∑_K∈𝒯_hΞ_s,K^2 + Ξ_f,K^2.




 §.§ Reliability of the a posteriori error estimator


First we prove preliminary results that will be key in showing the reliability of the global estimator. 

There exists C_1>0, such that
    
    -_h_;Ω + -_h_0,Ω + -_h_0,Ω≤ C_1 (_1+ +(_h)_0,Ω + _h-_h^t_0,Ω),

    where 

    _1 ():=a(-_h,) + b(,(-_h,-_h)),
    
with _1 (_h)=-c(_h,p-p_h) for all _h ∈^_h, and ||_1||=sup_≠∈_(;Ω)ℛ_1()/_;Ω. 


Using the properties of bilinear forms a and b, as outlined in equations (<ref>) and (<ref>), along with the insight from <cit.>, there exists C_1>0 depending on μ,  c_a, β such that

    e__;Ω +  e__0,Ω +  e__0,Ω   ≤ C_1 sup_≠(,,) 
      ∈_(;Ω)×^2(Ω)×a( e_,) + b(,( e_, e_)) + b( e_,(,))/_;Ω + _0,Ω + _0,Ω
       ≤ C_1  (sup_≠∈_(;Ω) _1()/_;Ω+ sup_≠(,)∈^2(Ω)×b( e_,(,))/_0,Ω + _0,Ω).

Then, recalling the definitions of the bilinear form b (cf. (<ref>)), using the equation (<ref>), along with the fact that ∫_Ω_h: =1/2∫_Ω (_h-_h^ t): for ∈, the following estimate holds

    |b( e_,(,))|≤ ( ||+ _h||_0,Ω + ||_h-_h^ t||_0,Ω)(_0,Ω + _0,Ω),

and this gives the asserted inequality. 



There exists C_2>0 such that

    -_h_÷;Ω + p-p_h_0,Ω≤ C_2 (_2+ g-c_0 p_h - αdλ+2μ_h - dα^2dλ+2μ p_h + ÷_h _0,Ω + γ-_h _;Ω),

where

    _2 ():= a_,p(-_h,) +  b(,p-p_h),

satisfies _2 (_h)=0 for all _h ∈^_h, and ||_2||=
    sup_≠∈_(÷;Ω) _2()/_÷;Ω.



Similarly to Lemma <ref>, using the properties of bilinear forms a_,p and b (as outlined in equations (<ref>), (<ref>) and (<ref>)), along with the insight from <cit.>, we establish that there exists C_2>0 depending on κ_1, κ_2, C_ a, γ, γ, β such that

    e__÷;Ω +  e_p_0,Ω   ≤ C_2 sup_≠(,q)∈_(÷;Ω)×^2(Ω)  a_,p( e_,) +  b(, e_p) +  b( e_,q) -  c( e_p,q)/_÷;Ω + q_0,Ω
       ≤ C_2  (sup_≠∈_(÷;Ω) _2()/_÷;Ω+ sup_0≠ q∈^2(Ω) b( e_,q) -  c( e_p,q)/q_0,Ω).

Hence, recalling the definitions of   b,  c, adding ±  c  (_h,p), and using (<ref>), we arrive at 

    | b( e_,q) -  c( e_p,q)|≤g-c_0 p_h - αdλ+2μ_h - dα^2dλ+2μ p_h + ÷_h _0,Ωq_0,Ω + γ-_h _;Ωq_0,Ω,

and therefore, we obtain the desired result. 


Throughout the rest of this section, we provide suitable upper bounds for _1 and _2. We begin by establishing the corresponding estimates for _1, which are based on a suitable Helmholtz decomposition of the space _(;Ω) (see <cit.> for details), along with the following two technical results.

There exists a positive constant C_3, independent of h, such that for each ∈^1(Ω) there holds 

    [                                                                    |_1(-_h())| ≤  C_3  γp-p_h_0,Ω_1,Ω; + C_3(∑_K∈𝒯_h h_K^2||^-1_h + αdλ+2μ p_h 𝕀 - _h +_h||_0,K^2 + ∑_E ∈ℰ_h(Γ_)h_E||_-_h||_0,E^2 )^1/2_1,Ω. ]



From the definition of _1 (cf. (<ref>)), adding ±  c  (-_h(),p_h), and using equation (<ref>), we have

    [      _1(-_h())= H(-_h()) - c(-_h(),p) -a(_h,-_h()) 
    -b(-_h(),(_h,_h)); = ⟨ (-_h()),_⟩_Γ_ - αdλ+2μ∫_Ω (p-p_h) (-_h()) - αdλ+2μ∫_Ω p_h (-_h());                  - ∫_Ω^-1_h:(-_h()) - ∫_Ω_h:(-_h())  - ∫_Ω_h·(-_h()), ]

then, applying a local integration by parts to the last term above, using the identity
∫_E _h()·=∫_E ·, for all ∈_k(E), for all edge E of 𝒯_h, the fact that _∈^2(Γ_), and the Cauchy-Schwarz inequality, we obtain

    [ _1(-_h()) = ∑_K∈𝒯_h∫_K (-^-1_h - αdλ+2μ p_h 𝕀 + _h -_h):(-_h());  +∑_E ∈ℰ_h(Γ_)⟨ (-_h()),_ - _h⟩_E   - αdλ+2μ∫_Ω (p-p_h) (-_h());              ≤∑_K∈𝒯_h^-1_h + αdλ+2μ p_h 𝕀 - _h +_h_0,K-_h()_0,K;         +∑_E ∈ℰ_h(Γ_)_ - _h_0,E-_h()_0,E + γp-p_h_0,Ω-_h()_0,Ω, ]

with γ defined as in (<ref>). Therefore, using the approximations properties of _h (see, e.g., <cit.>) and the Cauchy–Schwarz inequality, we obtain the desired result.


Let ∈^1_Γ_(Ω):= {∈^1(Ω):  =  on Γ_} and assume that _∈^1(Γ_). Then,  there exists C_4>0, independent of h, such that 

    [                                                                           |_1( ( - _h ) )| 
     ≤ C_4 γ p- p_h||_0,Ω_1,Ω; + C_4( ∑_K∈𝒯_h h_K^2||𝐜𝐮𝐫𝐥(^-1_h + αdλ+2μ p_h 𝕀  +_h)||_0,K^2 +∑_E∈ℰ_h(Ω)h_E ||[[(^-1_h + αdλ+2μ p_h 𝕀 +_h)]]||_0,E^2;                                                 +∑_E∈ℰ_h(Γ_)h_E ||(C^-1_h + αdλ+2μ p_h 𝕀  +_h) -d_d||_0,E^2)^1/2_1,Ω. ]




Similarly to Lemma <ref>, adding ±  c  ( ( - _h ),p_h), we have

    [                                                          _1(  ( - _h )); = H( ( - _h )) - c( ( - _h ),p) -a(_h, ( - _h )) - b( ( - _h ),(_h,_h));                   = ⟨ ( ( - _h )),_⟩_Γ_ - αdλ+2μ∫_Ω (p-p_h) ( ( - _h ));     - ∫_Ω^-1_h: ( - _h ) - ∫_Ω_h: ( - _h ) - αdλ+2μ∫_Ω p_h ( ( - _h )). ]

Then, applying a local integration by parts, using that _∈^1(Γ_), the identity ⟨(-_h )   ,_⟩_Γ_= -⟨-_h , d_/d⟩_Γ_, and the Cauchy–Schwarz inequality, we obtain

    [                       _1(  ( - _h )) = -∑_K∈𝒯_h∫_K (^-1_h + αdλ+2μ p_h 𝕀 +_h):( - _h ); +∑_E∈ℰ_h(Ω)∫_E [[(^-1_h + αdλ+2μ p_h 𝕀 +_h)]]·( - _h ) - αdλ+2μ∫_Ω (p-p_h) ( ( - _h ));                                 +∑_E∈ℰ_h(Γ_)∫_E (^-1_h + αdλ+2μ p_h 𝕀 +_h -_)·( - _h );                                       ≤∑_K∈𝒯_h(^-1_h + αdλ+2μ p_h 𝕀 +_h)_0,K - _h _0,K;      +∑_E∈ℰ_h(Ω)[[(^-1_h + αdλ+2μ p_h 𝕀 +_h)]]_0,E - _h _0,E +γp-p_h_0,Ω  ( - _h )_0,Ω;                            +∑_E∈ℰ_h(Γ_) (^-1_h + αdλ+2μ p_h 𝕀 +_h) -d_d_0,E - _h _0,E. ]

As in the previous result, the approximation properties of the  Clément interpolation (see, e.g., <cit.>) in conjunction with the Cauchy–Schwarz inequality, produces the desired result.


The following lemma establishes the desired estimate for _1.


Assume that there exists a convex domain B such that Ω⊆ B and Γ_⊆∂ B, and that _∈^1(Γ_). Then, there exists a constant C_5>0, independent of h, such that

    -_h_;Ω + -_h_0,Ω + -_h_0,Ω ≤  C_5 { ∑_T∈𝒯_hΞ_s,K^2}^1/2 + C_5γp-p_h_0,Ω.





Let ∈_(;Ω).  From <cit.>, there exist ∈^1(Ω) and ∈_Γ_^1(Ω), such that

    = + _1,Ω+_1,Ω≤   C_Helm_; ,

Using that _1(_h)=c(_h,p_h-p) for all _h∈^_h, we have

    _1()=_1(-_h) + c(_h,p-p_h)   ∀ _h∈^_h.

In particular, this holds for _h defined as _h=_h+(_h ), whence

    _1()=_1(-_h)+_1((-_h )) + c(_h,p-p_h).

Hence, the proof follows from Lemmas <ref>, <ref> and <ref>, and estimate (<ref>).


The following lemma establishes the estimate for _2.


Assume that there exists a convex domain B such that Ω⊆ B and Γ_⊆∂ B, and that p_∈^1(Γ_). Then, there exists a constant C_8>0, independent of h, such that

    -_h_÷;Ω + p-p_h_0,Ω ≤  C_8 { ∑_T∈𝒯_hΞ_f,K^2}^1/2 + C_8(_h_0,Ωp-p_h_0,Ω + γ-_h _0,Ω).




It follows the steps of Lemma <ref>. From <cit.>, we have that for all ∈_(÷;Ω) there exist ∈^1(Ω) and ϕ∈_Γ_^1(Ω), such that = + ϕ and _1,Ω+ϕ_1,Ω≤ C_Helm_;.
Thus, proceeding similarly to Lemmas <ref> and <ref>, applying local integration by parts and the approximation properties of Π_h and I_h along with the Cauchy–Schwarz inequality, we can obtain the following estimates

    [                                                   |_2(-Π_h())| ≤  C_6  _h_0,Ωp-p_h_0,Ω_1,Ω; + C_6(∑_K∈𝒯_h h_K^2κ(_h,p_h)^-1_h- ∇ p_h_0,K^2 + ∑_E ∈ℰ_h(Γ_)h_E p_ - p_h_0,E^2 )^1/2_1,Ω,;                                           |_2( (ϕ - I_h ϕ) )| 
     ≤ C_7 _h_0,Ωp-p_h_0,Ωϕ_1,Ω;      + C_7( ∑_K∈𝒯_h h_K^2(κ(_h,p_h)^-1_h)_0,K^2 +∑_E∈ℰ_h(Ω)h_E [[(κ(_h,p_h)^-1_h)·]]_0,E^2;                                   +∑_E∈ℰ_h(Γ_)h_E κ(_h,p_h)^-1_h· -d p_/d_0,E^2)^1/2ϕ_1,Ω. ]

Then, noting that _2(_h)=0 for all _h∈^_h, and defining _h as _h=Π_h+(I_h ϕ), we have

    _2()=_2(-_h)=_2(-Π_h)+_2((ϕ-I_h ϕ)).

Hence, the proof follows from Lemma <ref>, estimates (<ref>), and the Helmholtz decomposition of _(÷;Ω).


Finally we state the main reliability bound for the proposed estimator.
 
Assume that the hypotheses stated in Theorem <ref> and Lemmas <ref>-<ref> are satisfied. Let (,,,,p)∈_(;Ω)×^2(Ω)××_(÷;Ω)×^2(Ω) and (_h,_h,_h,_h,p_h)∈^_h×^_h×^_h×^_h×^p_h be the unique solutions of 
(<ref>) and (<ref>), respectively.
Assume further that

    (C_5+C_8)γ + C_8 C^* (g_0,Ω + p__1/2,Γ_ + __1/2,Γ_ + _0,Ω + γ r ) ≤1/2.

Then, there exists C_rel>0, independent of h, such that 

    __;Ω + __0,Ω + __0,Ω + __÷;Ω + _p_0,Ω≤ C_rel Ξ.




It follows directly from  Lemmas <ref> and <ref>, using the fact that _h satisfies the estimate (<ref>), and the assumption (<ref>).




 §.§ Efficiency of the a posteriori error estimator



In this section we derive the efficiency estimate of the estimator defined in (<ref>). The main result of this section is stated  as follows. 

There exists C_eff>0, independent of h, such that

    C_eff Ξ ≤ -_h_;Ω + -_h_0,Ω + -_h_0,Ω + -_h_÷;Ω + p-p_h_0,Ω +  h.o.t,

where h.o.t. stands for one or several terms of higher order.

We begin with the estimates for the zero order terms appearing in the definition of Ξ_s,K (cf. (<ref>)).


For all K∈_h there holds

    +_h _0,K≲-_h_;K _h-^ t_h_0,K≲-_h_;K.




By employing the same arguments as in <cit.>, we can conclude that = -, which, together with the symmetry of ,  implies the desired result. Further details are omitted.


In order to derive the upper bounds for the remaining terms
defining the error estimator Ξ_s,K, we use results from <cit.>, inverse inequalities, and the localisation technique based on element-bubble and edge-bubble functions. The main properties that we will use can be found in <cit.>.


For all K∈_h there holds

    [           h_K𝒞^-1_h + αdλ+2μ p_h 𝕀 - _h +_h_0,K; ≲ h_K(-_h_;K  + -_h_0,K + p-p_h_0,K) + -_h_0,K. ]




It follows from an application of <cit.> with = 𝒞^-1_h + αdλ+2μ p_h 𝕀 - _h +_h, using that ^-1 + αdλ+2μ p 𝕀 =  - and 
<cit.>. We refer to <cit.> for further details.



For all K∈_h and E∈ℰ_h(Ω), there holds

    [           h_K𝐜𝐮𝐫𝐥(𝒞^-1_h + αdλ+2μ p_h 𝕀  +_h)_0,K≲-_h_;K + -_h_0,K + p-p_h_0,K,; h_E^1/2[[(𝒞^-1_h + αdλ+2μ p_h 𝕀  +_h)]]_0,E≲-_h_;ω_E + -_h_0,ω_E + p-p_h_0,ω_E, ]

where the patch of elements sharing the edge E is denoted as ω_E:=∪{K'∈𝒯_h E∈ℰ_h(K')}.



It suffices to apply <cit.> with :=^-1 + αdλ+2μ p 𝕀 + = and _h:=^-1_h + αdλ+2μ p_h 𝕀 +_h.



Assume that _ is piecewise polynomial. Then, for all E∈ℰ_h(Γ_), there holds

    [ h_E^1/2(
    𝒞^-1_h + αdλ+2μ p_h 𝕀 +_h) -d_/d_0,E≲-_h_;K_E + -_h_0,K_E + p-p_h_0,K_E,;         h_E^1/2_-_h_0,E≲  h_K_E(-_h_;K_E  + -_h_0,K_E + p-p_h_0,K_E) + -_h_0,K_E, ]

where K_E is a triangle in 𝒯_h  that contains  E on its boundary.



The first estimate follows as in <cit.>, defining  and _h as in Lemma <ref>. On the other hand, the second estimate follows from an application of the discrete trace inequality (see <cit.>), using that ^-1 + αdλ+2μ p 𝕀 =  -, and the fact that =_ on Γ_. See also  <cit.>.


A direct application of Lemmas <ref>-<ref> yields 

    ∑_K∈𝒯_hΞ_s,K≲__;Ω + __0,Ω + __0,Ω + __÷;Ω + _p_0,Ω.

Similarly, using the same arguments as in Lemmas <ref>-<ref>, along with algebraic manipulations as in Section <ref>, assuming that p_ is piecewise polynomial, together with the Lipschitz continuity of κ (cf. (<ref>)), we can bound each of the terms that appear in the estimator Ξ_f,K and obtain the following result

    ∑_K∈𝒯_hΞ_f,K≲__;Ω + __0,Ω + __0,Ω + __÷;Ω + _p_0,Ω.

We remark that the efficiency of Ξ (cf. (<ref>)) in Theorem <ref> is now a straightforward consequence of estimates (<ref>) and (<ref>). In turn, we emphasize that the resulting constant, denoted by C_eff>0 is independent of h.

For simplicity, we have assumed that _ and p_ are piecewise polynomial in the derivation of (<ref>) and (<ref>). However, similar estimates can also be obtained by assuming _ and p_ are sufficiently smooth (taking, for example, _∈^1(Γ_) and p_∈^1(Γ_), as in Lemmas <ref>-<ref>), and proceeding as in <cit.>. In such cases, higher-order terms, stemming from errors in the polynomial approximations, would appear in (<ref>) and (<ref>), accounting for the presence of h.o.t in (<ref>).


We conclude this section by noting that the a posteriori error estimation analysis developed here can be adapted to the three-dimensional case. In particular, in <cit.> and <cit.>, one can find the suitable Helmholtz decompositions for the spaces _(;Ω) and _(÷;Ω), respectively.





§ NUMERICAL RESULTS

The computational examples in this section  verify the theoretical properties (optimal convergence, conservativity,  and robustness) of the proposed schemes. The implementation has been carried out using the FE library FEniCS <cit.>. The nonlinear  systems were solved with Newton–Raphson's method with  a residual tolerance of 10^-7, and the linear systems were solved using the sparse LU factorisation of MUMPS <cit.>. 



 §.§ Optimal convergence to smooth solutions and conservativity in 2D

We first consider a simple planar problem setup with manufactured exact solution. We take the unit square domain Ω = (0,1)^2, the bottom and left segments represent Γ_ and the top and right sides are Γ_. We choose the body load f, mass source g, boundary displacement _, boundary pressure p_, as well as (not necessarily homogeneous,  but standard arguments can be used to extend the theory to the inhomogeneous case) boundary data · = φ_ and = σ_, such that the exact displacement and fluid pressure are 

    (x,y) = 1/20[ cos[3π/2(x+y)]; sin[3π/2(x-y)] ],    p(x,y) = sin(π x) sin (π y).

These exact primary variables are used to construct exact mixed variables of stress, rotation, and discharge flux. We choose the second constitutive relation for the permeability in (<ref>) and choose the following arbitrary model parameters (all adimensional)
k_0 = k_1 = c_0 = α = 0.1,   λ = μ = μ_f = 1.
These values indicate a mild permeability variation and it is expected that the nonlinear solver (in this case, Newton–Raphson) converges in only a few iterations. We construct six levels of uniform mesh refinement of the domain, on which we compute approximate solutions and the associated errors for each primal and mixed variable in their natural norms. Convergence rates are calculated as usual: 

    rate =log(e/e)[log(h/h)]^-1 ,

where e and e denote errors produced on two consecutive meshes of sizes h and h, respectively. Table <ref> reports on this error history focusing on the methods defined by the PEERS_k family with k=0 and k=1, showing a O(h^k+1) convergence for all unknowns as expected from the theoretical error bound of Theorem <ref> (except for the rotation approximation that shows a slight superconvergence for the lowest-order case and only in 2D – a well-known phenomenon associated with PEERS_k elements). With the purpose of illustrating the character of the chosen manufactured solution and the parameter regime, we show sample discrete solutions in Figure <ref>. 









We also exemplify the momentum and mass conservativity of the formulation. To do so we 





represent the loss of momentum and mass as 

    _h := 𝒫_h[ (_h) + f]_ℓ^∞,     _h := 𝒫_h[(c_0+ dα^2/dλ + 2μ) p_h + α/dλ + 2μ_h + ÷(_h) + g]_ℓ^∞,

where 𝒫_h: L^2(Ω)→P_k(𝒯_h) is the scalar version of 𝒫_h. 
They are computed at each refinement level and tabulated in Table <ref> together with the total error := e() + e() + e() + e() + e(p), and its experimental convergence rate. We report on the nonlinear iteration count as well. The expected optimal convergence of the total error, and the announced local conservativity are confirmed. We also note that, at least for this parameter regime, for all the refinements and polynomial degrees  the nonlinear solver takes  three iterations to get a residual below the tolerance. In the last column of the same table we report on the efficiency of the global a posteriori error estimator designed in Section <ref> (Ξ) = /Ξ, which – in this case of smooth solutions – is asymptotically constant (approximately 0.98 for k=0 and 1.52 for k=1).








 §.§ Convergence  in 3D using physically relevant  parameters
 Next we investigate the behaviour of the proposed numerical methods in a 3D setting and taking model parameters more closely related to applications in tissue poroelasticity. We still use manufactured solutions to assess the accuracy of the formulation, but take an exact  displacement that satisfies ÷→ 0 as λ→∞. The domain is the 3D box Ω=(0,L)×(0,L)× (0,2L) with L=0.01 m,  and mixed boundary conditions were taken analogously as before, separating the domain boundary between Γ_ defined as the planes x=0, y=0 and z=0,  and Γ_ as the remainder of the boundary. The manufactured displacement and pressure head are 

    (x,y,z) = L/4[   sin(x/L)cos(y/L)sin(z/(2L)) + x^2/λ; -2cos(x/L)sin(y/L)cos(z/(2L)) + y^2/λ; 2cos(x/L)cos(y/L)sin(z/(2L)) - 2z^2/λ ],    p(x,y,z) = sin(x/L)cos(y/L)sin(z/(2L)).







First we set again the model parameters to mild values λ = μ = c_0=k_0=α=μ_f =1, k_1=k_2=0.1, and we compare them against the following values (from, e.g., <cit.>) 

    k_0 = 2.28× 10^-11 m^3,   k_1 = 5×10^-12 m^3,   λ =1.44×10^6 Pa,  μ = 9.18×10^3 Pa,   μ_f = 7.5×10^-4 Pa·s,

and c_0 = 0, α = 0.99. Table <ref> reports on the convergence of the method. While the magnitude of the stress  errors is much higher for the second parameter regime, the discharge flux error magnitude is smaller than in the first case and the displacement, rotation, and fluid pressure errors remain roughly of the same magnitude. In any case, the table confirms that the optimal slope of the error decay is not affected by a vanishing storativity nor large Lamé constants. The numerical solutions are displayed in Figure <ref>.







 §.§ Convergence in the case of adaptive mesh refinement

We continue with a test targeting the recovery of optimal convergence through adaptive mesh refinement guided by the a posteriori error estimator proposed in Section <ref>.  

We employ the well-known adaptive mesh refinement approach  of solving, then computing the estimator, marking, and refining. Marking is done as follows <cit.>: a given K∈𝒯_h is added to the marking set ℳ_h⊂𝒯_h  whenever the local error indicator Ξ_K satisfies 

    ∑_K ∈ℳ_hΞ^2_K ≥ζ∑_K∈𝒯_hΞ_K^2,

where ζ is a user-defined bulk density. All  elements in ℳ_h are marked for refinement and also some neighbours are marked for the sake of closure. For convergence rates we use the alternative form 

    rate = -2log(e/e)[log(DoFs/ DoF)]^-1.

 
  Let us consider the non-convex rotated L-shaped domain Ω = (-1,1)^2∖ (-1,0)^2 and use manufactured displacement and fluid pressure with sharp gradients near the domain re-entrant corner (see, e.g., <cit.> for the displacement and <cit.> for the fluid pressure)

    (r,θ)  =r^χ/2μ[ -(χ + 1) cos([χ+1]θ) + (M_2 - χ - 1)M_1cos([χ-1]θ);  (χ + 1) sin([χ+1]θ) + (M_2 + χ - 1)M_1sin([χ-1]θ) ],      
    p(r,θ)  =  r^1/3sin(1/3(π/2 + θ)),

with polar coordinates r = √(x_1^2 + x_2^2), θ= arctan(x_2,x_1), χ≈ 0.54448373, M_1 = -cos([χ+1]ω)/cos([χ-1]ω), and M_2 = 2(λ + 2μ)/(μ + λ). The boundary conditions (taking as Γ_ the segments at x=±1 and y =± 1 and Γ_ the remainder of the boundary) and forcing data are constructed from these solutions and the model parameters are  λ = 10^3, μ = 10, k_0 = 1/2, μ_f =c_0=k_1= 0.1, α = 1/4, where for this test we consider a Kozeny–Carman permeability form. As in <cit.>, a sub-optimal rate of convergence is expected for the mixed elasticity sub-problem in its energy norm. Note that since the exact pressure is in ^4/3-ϵ(Ω)  for any
ϵ > 0 (cf. <cit.>), it is still regular enough to have optimal convergence. However its  gradient (and therefore also the exact discharge flux ) has a singularity located at the reentrant corner and therefore we expect an order of convergence of approximately O(h^1/3).



The numerical results of this test are reported in Table <ref>. We observe the expected sub-optimal convergence under an uniform mesh refinement while the optimal convergence in all variables is attained as the mesh is locally refined (the first three rows are very similar as most of the elements are refined in the first three steps. This can be controlled by the bulk density, here taken as  ζ = 9.5·10^-5). We also note that the individual errors are approximately of the same magnitude in the last row of each section of the table, but for the adaptive case this is achieved using approximately 5.5% of the number of degrees of freedom needed in the uniformly refined case. The last column of the table again confirms the reliability and efficiency of the a posteriori error estimator. Note that for this case we compute the divergence part of the error norm in the stress and fluxes as projections of the momentum and mass residuals onto the displacement and pressure discrete spaces, respectively. We plot in Figure <ref> the approximate displacement and pressure as well as sample triangulations obtained after a few adaptive refinement steps that confirm the expected agglomeration of vertices near the reentrant corner.  









 §.§ Adaptive computation of cross-sectional flow and deformation  in a soft tissue specimen
 Finally, we apply the proposed methods to  simulate the localisation of stress, deformation, and flow patterns in a multi-layer cross-section of cervical spinal cord. We follow the setup in <cit.>.  The geometry and unstructured mesh have been generated using GMSH <cit.> from the  images in <cit.>. The heterogeneous porous material consists of white and grey matter  surrounded by the pia mater (a thin layer, also considered poroelastic. See Figure <ref>, top two left panels). All components are assumed fully saturated with cerebospinal fluid. The transversal cross-section  has 1.3 cm in maximal diameter and  the indentation region is  a curved subset of the anterior part of the pia mater (a sub-boundary of length 0.4 cm). Boundary conditions are of mixed load-traction type, but slightly different than the ones analysed in the previous section. We conduct an indentation test applying a traction = (0,-P)^- t, with P a constant solid pressure of 950 dyne/cm^2. The posterior part of the pia mater  acts as a rigid posterior support where we prescribe zero displacement. The remainder of the boundary  of the pia mater  is   stress-free. For the fluid phase we impose a constant inflow pressure of cerebrospinal fluid of 1.1 dyne/cm^2 and zero outflow pressure at the stress-free sub-boundary, as well as zero normal discharge flux at the posterior support.  For the three different layers of the domain we use the following values for   Young modulus, Poisson ratio, and lower bound for permeability (some values from <cit.>) 

E^pia = 23'000 dyne/cm^2,  ν^pia = 0.3, E^white = 8'400 dyne/cm^2,  ν^white = 0.479,  E^grey = 16'000 dyne/cm^2,  ν^grey = 0.49,  k_0^grey = 1.4·10^-9 dyne/cm^2,  k_0^white = 1.4·10^-6 dyne/cm^2, k_0^pia = 3.9·10^-10 dyne/cm^2. Further, we take =0, g=0, μ_f = 70 dyne/cm^2· s (for cerebrospinal fluid at 37^∘), k_1 = 1/2k_0, α = 1/4, and c_0 = 10^-3. 


The initial and the final adapted mesh, together with samples of  solutions are shown in Figure <ref>, where we have used the mesh density parameter ζ = 5.5·10^-4. After each adaptation iteration guided by the a posteriori error indicator (<ref>), a mesh smoothing step was included. The figure indicates that most of the refinement occurs near the interface between the heterogeneous components of the porous media, and the plots also confirm a flow pattern moving slowly from top to bottom, consistently with a typical indentation test. 






siam

