Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelastic models

The parameters in the governing system of partial differential equations of multicompartmental poroelastic models typically vary over several orders of magnitude making its stable discretization and efficient solution a challenging task. In this paper, inspired by the approach recently presented by Hong and Kraus~[Parameter-robust stability of classical three-field formulation of Biot's consolidation model, ETNA (to appear)] for the Biot model, we prove the uniform stability, and design stable disretizations and parameter-robust preconditioners for flux-based formulations of multiple-network poroelastic systems. Novel parameter-matrix-dependent norms that provide the key for establishing uniform inf-sup stability of the continuous problem are introduced. As a result, the stability estimates presented here are uniform not only with respect to the Lam\'e parameter $\lambda$, but also with respect to all the other model parameters such as permeability coefficients $K_i$, storage coefficients $c_{p_i}$, network transfer coefficients $\beta_{ij}, i,j=1,\cdots,n$, the scale of the networks $n$ and the time step size $\tau$. Moreover, strongly mass conservative discretizations that meet the required conditions for parameter-robust stability are suggested and corresponding optimal error estimates proved. The transfer of the canonical (norm-equivalent) operator preconditioners from the continuous to the discrete level lays the foundation for optimal and fully robust iterative solution methods. The theoretical results are confirmed in numerical experiments that are motivated by practical applications.


Introduction
Multiple-network poroelastic theory (MPET) has been introduced into geomechanics [8] to describe mechanical deformation and fluid flow in porous media as a generalization of Biot's theory [9,10]. The deformable elastic matrix is assumed to be permeated by multiple fluid networks of pores and fissures with differing porosity and permeability.
During the last decade, MPET has acquired many important applications in medicine and biomechanics and therefore become an active area of scientific research. The biological MPET model captures flow across scales and networks in soft tissue and can be used as an embedding platform for more specific models, e.g. to describe water transport in the cerebral environment and to explore hypotheses defining the initiation and progression of both acute and chronic hydrocephalus [50]. In [52,51] multicompartmental poroelastic models have been proposed to study the effects of obstructing cerebrospinal fluid (CSF) transport within an anatomically accurate cerebral environment and to demonstrate the impact of aqueductal stenosis and fourth ventricle outlet obstruction (FVOO). As a consequence, the efficacy of treating such clinical conditions by surgical procedures that focus on relieving the buildup of CSF pressure in the brain's third or fourth ventricle could be explored by means of computer simulations, which can also assist in finding medical indications of oedema formation [16].
Recently, the MPET model has also been used to better understand the influence of biomechanical risk factors associated with the early stages of Alzheimer's disease (AD), the most common form of dementia [24]. Modeling transport of fluid within the brain is essential in order to discover the underlying mechanisms that are currently being investigated with regard to AD, such as the amyloid hypothesis according to which the accumulation of neurotoxic amyloid-β (Aβ) into parenchymal senile plaques or within the walls of arteries is a root cause of this disease.
Biot's and multiple-network poroelastic models are challenging from a computational point of view in that the physical parameters for different practical applications exhibit extremely large variations. For instance, permeabilities in geophysical applications typically range from 10 −9 to 10 −21 m 2 while Young's modulus is of the order of GPa and the Poisson ratio in the range 0.1 − 0.3, see [53,39,18]. Permeabilities in biological applications typically range from 10 −14 to 10 −16 m 2 . Young's modulus of soft tissues is in the order of kPa and the Poisson ratio in the range 0.3 to almost 0.5, see, e.g., [48,49]. For that reason it is important that the problem is well posed and the numerical methods for its solution are stable over the whole range of values of the physical (model) and discretization parameters.
The stability of the time discretization and space discretization by finite difference or finite volume methods have been studied in [5,23,22,43] and will not be addressed here. Instead we focus on the issue of uniform inf-sup stable finite element discretizations of the static multiple-network poroelastic problem. It is well known that the well-posedness analysis of saddle-point problems in their weak formulation, apart from the boundedness and definiteness of the underlying bilinear form, relies on a stability estimate that is often referred to as Ladyzenskaja-Babuska-Brezzi (LBB) condition [11,19]. The LBB condition, see [6,14], is also crucial in the analysis of stable discretizations and in the derivation of a priori error etsimates for mixed problems. Inf-sup stability for the Darcy problem as well as for the Stokes and linear elasticity problems have been established under rather general conditions and various stable mixed discretizations of either of these problems have been proposed over the years, see, e.g. [11] and the references therein.
Biot's model of poroelasticity combines these equations and the parameter-robust stability of its classical three-field formulation has been established only recently in [34]. Alternative formulations that can be proven to be stable include a two-field formulation for the displacement and the pore pressure [12,1] and a new three-field formulation based on introducing the total pressure as a weighted sum of fluid and solid pressure as the third unknown besides the displacement and fluid pressure [44,39]. Contrary to this new three-field formulation as analyzed in [39], the classic three-field formulation of Biot's consolidation model considered in [34] builds on Darcy's law in order to guarantee fluid mass conservation, a property that the discrete models studied in this paper maintain. Aside from two-and three-field formulations, a four-field formulation has been considered for the Biot model in which the stress tensor is kept as a variable in the system, see [38]. The error analysis in the latter work is robust with respect to the Lamé parameter λ, but not uniform with respect to the other model parameters such as K. Another formulation for Biot's model has recently been proposed and analyzed in [7]. The authors use mixed methods based on the Hellinger-Reissner variational principle for the elasticity part of the system, and impose weakly the symmetry of the stress tensor σ, resulting in a saddle point problem for σ, u, p, and a Lagrange multiplier. They prove the parameter-robust stability of the resulting four-field formulation.
The first attempt to design parameter-robust discretizations and analyze their stability for the MPET model is presented in [37]. Motivated by [44,39], the authors of [37] propose a mixed finite element formulation based on introducing an additional total pressure variable. Utilizing energy estimates for the solutions of the continuous problem and a priori error estimates for a family of compatible semidiscretizations, they show that the formulation is robust in the limits of incompressibility, vanishing storage coefficients, and vanishing transfer between networks. The robustness with respect to the permeability coefficients remains an open question in [37].
There are various discretizations for the classic three-field formulation of Biot's model that meet the conditions for the proof of full parameter-robust stability that has been presented in [34]. In general, whenever a discretization is based on a Stokes-stable pair of finite element spaces for the displacement and pressure and a Poisson-stable pair of finite element spaces for the flux and pressure unknowns, it is possible to define a parameter-dependent norm (which in general is not uniquely determined) such that the constant in the inf-sup condition for the Biot problem does not depend on any of the model or disretization parameters. For example, the triplets CR l /RT l−1 /P l−1 (l = 1, 2) together with the stabilization techniques suggested in [25,31], see also [21], or the triplets P 2 /RT 0 /P 0 (in 2D) and P stab 2 /RT 0 /P 0 (in 3D), or P 2 /RT 1 /P 1 , or the stabilized discretization that has recently been advocated in [45], or the finite element methods proposed in [36] would qualify for such parameter-robustness. However, the above-mentioned finite element methods do not have the property of strong mass conservation in the sense of satisfying the mass balance equation pointwise and therefore locally and globally on a discrete level.
A priori error estimates for the continuous-in-time scheme and the discontinuous Galerkin (DG) spatial discretization (similar to [34]) have been presented in [32] for the Biot model. Inspired by the approach proposed in [34] in context of the static Biot problem, we make use of the DG technology in the present work for solving the MPET system by introducing novel parameter-matrix-dependent norms.
The aim of this work is to establish the results regarding the parameter-robust stability of the weak formulation of the continuous problem as well as the stability of strongly mass conservative discretizations, corresponding error estimates and parameter-robust preconditioners for the multiple-network (MPET) model. The presented stability results and error estimates and preconditioners are independent of all model and discretization parameters including the Lamé parameter λ, permeability coefficients K i , arbitrary small or even vanishing storage coefficients c pi , network transfer coefficients β ij , i, j = 1, · · · , n, the scale of the networks n, the time step size τ and mesh size h. To our knowledge, these are the first fully parameterrobust stability results for the MPET model in a flux-based formulation.
The paper is organized as follows. In Section 2 the multiple-network poroelastic model is stated in a flux-based formulation, which can be considered as an extension of the classical three-field formulation considered in [34]. The governing partial differential equations are then rescaled and the static boundaryvalue problem resulting from semi-discretization in time by the implicit Euler method is presented in its weak formulation in the beginning of Section 3. The proofs of the uniform boundedness and the parameter-robust inf-sup stability of the underlying bilinear form are the main results that follow in this section. Section 4 then discusses a class of uniformly stable and strongly mass conservative mixed finite element discretizations that are based on H(div)-conforming discontinuous Galerkin approximations of the displacement field. Uniform boundedness and inf-sup stability are proved to be independent of all model and discretization parameters and the corresponding parameter-robust preconditioners are provided. Next, in Section 5, optimal parameter-robust error estimates are proved. Finally, Section 6 is devoted to the validation and illustration of the theoretical results in this work and Section 7 provides a brief conclusion. model are as follows: where σ = 2µǫ(u) + λdiv(u)I, (2a) In equation (2a), λ and µ denote the Lamé parameters defined in terms of the modulus of elasticity (Young's modulus) E and the Poisson ratio ν ∈ [0, 1/2) by The constants α i appearing in (1a) couple n pore pressures p i with the displacement variable u and are known in the literature as Biot-Willis parameters. The corresponding right hand side f describes the body force density. Each fluid flux v i is related to a specific negative pressure gradient −∇p i via Darcy's law in (1b). The tensors K i denote the hydraulic conductivities which give an indication of the general permeability of a porous medium. In (1c)u andṗ i express the time derivatives of the displacement u and the pressure variables p i . The constants c pi are referred to as the constrained specific storage coefficients and are connnected to compressibility of each fluid, for more see e.g. [47] and the references therein. The parameters β ij are the network transfer coefficients coupling the network pressures [50], hence β ij = β ji . The source terms g i in (1c) represent forced fluid extractions or injections into the medium.
It is assumed that the effective stress tensor σ satisfies Hooke's law (2a) where the effective strain tensor ǫ(u) is given by the symmetric part of the gradient of the displacement field, see (2b). Here I is used to denote the identity tensor.
The following boundary and initial conditions guarantee the well posedness of system (1): where for i = 1, . . . , n it is fulfilled Γ pi,D ∩ Γ pi,N = ∅, Γ pi,D ∪ Γ pi,N = Γ = ∂Ω and Γ u,D ∩ Γ u,N = ∅, Γ u,D ∪ Γ u,N = Γ. Initial conditions at the time t = 0 to complement the boundary conditions (3), have to satisfy (1a), and are given by The stress variable σ is eliminated from the MPET system by substituting the constitutive equation (2a) in (1a) thus obtaining the classical flux-based formulation of the MPET model.
To solve numerically the time-dependent problem, the backward Euler method is employed for time discretization resulting in the following system of time-step equations: where The unknown time-step functions u k , v k i , p k i for i = 1, . . . , n at any given time t = t k = t k−1 + τ are defined as whereas the right hand side time-step functions are . . , n. Later, the static problem (5)-(6) is considered and, for convenience, the superscript for the time-step functions is dropped, that is, u k , v k i and p k i will be denoted by u, v i and p i , respectively. The considered function spaces are as follows: • L 2 (Ω) is the space of square Lebesgue integrable functions equipped with the standard L 2 norm · ; • H 1 (Ω) d denotes the space of vector-valued H 1 -functions equipped with the norm · 1 for which u 2 1 := u 2 + ∇u 2 ; • H(div; Ω) := {v ∈ L 2 (Ω) d : div v ∈ L 2 (Ω)} with norm · div defined by v 2 div := v 2 + div v 2 . When the case Γ u,D = Γ pi,N = Γ and u D = 0, q i,N = 0 is considered, the notations U = H 1 0 (Ω) d and V i = H 0 (div, Ω), i = 1, . . . , n are used. To guarantee the uniqueness of the solution for the pressure variables p i , we set P i = L 2 0 (Ω) := {p ∈ L 2 (Ω) : Ω p dx = 0} for i = 1, . . . , n.

Stability analysis
First the parameter µ is eliminated from the system by dividing equations (5)-(6) by 2µ and making the substitutions:

Equation (5) then becomes
Next, equation (7b) is multiplied by α i τ −1 , equation (7c) is multiplied by α −1 i so that the substitutions For convenience, the "tilde" symbol is skipped and system (8) is written as: Further, we denote and make the rather general and reasonable assumptions that . . , n, and α ij ≥ 0 for i, j = 1, . . . , n.
Making use of these substitutions, without loss of generality, system (5) becomes where is the scaled operator from (6) . . , q n ) and V = V 1 × · · · × V n , P = P 1 × · · · × P n . Taking into account the boundary conditions, system (10) has the following weak formulation: Find (u; v; p) ∈ U × V × P , such that for any (w; z; q) ∈ U × V × P there holds Following [40], we first consider the following Hilbert spaces and weighted norms System (13), however, is not uniformly stable with respect to the parameters R −1 i under these norms as shown in [34]. Therefore, proper parameter-dependent norms for the spaces U , V i , P i , i = 1, . . . , n, have to be introduced that allow to establish the parameter-robust stability of the MPET model (13) for parameters in the ranges From experience, we know that the largest of the values R −1 i , i = 1, . . . , n is important to us, and we note that the term (ǫ(u), ǫ(w)) dominates in the elasticity form when λ ≪ 1. Hence, we define Again by trial and error, we find that we have to deal with the parameters in a "matrix" format. Therefore, we define the following n × n matrices j =i β ij and β ij = β ji , it is obvious that Λ 1 is symmetric positive semidefinite (SPSD). Since α pi ≥ 0, we have that Λ 2 is SPSD. Noting that R > 0, it follows that Λ 3 is symmetric positive definite (SPD). Moreover, it is obvious that Λ 4 is a rank-one matrix with eigenvalues λ i = 0, i = 1, . . . , n − 1 and λ n = n λ .
|Ω| Ω gdx and Λ g = [Λ 1 + Λ 2 , g c ] be the matrix that is obtained by augmenting Λ 1 + Λ 2 with the column g c . In general, we assume that Ω gdx = 0. When Λ 1 + Λ 2 is the zero matrix, this assumption is a "(classical) consistency condition". If Λ 1 +Λ 2 is nonzero and Ω gdx = 0, then g has to satisfy the "general consistency condition" rank(Λ 1 + Λ 2 ) = rank(Λ g ), where rank(X) denotes the rank of a matrix X. In this case, there must be p T c = (p 1,c , · · · , p n,c ) ∈ R n such that (Λ 1 + Λ 2 )p c = g c (in many applications, |Ω| Ω gdx, and thus Ω g 0 dx = 0. Then the solution (u; v; p) can be decomposed according to (Ω) and p c is a basic solution of (Λ 1 +Λ 2 )p c = g c . Therefore we only need to consider the case when Ω gdx = 0.
As we will see, it will play an important role in the definition of proper norms and the splitting (19) in our analysis. The crucial idea is that we equip the Hilbert spaces U , V , P with parameter-matrix-dependent norms · U , · V , · P induced by the following inner products: . It is easy to show that (20a)-(20c) are indeed inner products on U , V , P respectively. It should be noted that Div v, Div z and p, q are vectors and the SPD matrix Λ is used to define the norms. These novel parameter-matrix-dependent norms play a key role in the analysis of the uniform stability for the MPET model. We further point out that for n = 1, the norms defined by (20) are slightly different, but equivalent to the norms that were used in [34] to establish the parameter-robust inf-sup stability of the three-field formulation of Biot's model of consolidation.
The main result of this section is a proof of the uniform well-posedness of problem (13) under the norms induced by (20). Firstly, directly related to problem (13), we introduce the bilinear form which, in view of the definition of the matrices Λ 1 and Λ 2 , can be written in the form Then the following theorem shows the boundedness of A((·; ·; ·), (·; ·; ·)) in the norms induced by (20).
There exists a constant C b independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n and the network scale n, such that for any Proof. From the definition of the bilinear form, by using Cauchy's inequality, we obtain Then, another application of Cauchy's inequality, in view of the definition of Λ 4 , yields Before we study the uniform inf-sup condition for the MPET equations, we recall the following well known results, see, e.g. [14,11]: Furthermore, we summarize some useful properties of the matrix Λ in the following lemma.
thenΛ is SPD and for any n-dimensional vector x, we have Also, Proof. From the definitions of Λ 3 , Λ 4 , noting that Λ 3 is SPD and Λ 4 is SPSD, it is obvious thatΛ is SPD. From the definition of Λ, noting that Λ 1 and Λ 2 are SPSD, we infer the estimates From the definitions of Λ 3 , Λ 4 andΛ, we havẽ Now, using the Sherman-Morrison-Woodbury formula, we find Further, noting that where I n×n is the n-th order identity matrix, we obtain Now we can calculateΛ −1 as follows: We are ready to prove the uniform inf-sup condition for A((·; ·; ·), (·; ·; ·)) in the norms induced by (20).
Remark 5. We want to emphasize that the parameter ranges as specified in (17) are indeed relevant since the variations of the model parameters are quite large in many applications. For that reason, Theorem 2 and Theorem 3 are very important fundamental results that provide the parameter-robust stability of the model (13a)-(13c). We also point out that the matrix technique plays an interesting role for proving the uniform stability.
where Inferring from the theory presented in [42], Theorems 2 and 3 imply that the operator B defined in (36) is a uniform norm-equivalent (canonical) block-diagonal preconditioner for the operator A in (12), robust in all model and discretization parameters, i.e., κ(BA) = O(1).

Uniformly stable and strongly mass conservative discretizations
There are various discretizations that meet the requirements for the proof of full parameter-robust stability as presented in this section. They include conforming as well as nonconforming methods. In general, if is a Stokes-stable pair and V i,h /P i,h satisfy the H(div) inf-sup condition for i = 1, . . . , n, see (55), then the norm that we have proposed in Section 3 allows for the proof of full parameter-robust stability using similar arguments as in the proof of Theorem 10. To give a few examples, the triplets U h /V i,h /P i,h = CR l /RT l−1 /P l−1 (l = 1, 2) together with stabilization [25,21] results in a parameterrobust stable discretization of the MPET model if the norms are defined as in Section 3. The same is true for the conforming discretizations based on the spaces P 2 /RT 0 /P 0 (in 2D), P stab 2 /RT 0 /P 0 (in 3D), or P 2 /RT 1 /P 1 . However, the above-mentioned finite element methods do not have the property of strong mass conservation in the sense of Proposition 8 although they result in parameter-robust inf-sup stability under the norms we proposed in Section 3.
In recent years, DG methods have been developed to solve various problems [3,15,4,17,26] and some unified analysis for finite element including DG methods has recently been presented in [29,30]. In this section, motivated by the works [46,28,27], we propose discretizations of the MPET model problem (13). These discretizations preserve the divergence condition (namely equation (10c)) pointwise, which results in a strong conservation of mass, see Proposition 8. Furthermore, they are also locking-free when the Lamé parameter λ tends to ∞.

Preliminaries and notation
By T h we denote a shape-regular triangulation of mesh-size h of the domain Ω into triangles {K}. We further denote by E I h the set of all interior edges (or faces) of T h and by E B h the set of all boundary edges (or faces); we set E h = E I h ∪ E B h . For s ≥ 1, we define H s (T h ) = {φ ∈ L 2 (Ω), such that φ| K ∈ H s (K) for all K ∈ T h }.
As we consider discontinuous Galerkin (DG) discretizations, we also define some trace operators. Let e = ∂K 1 ∩ ∂K 2 be the common boundary (interface) of two subdomains K 1 and K 2 in T h , and n 1 and n 2 be unit normal vectors to e pointing to the exterior of K 1 and K 2 , respectively. For any edge (or face) e ∈ E I h and a scalar q ∈ H 1 (T h ), vector v ∈ H 1 (T h ) d and tensor τ ∈ H 1 (T h ) d×d , we define the averages and jumps where v ⊙ n = 1 2 (vn T + nv T ) is the symmetric part of the tensor product of v and n. When e ∈ E B h then the above quantities are defined as If n K is the outward unit normal to ∂K, it is easy to show that Also, for τ ∈ H 1 (Ω) d×d and for The finite element spaces we consider are denoted by U h = {u ∈ H(div; Ω) : u| K ∈ U (K), K ∈ T h ; u · n = 0 on ∂Ω}, ; v · n = 0 on ∂Ω}, i = 1, . . . , n, The discretizations that we analyze in the present context define the local spaces U (K)/V i (K)/Q i (K) via the triplets BDM l (K)/ RT l−1 (K)/P l−1 (K), or BDF M l (K)/RT l−1 (K)/P l−1 (K) for l ≥ 1. Note that for each of these choices, the important condition div We recall the following basic approximation properties of these spaces: For all K ∈ T h and for all u ∈ H s (K) d , there exists u I ∈ U (K) such that

DG discretization
We note that according to the definition of U h , the normal component of any u ∈ U h is continuous on the internal edges and vanishes on the boundary edges. Therefore, by splitting a vector u ∈ U h into its normal and tangential components u n and u t , we have implying that A direct computation shows that Similar to the continuous problem, we denote ,h × · · · × P n,h . With this notation at hand, the discretization of the variational problem (13) is given as follows: , and η is a stabilization parameter independent of parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.

Remark 7. Consider the general rescaled boundary conditions
Usually, it is assumed that the measure of Γ u,D is nonzero to guarantee the discrete Korn's inequality [13].
The standard way to incorporate the boundary conditions (46) is to modify the trial spaces according to the boundary conditions, i.e., to seek the solution in the spaces and use the test spaces given by and and η is again a stabilization parameter which is independent of λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h. If Γ u,D = Γ p,N = Γ and u D = 0, q N = 0, then (47) reduces to (44) which will be analyzed in the remainder of this paper. If the measure of Γ u,N is nonzero, then the analysis is similar. If Γ u,D = Γ and the measure of any Γ pi,D , i = 1, · · · , n, is nonzero, then one has to modify the norms according to Remark 3.1 in [34]. This part of the analysis is left as future work.
where Q i,h denotes the L 2 -projection on P i,h .
For any u ∈ U h , we introduce the mesh dependent norms: Next, for u ∈ U h , we define the "DG"-norm and, finally, the mesh-dependent norm · U h by We now summarize several results on well-posedness and approximation properties of the DG formulation, see, e.g. [28,27]: • From the discrete version of Korn's inequality we have that the norms · DG , · h , and · 1,h are equivalent on U h , namely, • The bilinear form a h (·, ·), introduced in (45) is continuous and we have • For our choice of the finite element spaces U h , V h and P h we have the following inf-sup conditions, see, e.g. [46]: where β sd and β dd are positive constant independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.
• We also have that a h (·, ·) is coercive, and the proof of this fact parallels the proofs of similar results: where α a is a positive constant independent of parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.

Related to the discrete problem (44a)-(44c) we introduce the bilinear form
In view of the definitions of the norms · U h , · V and · P , the boundedness of the bilinear form A h ((u h ; v h ; p h ), (w h ; z h ; q h )) is obvious, i.e., the following theorem holds.
Theorem 9. There exists a constant C bd independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h, such that for any (u We come to our second main result.
Theorem 10. There exits a constant β 0 > 0 independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h, such that The proof of this theorem can be obtained by following the proof of Theorem 3 and using the technique shown in [34].
From the above theorem, we get the following stability estimate.
be the solution of (44a)-(44c), then we have the estimate q h P and C 2 is a constant independent of λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.

Remark 12. Define
where Then due to the theory presented in [42], Theorems 9 and 10 imply that the norm-equivalent (canonical) block-diagonal preconditioner B h for the operator induced by the bilinear form (57) is uniform with respect to variation of the model and dicretization parameters. This means that the condition number κ(B h A h ) is uniformly bounded with respect to the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n in the ranges specified in (17), the network scale n and the mesh size h.
To apply the preconditioner B h , one has to solve an elasticity system discretized by an H(div)-conforming discontinuous Galerkin method [28] and n coupled elliptic H(div) problems discretized by RT elements. In the lowest order case and for n = 1, optimal solvers for this task have been proposed in [35].

Error estimates
In this section, we derive the error estimates that follow from the results presented in Section 4. Let Π div B : H 1 (Ω) d → U h be the canonical interpolation operator. We also denote the L 2 -projection on P i,h by Q i,h . The following Lemma, see [28], summarizes some of the properties of Π div B and Q i,h needed for our proof.
Theorem 13. Let (u; v; p) be the solution of (13) and (u h ; v h ; p h ) be the solution of (44a)-(44c). Then the error estimates and hold, where C e,u , C e,p are constants independent of λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.
Proof. Subtracting (44a)-(44c) from (13a)-(13c) and noting the consistency of a h (·, ·), we have that for Let by the stability result (58) for the discrete problem (44a)-(44c), we obtain Using the boundedness of a h (·, ·), the second inequality in Lemma 4, the triangle inequality and noting that v I is arbitrary and (Div and Remark 14. From the above theorem, we can see that the discretizations are locking-free.

Numerical Experiments
The following numerical experiments are for three widely applied MPET models, namely the one-network, two-network and four-network models. We suppose that the domain Ω is the unit square in R 2 and during the discretization it has been partitioned as bisections of 2N 2 triangles with mesh size h = 1/N . To discretize the pressure variables we use discontinuous piecewise constant elements, the fluxes are discretized employing the lowest-order Raviart-Thomas space and the displacement we approximate with the Brezzi-Douglas-Marini elements of lowest order. All the numerical tests included in this section have been carried out in FEniCS, [2,41]. The aim of these experiments is: (i) to validate the convergence of the error estimates in the derived parameter-dependent norms; (ii) to test the robustness of the proposed block-diagonal preconditioners by using it within the MinRes algorithm.
Then the exact solution of system (10) with boundary conditions u| ∂Ω = 0, v · n| ∂Ω = 0 is given by u = ∂φ1 ∂y , − ∂φ1 ∂x , p = φ 2 − 1, v = −R 1 ∇p and p ∈ L 2 0 (Ω). We performed experiments with different sets of input parameters. In Tables 1-3 we report the error of the numerical solution in the introduced parameter-dependent norms · P , · V , · U h . Additionally, we list the number of MinRes iterations n it and average residual convergence factor with the proposed blockdiagonal preconditioner where the stopping criterion is residual reduction by 10 8 in the norm induced by the preconditioner. The robustness of the method is validated with respect to variation of the parameters λ, R −1 1 , α p1 , as introduced in (10), and the discretization parameter h. As can be seen from Tables 1-3 the error in the considered parameter-dependent norms decreases by a factor 2 when decreasing the mesh size by the same factor independently of the model parameters. The results in Table 4 suggest that the number of MinRes iterations required to achieve a prescribed solution accuracy is bounded by a constant independent of λ, R −1 1 , α p1 and h while the average residual reduction factor always remains smaller than 0.70. Note that in this table the authors have tried to present the most unfavourable setting of input parameters in order to stress test the proposed method.

The two-network model
The governing partial differential equations of the Biot-Barenblatt model in which the flux-based MPET system involves two pressures and two fluxes are given by We consider here the cantilever bracket benchmark problem proposed by the National Agency for Finite Element Methods and Standards in [20] with f = 0, g 1 = 0 and g 2 = 0. The boundary of the domain Ω = [0, 1] 2 is split into Γ 1 , Γ 2 , Γ 3 and Γ 4 denoting the bottom, right, top and left boundaries respectively, and the boundary conditions u = 0 on Γ 4 , (σ − p 1 I − p 2 I)n = (0, 0) T on Γ 1 ∪ Γ 2 , (σ − p 1 I − p 2 I)n = (0, −1) T on Γ 3 , p 1 = 2 on Γ p 2 = 20 on Γ are imposed.
The base values of the model parameters are taken from [33] and are presented in Table 5. The computed numerical results in Table 6 show robust behaviour with respect to mesh refinements and variation of the parameters including high contrasts of the hydraulic conductivities. The parameter K 2 has been varied over a wider range than K 1 as it appeared to be the more interesting case.
Here, it is important to note that the authors have attempted to present again the least optimal choice of parameters for testing their implementation.

Conclusions
In this paper, motivated by the approach recently presented by Hong and Kraus [Parameter-robust stability of classical three-field formulation of Biot's consolidation model, ETNA (to appear)] for the Biot model, we establish the uniform stability, design stable disretizations and a parameter-robust preconditioners for flux-     216 Nm −2 c p1 = c p2 = c p3 = c p4 4.5 · 10 −10 m 2 N −1 α 1 = α 2 = α 3 = α 4 0.99 β 12 = β 24 1.5 · 10 −19 m 2 N −1 s −1 β 23 2.0 · 10 −19 m 2 N −1 s −1 β 34 1.0 · 10 −13 m 2 N −1 s −1 K 1 = K 2 = K 4 = K (1.0 · 10 −10 )/(2.67 · 10 −3 ) m 2 /Nsm −2 K 3 (1.4 · 10 −14 )/(8.9 · 10 −4 ) m 2 /Nsm −2 based formulations of multiple-network poroelastic systems. Novel proper parameter-matrix-dependent norms that provide the key for establishing uniform inf-sup stability of the continuous problems are introduced. The stability results that could be obtained using the presented matrix technique are uniform not only with respect to the Lamé parameter λ but also with respect to all the other model parameters such as small or large permeability coefficients K i , arbitrary small or even vanishing storage coefficients c pi , arbitrary small or even vanishing network transfer coefficients β ij , i, j = 1, · · · , n, the scale of the networks n, and the time step size τ . Moreover, strongly mass conservative and uniformly stable discretizations are proposed and corresponding uniform and optimal error estimates proved which are also independent of the Lamé parameter λ, the permeability coefficients K i , the storage coefficients c pi , the network transfer coefficients β ij , i, j = 1, · · · , n, the scale of the networks n, the time step size τ and the mesh size h. The transfer of the canonical (norm-equivalent) operator preconditioners from the continuous to the discrete level lays the foundation for optimal and fully robust iterative solution methods. Numerical experiments that are motivated by practical applications are presented confirming both the uniform and optimal convergence of the proposed finite element methods and the uniform robustness of the norm-equivalent preconditioners.