A mixed finite element method for nearly incompressible multiple-network poroelasticity

In this paper, we present and analyze a new mixed finite element formulation of a general family of quasi-static multiple-network poroelasticity (MPET) equations. The MPET equations describe flow and deformation in an elastic porous medium that is permeated by multiple fluid networks of differing characteristics. As such, the MPET equations represent a generalization of Biot's equations, and numerical discretizations of the MPET equations face similar challenges. Here, we focus on the nearly incompressible case for which standard mixed finite element discretizations of the MPET equations perform poorly. Instead, we propose a new mixed finite element formulation based on introducing an additional total pressure variable. By presenting energy estimates for the continuous solutions and a priori error estimates for a family of compatible semi-discretizations, we show that this formulation is robust in the limits of incompressibility, vanishing storage coefficients, and vanishing transfer between networks. These theoretical results are corroborated by numerical experiments. Our primary interest in the MPET equations stems from the use of these equations in modelling interactions between biological fluids and tissues in physiological settings. So, we additionally present physiologically realistic numerical results for blood and tissue fluid flow interactions in the human brain.

1. Introduction.In this paper, we consider a family of quasi-static multiplenetwork poroelasticity (MPET 1 ) equations reading as follows: for a given number of networks A ∈ N, find the displacement u and the network pressures p j for j = 1, . . ., A such that − div Cε(u) + j α j ∇ p j = f, (1.1a) c j ṗj + α j div u − div K j ∇ p j + S j = g j , 1 ≤ j ≤ A, (1.1b) where u = u(x, t) and p j = p j (x, t), 1 ≤ j ≤ A for x ∈ Ω ⊂ R d (d = 1, 2, 3) and for t ∈ [0, T ].
In our context, (1.1) originates from balance of mass and momentum in a porous, linearly elastic medium permeated by A segregated viscous fluid networks.The operators and parameters are as follows: C is the elastic stiffness tensor, each network j is associated with a Biot-Willis coefficient α j ∈ (0, 1], storage coefficient c j ≥ 0, and hydraulic conductivity tensor K j = κ j /µ j > 0 (where κ j and µ j represent the network permeability and the network fluid viscosity, respectively).In (1.1a), ∇ denotes the gradient, ε is the symmetric (row-wise) gradient, div denotes the row-wise divergence.In (1.1b), ∇ and div are the standard gradient and divergence operators, and the superposed dot denotes the time derivative.Further, f represents a body force and g j represents sources in network j for j = 1, . . ., A, while S j represents transfer terms out of network j.
In this paper, we consider the case of an isotropic stiffness tensor for which (1.2) Cε(u) = 2µε(u) + λ div uI where µ, λ are the standard non-negative Lamé parameters and I denotes the identity tensor.Moreover, we will consider the case where the transfer terms S j , quantifying the transfer out of network j into the other fluid networks, are proportional to pressure differences between the networks.More precisely, we assume that S j takes the form: S j = S j (p 1 , . . ., p A ) = A i=1 ξ j←i (p j − p i ), (1.3) where ξ j←i are non-negative transfer coefficients for i, j = 1, . . ., A. We will also assume that these transfer coefficients are symmetric in the sense that ξ j←i = ξ i←j , and note that ξ j←j is arbitrary.
The MPET equations have an abundance of both geophysical and biological applications.In the case A = 1, (1.1) reduces to the well-known quasi-static Biot equations.While the Biot equations have been studied extensively, see e.g.[32,25,29,2,28,22,38]; to the best of our knowledge, the general multiple-network poroelasticity equations have received much less attention, especially from the numerical perspective.The case A = 2 is known as the Barenblatt-Biot model, and we note that Showalter and Momken [33] present an existence analysis for this model, while Nordbotten and co-authors [27] present an a posteriori error analysis for an approximation of a static Barenblatt-Biot system.
Our interest in the multiple-network poroelasticity equations primarily stems from the use of these equations in modelling interactions between biological fluids and tissue in physiological settings.As one example, Tully and Ventikos [35] considers (1.1) with four different networks (A = 4) to model fluid flows, network pressures and elastic displacement in brain tissue.The fluid networks represent the arteries, the arterioles/capillaries, the veins and the interstitial fluid-filled extracellular space, each network with e.g. a different permeability κ j and different transfer coefficients ξ j←i .
A particularly important motivation for the current work is the recently proposed theory of the glymphatic system which describes a new mechanism for waste clearance in the human brain [18,19,1].This mechanism is proposed to take the form of a convective flow of water-like fluid through (a) spaces surrounding the cerebral vasculature (paravascular spaces) and (b) through the extracellular spaces, driven by a hydrostatic pressure gradient between the arterial and venous compartments.Compared to diffusion only, such a convective flow would lead to enhanced transport of solutes through the brain parenchyma and, in particular, contribute to clearance of metabolic waste products such as amyloid beta.The accumulation of amyloid beta frequently seen in patients with Alzheimer's disease is as such seen as a malfunction of the glymphatic system.In this context, the original system of [35] represents a macroscopic model of interaction between the different fluid networks in the brain.
Discretization of Biot's equations is known to be challenging, in particular because of so-called poroelastic locking.Poroelastic locking has two main characteristics: 1) underestimation of the solid deformation if the material is close to being incompressible and 2) nonphysical pressure oscillations, in particular in the areas close to jumps in the permeabilities or to the boundary.Several recent (and not so recent) studies, see e.g.[29,6,4,17,31,38], focus on a three-field formulation of Biot's model, involving the elastic displacement, fluid pressure and fluid velocity.Four-field formulations where also the elasticity equation is in mixed form, designed to provide robust numerical methods for nearly incompressible materials, have also been studied [37,20,21].
In biological tissues, any jumps in the permeability parameters are typically small in contrast to geophysical applications.The challenge in the biomedical applications is rather that the tissues in our body mostly consist of water and as such should be close to be incompressible (for short time-scales and normal physiological pressures).Therefore, it may be crucial for accurate modeling of the interaction of the different network pressures in (1.1) to allow for an elastic material that is almost incompressible and/or with (nearly) vanishing storage coefficients, i.e. for 1 λ < +∞ and 0 < c j 1 in (1.1).Standard two-field mixed finite element discretizations of the Biot model, approximating the displacement and the fluid pressure only using Stokes-stable elements, are well-known to perform poorly in the incompressible limit, see e.g.[22] and references therein.Moreover, we can easily demonstrate a suboptimal convergence rate for the corresponding standard mixed finite element discretization of the MPET equations, see Example 1.1 below.On the other hand, two-field approximations are computationally inexpensive compared to three-field approximations in the sense that only one unknown, the network pressure, is involved in each network.
To discretize (1.1), we consider a Crank-Nicolson discretization in time and a standard mixed finite element discretization in space in this example.More precisely, we approximate the displacement u using continuous piecewise quadratic vector fields (and denote the approximation by u h ) and the fluid pressures p j for j = 1, 2 using continuous piecewise linears defined relative to a uniform mesh of Ω of mesh size h.
The resulting approximation errors for u(T ) in the L 2 (Ω) and H 1 (Ω) norms are listed in Table 1 for a series of meshes generated by nested uniform refinements, together with the corresponding rates of convergence.We observe that the convergence rates are one order sub-optimal for this choice of spatial discretization.The primary objective of this paper is to propose and analyze a new variational formulation and a corresponding spatial discretization of the MPET equations that are robust with respect to a nearly incompressible poroelastic matrix; i.e. the implicit constants in the error estimates are uniformly bounded for arbitrarily large λ > 0. To this end, we introduce a formulation with one additional scalar field unknown.For the MPET equations (1.1) with potentially multiple networks, the additional computational cost is thus small.Instead of taking the "solid pressure" λ div u as a new unknown, we take the total pressure, which is defined as a weighted sum of the network pressures and the solid pressure, as the new unknown.Such a formulation has previously been shown to be advantageous in the context of parameter-robust preconditioners for the Biot model [23].Here, we focus on stability and error estimates of the total pressure formulation for the more general MPET equations.The construction of preconditioners for the MPET equations will be addressed in a forthcoming paper.
Our new theoretical results include an energy estimate for the continuous variational formulation that is robust in the relevant parameter limits, in particular, that is uniform in the Lamé parameter λ, storage coefficients c j for j = 1, . . ., A, and transfer coefficients ξ j←i for i, j = 1, . . ., A, and a robust a priori error estimate for a class of compatible semi-discretizations of the new formulation.These theoretical results are supported by numerical experiments.Finally, we also present new numerical MPET simulations modelling blood and tissue fluid interactions in a physiologically realistic human brain.
This paper is organized as follows.Section 2 presents notation and general preliminaries.In Section 3, we introduce a total-pressure-based variational formulation (3.6) for the quasi-static MPET equations (1.1), together with a robust energy estimate in Theorem 3.3.We continue in Section 4 by proposing a general class of compatible semi-discretizations (4.1) of this formulation, and estimate the a priori discretization errors in Proposition 4.1 and the semi-discrete errors for a specific choice of finite element spaces in Theorem 4.2 and Proposition 4.4.These theoretical results are corroborated by synthetic numerical convergence experiments in Section 5.In Section 6, we present a more physiologically realistic numerical experiment using a 4-network MPET model to investigate blood and tissue fluid flow in the human brain.Some conclusions and directions of future research are highlighted in Section 7.
2. Notation and preliminaries.Throughout this paper we use X Y to denote the inequality X ≤ CY with a generic constant C > 0 which is independent of mesh sizes.If needed, we will write C explicitly in inequalities but it can vary across expressions.

Sobolev spaces.
Let Ω be a bounded polyhedral domain in R d (d = 1, 2, or 3) with boundary ∂Ω.We let L 2 (Ω) be the set of square-integrable real-valued functions on Ω.The inner product of L 2 (Ω) and the induced norm are denoted by •, • and • , respectively.For a finite-dimensional inner product space X, typically X = R d , let L 2 (Ω; X) be the space of X-valued functions such that each component is in L 2 (Ω).The inner product of L 2 (Ω; X) is naturally defined by the inner product of X and L 2 (Ω), so we use the same notation •, • and • to denote the inner product and norm on L 2 (Ω; X).For a non-negative real-valued function on Ω (or symmetric positive semi-definite tensor-valued function on Ω) w, we also introduce the short-hand notations noting that the latter is a norm only when w is strictly positive a.e. on Ω (or is positive definite a.e. on Ω).For a non-negative integer m, H m (Ω) denotes the standard Sobolev spaces of real-valued functions based on the L 2 -norm, and H m (Ω; X) is defined similarly based on L 2 (Ω; X).To avoid confusion with the weighted L 2 -norms cf.(2.1) we use • H m to denote the H m -norm (both for H m (Ω) and H m (Ω; X)).For m ≥ 1, we use H m 0,Γ (Ω) to denote the subspace of H m (Ω) with vanishing trace on Γ ⊂ ∂Ω, and H m 0,Γ (Ω; X) is defined similarly [14].For Γ = ∂Ω, we write H m 0 (Ω) and analogously H m 0 (Ω; X). 2.2.Spaces involving time.We will consider an interval [0, T ], T > 0. For a reflexive Banach space X , let C 0 ([0, T ]; X ) denote the set of functions f : [0, T ] → X that are continuous in t ∈ [0, T ].For an integer m ≥ 1, we define where ∂ i f /∂t i is the i-th time derivative in the sense of the Fréchet derivative in X (see e.g.[39]).
For a function f : [0, T ] → X , we define the space-time norm We define the space-time Sobolev spaces W k,r ([0, T ]; X ) for a non-negative integer k and 1 ≤ r ≤ ∞ as the closure of 3. Finite element spaces.Let T h be an admissible, conforming, simplicial tessellation of the domain Ω.For any integer k ≥ 1, we let P k (T h ) denote the space of continuous piecewise polynomials of order k defined relative to T h , and P d k (T h ) as the space of d-tuples with components in P k .We will typically omit the reference to T h when context allows.We let Pk denote the restriction of these piecewise polynomial spaces to conform with given essential homogeneous boundary conditions.

Parameter values.
Based on physical considerations and typical applications, we will make the following assumptions on the material parameter values.First, we assume that the Biot-Willis coefficients α j ∈ (0, 1], j = 1, . . ., A, and the storage coefficients c j > 0 are constant in time for j = 1, . . ., A. In the analysis, we will pay particular attention to robustness of estimates with respect to arbitrarily large λ and arbitrarily small (but not vanishing) c j 's.We also comment on the case c j = 0 in Remark 4.3.
We will assume that the hydraulic conductivities K j are constant in time, but possibly spatially-varying and that these satisfy standard ellipticity constraints: i.e. there exist positive constants K − j and K + j such that We assume that the transfer coefficients ξ j←i are constant in time and non-negative: i.e. ξ j←i (x) ≥ 0 for 1 ≤ i, j ≤ A, x ∈ Ω.
2.5.Boundary conditions.We will consider (1.1) augmented by the following standard boundary conditions.First, we assume that the boundary decomposes in two parts: where |Γ| is the Lebesgue measure of Γ.We use n to denote the outward unit normal vector field on ∂Ω.Relative to this partition, we consider the homogeneous boundary conditions p j = 0 on ∂Ω for j = 1, . . ., A. (2.2c) The subsequent formulation and analysis can easily be extended to cover inhomogeneous and other types of boundary conditions.

Key inequalities.
For the space V = H 1 0,Γ D (Ω), Korn's inequality [9, p. 288] holds; i.e. there exists a constant C > 0 depending only on Ω and Γ D such that Furthermore, for the combination of spaces V and Q 0 = L 2 (Ω), the following (continuous Stokes) inf-sup condition holds: there exists a constant C > 0 depending only on Ω and Γ D such that sup Our discretization schemes will also satisfy corresponding discrete versions of Korn's inequality and the inf-sup condition with constants independent of the discretization.2.7.Initial conditions.The MPET equations (1.1) must also be complemented by appropriate initial conditions.In particular, in agreement with the assumption that c j > 0 for j = 1, . . ., A, we assume that initial conditions are given for all p j : (2.5) p j (x, 0) = p 0 j (x), x ∈ Ω, j = 1, . . ., A. Given such p 0 j , we note that we may compute u(x, 0) = u 0 (x) from (1.1a), which in particular yields a div u(x, 0) = div u 0 (x) for x ∈ Ω.In the following, we will assume that any initial conditions given are compatible in the sense described here.

A new formulation for multiple-network poroelasticity.
In this section, we introduce a new variational formulation for the quasi-static multiple-network poroelasticity equations targeting the incompressible and nearly incompressible regime.Inspired by [28,23], we introduce an additional variable, namely the total pressure.In the subsequent subsections, we present the augmented governing equations, introduce a corresponding variational formulation, and demonstrate the robustness of this formulation via an energy estimate.
3.1.Governing equations introducing the total pressure.Let u and p j for j = 1, . . ., A be solutions of (1.1) with boundary conditions given by (2.2), initial conditions given by (2.5) and recall the isotropic stiffness tensor assumption, cf.(1.2).Additionally, we now introduce the total pressure p 0 defined as Defining α 0 = 1 for the purpose of short-hand, and rearranging, we thus have that For simplicity, we denote α = (α 0 , α 1 , . . ., α A ) and p = (p 0 , p 1 , . . ., p A ), and we can thus write 2) and its time-derivative into (1.1b),we obtain an augmented system of quasi-static multiple-network poroelasticity equations: for t ∈ (0, T ], find the displacement vector field u and the pressure scalar fields p i for i = 0, . . ., A such that We note that p 0 (x, 0) can be computed from (2.5) and (3.1).Remark 3.1.In the limit λ = ∞, the equations for the displacement u and total pressure p 0 , and the network pressures p i decouple, and (3.3) reduces to a Stokes system for (u, p 0 ) and a system of parabolic equations for p j : We next present and study a continuous variational formulation based on the total pressure formulation (3.3) of the quasi-static multiple-network poroelasticity equations.

Variational formulation.
With reference to the notation for domains and Sobolev spaces as introduced in Section 2, let 3) by test functions and integrating by parts with boundary conditions given by (2.2) and initial conditions given by (2.5) yield the following variational formulation: given compatible u 0 and p 0 j , f and g j for j = 1, . . ., A, find for j = 1, . . ., A and such that u(•, 0) = u 0 (•) and p j (•, 0) = p 0 j (•) for j = 1, . . ., A. The following lemma is a modified version of Lemma 3.1 in [21] and will be used in the energy estimates below.For the sake of completeness, we present its proof here.
Proof.It suffices to show the estimate for the smallest t such that By this assumption, X (t) = max s∈[0,T ] X (s) and X (s) < X (t) for all 0 ≤ s < t.We now consider two cases: either Dividing both sides by X (t) yields (3.8) because X (t) ≥ X (0).
On the other hand, if (3.10) is the case, then (3.7) gives and taking the square roots of both sides gives (3.8).
Theorem 3.3 below establishes a basic energy estimate for solutions of (3.6), but also for solutions with an additional right-hand side (for the sake of reuse in the a priori error estimates).Theorem 3.3 (Energy estimate for quasi-static multiple-network poroelasticity).
Proof.The result follows using standard techniques.Note that the time derivative of (3.11b) reads as Taking v = u in (3.11a), q j = p j for 1 ≤ j ≤ A in (3.11c) and q 0 = −p 0 in (3.15), summing the equations, and rearranging some constants (recalling that α 0 = 1), we obtain: By definition (1.3), and the assumption that ξ j←i = ξ i←j , it follows that (3.17) Combining (3.16) and (3.17), and pulling out the time derivatives, we find that Integrating in time from 0 to t gives using Young's inequality (with ) for any 0 > 0. Again using Young's inequality with , Poincare's inequality on Q j and the assumption of uniform positivity of K j on the last terms on the right hand side of (3.18), we have that for each j = 1, . . ., A and any j > 0: with the last inequality depending on K j .Choosing j for j = 0, 1, . . ., A appropri-ately and transferring terms thus give Finally, the Cauchy-Schwarz inequality combined with Lemma 3.2, taking , and give the desired estimate.The bound for p 0 immediately follows from an inf-sup type argument: by the choice of V and Q 0 , the inf-sup condition (see e.g.[9]), by (3.6a), and Korn's inequality, we obtain that for any t ∈ (0, T ]: holds with constant depending on µ. We remark that Theorem 3.3 gives a uniform bound on u in L ∞ (0, T ; V ), p 0 ∈ L ∞ (0, T ; Q 0 ), and p j in L 2 (0, T ; Q j ) for j = 1, . . ., A, for arbitrarily large λ and arbitrarily small c j > 0 for j = 1, . . ., A in particular.
4. Semi-discretization of multiple network poroelasticity.In this section, we present a finite element semi-discretization of the total pressure variational formulation (3.3) of the quasi-static multiple-network poroelasticity equations.We introduce both abstract compatibility assumptions (A1 and A2 below) and a specific choice of conforming, mixed finite element spaces.We end this section by an a priori error estimate for the discretization error in the abstract case, and an a priori semi-discrete error estimate for a specific family of mixed finite element spaces.

Finite element semi-discretization.
Let T h denote a conforming, shaperegular, simplicial discretization of Ω with discretization size h > 0. Relative to T h , we define finite element spaces V h ⊂ V and Q i,h ⊂ Q i for i = 0, . . ., A. We assume that V h and Q i,h , i = 0, . . ., A satisfy two compatibility assumptions (A1, A2) as follows: A1: V h × Q 0,h is a stable (in the Brezzi [11] sense) finite element pair for the Stokes equations.A2: Q j,h is an H 1 -conforming finite element space for j = 1, . . ., A. We also denote With reference to these element spaces, we define the following semi-discrete total pressure-based variational formulation of the quasi-static multiple-network poroelasticity equations: for t ∈ (0, T ], find u h (t) ∈ V h and p i,h (t) ) c j ṗj,h + α j λ −1 α • ṗh + S j,h , q j + K j ∇ p j,h , ∇ q j = g j , q j ∀ q j ∈ Q j,h , (4.1c) for j = 1, . . ., A. Here S j,h = A i=1 ξ j←i (p j,h −p i,h ) cf. (1.3) and p h = (p 0,h , . . ., p A,h ).4.2.Auxiliary interpolation operators.As a preliminary step for the a priori error analysis of the semi-discrete formulation, we introduce a set of auxiliary interpolation operators.In particular, we define interpolation operators as follows.
First, for any (u, p 0 ) ∈ V ×Q 0 , we define its interpolant (Π V h u, Π Q0 h p 0 ) ∈ V h ×Q 0,h as the unique discrete solution to the Stokes-type system of equations: The interpolant is well-defined and bounded by assumption A1 and the given boundary conditions.
Second, for j = 1, . . ., A, we define the interpolation operators Π Qj h as a weighted elliptic projection: i.e. for any p j ∈ Q j , we define its interpolant Π Qj h p j ∈ Q j,h as the unique solution of (4.3) This interpolant is well-defined and bounded by assumption A2 and the given boundary conditions.
4.3.Specific choice of finite element spaces: a family of Taylor-Hood type elements.In this paper, we will pay particular attention to one specific family of mixed finite element spaces for the total pressure-based semi-discretization of the multiple-network poroelasticity equations, namely a family of Taylor-Hood type element spaces [34,5].More precisely, we note that assumptions A1 and A2 are easily satisfied by the conforming mixed finite element space pairing: for polynomial degrees l ≥ 1 and l j ≥ 1 for j = 1, . . ., A. We will refer to the spaces (4.4) as Taylor-Hood type elements of order l and l j .The superimposed ring in (4.4) denotes the restriction of the piecewise polynomial spaces to conform to the given essential boundary conditions.For this choice of finite element spaces, in particular, for the Taylor-Hood elements of order l, the following error estimate holds for the Stokes-type interpolant defined by (4.2) (see e.g.[12,7,8]).For 1 ≤ m ≤ l + 1, if u ∈ H m+1 0,Γ D (Ω) and p 0 ∈ H m , then Moreover, the following error estimate holds for the elliptic interpolants defined by (4.3) (see e.g.[10, Chap.5]): For j = 1, . . ., A, for 1 ≤ m ≤ l j , if p j ∈ H m+1 0 , it holds that (4.6) and under the full elliptic regularity assumption of Ω, (4.7) In the next subsection, we show optimal error estimates of semi-discrete solutions assuming that both of the above estimates hold.

4.4.
Semi-discrete a priori error analysis.Assume that (u, p) is a solution of the continuous quasi-static multiple-network poroelasticity equations (3.6) and that (u h , p h ) solves the corresponding semi-discrete problem (4.1).We introduce the semidiscrete (approximation) errors (4.8) e u (t) ≡ u(t) − u h (t), e pj (t) ≡ p j (t) − p j,h (t) j = 0, . . ., A, and denote e p = (e p0 , . . ., e p A ).We also introduce the standard decomposition of the errors into interpolation (superscript I) and discretization (superscript h) errors: e pj ≡ e I pj + e h pj , e I pj ≡ p j − Π Qj h p j , e h pj ≡ Π Qj h p j − p j,h j = 0, . . ., A. (4.9b) Proposition 4.1 below provides estimates for the discretization errors that are robust with respect to c j and λ.In particular, the implicit constants in the estimates are uniformly bounded for arbitrarily large λ and arbitrarily small c j > 0 for j = 1, . . ., A. We also note that the discretization errors of u in the L ∞ (0, T ; V )-norm and p j in the L 2 (0, T ; Q j )-norms for j = 1, . . ., A converge at a higher rate than the corresponding interpolation errors, as the discretization errors are bounded essentially by the initial discretization error of u in the V -norm, by the initial discretization error of p i in the L 2 -norm for i = 0, . . ., A and by the interpolation error of p i in the L 2 (0, T ; L 2 )-norm.Proposition 4.1.Assume that (u, p) ∈ C 1 (0, T ; V ) × C 1 (0, T ; Q) solves the total pressure-based variational formulation of the MPET equations (3.6) for given f and g j for j = 1, . . ., A. Assume that V h × Q h satisfies assumptions A1-A2, that (u h , p h ) ∈ C 1 (0, T ; V h ) × C 1 (0, T ; Q h ) solves the corresponding finite element semidiscrete problem (4.1), and that the discretization errors e h u and e h p are defined by (4.9).
Then, the following estimate holds for all t ∈ (0, T ]: , with an implicit constant independent of h, T , λ, c j and ξ j←i for i, j = 1, . . ., A where S j (e p ) = A i=1 ξ j←i (e pj − e pi ) and e h p0 (t) ε(e h u (t)) 2µ .Proof.A standard subtraction of (4.1) from (3.6) gives that the errors e u and e p satisfy the error equations: ) c j ėpj + α j λ −1 α • ėp + S j (e p ), q j + K j ∇ e pj , ∇ q j = 0 ∀ q j ∈ Q j,h , (4.13c) for j = 1, . . ., A with S j (e p ) = A i=1 ξ j←i (e pj − e pi ).By the definition of the interpolation operators Π h , we obtain the reduced error representations: with constant depending on µ, thus yielding (4.12).
We now consider error estimates associated with the specific choice of Taylor-Hood type finite element spaces as introduced in Section 4.3.Theorem 4.2 below presents a complete semi-discrete error estimate for this case, and is easily extendable to other elements satisfying A1 and A2.Theorem 4.2.Assume that (u, p) and (u h , p h ) are defined as in Proposition 4.1 over Taylor-Hood type elements of order l and l j for j = 1, . . ., A as defined by (4.4), and that (e u , e p ) is defined by (4.8).Assume that (u, p) is sufficiently regular.Then the following three estimates hold for all t ∈ (0, T ] with implicit constants independent of h, T , λ, c j and ξ j←i for i, j = 1, . . ., A. First, holds with E h 0 defined in (4.11), and ṗj , p j L 2 (0,t;H l j +1 ) ≡ ṗj L 2 (0,t;H l j +1 ) + p j L 2 (0,t;H l j +1 ) .
In addition, (4.17) Proof.Let (u, p), (u h , p h ) and (e u , e p ) be as stated.By the triangle inequality, the definition of e h u , Korn's inequality, and (4.5) for any t ∈ (0, T ], we have that , with inequality constant depending on Ω and µ.Further, Proposition 4.1 gives for any t ∈ (0, T ] that , where E h 0 is defined by (4.11).Applying (4.5) and (4.7), we note that for any t ∈ (0, T ] Similarly, by (4.7) and the definition of S j , we have that (4.21) Combining the above estimates and rearranging terms yield (4.16).Turning to the pressures p j , analogously using the triangle inequality, (4.6), the Poincaré inequality, and the assumptions on K j , we have for any t ∈ (0, T ] and any j = 1, . . ., A that , where the constant in the second inequality depends on Ω and the lower bound on K j .Using Proposition 4.1 together with (4.20) and (4.21), we thus obtain the estimate given by (4.17).Finally, (4.18) follows from , and (4.12).

Numerical convergence experiments.
In this section, we present a set of numerical examples to illustrate the theoretical results presented.In particular, we examine the convergence of the numerical approximations for test cases with smooth solutions.All numerical simulations in this section and in the subsequent Section 6 were run using the FEniCS finite element software [3] (version 2018.1+), and the simulation and post-processing code is openly available [30].2: Approximation errors and convergence rates for the total pressure-based mixed finite element discretization for the smooth manufactured test case for a nearly incompressible material introduced in Example 1.1.We observe that the optimal convergence is restored for the total pressure-based scheme.This is in contrast to the sub-optimal rates observed with the standard scheme (cf.Table 1).The coarsest mesh size H corresponds to a uniform mesh constructed by dividing the unit square into 4 × 4 squares and dividing each square by a diagonal.
5.1.Convergence in the nearly incompressible case.We consider the manufactured solution test case introduced in Example 1.1.As before, we consider a series of uniform meshes of the computational domain.The coarsest mesh size H corresponds to a uniform mesh constructed by dividing the unit square into 4 × 4 squares and dividing each square by a diagonal.
We let V h ×Q h be the lowest-order Taylor-Hood-type elements, as defined by (4.4) with l = 1 and l j = 1 for j = 1, . . ., A, for the semi-discrete total pressure variational formulation (4.1).For this experiment, we used a Crank-Nicolson discretization in time with time step size ∆t = 0.125 and T = 0.5.Since the exact solutions are linear in time, we expected this choice of temporal discretization to be exact.Indeed, we tested with multiple time step sizes and found that the errors did not depend on the time step size.
We computed the approximation error of u h (T ) and p h (T ) in the L 2 and H 1 -norms.The resulting errors for u h , p 0,h , and p 1,h are presented in Table 2, together with computed convergence rates.The errors and convergence rates of p 2,h were comparable and analogous to those of p 1,h and, for this reason, not reported here.From Theorem 4.2 and Proposition 4.4, we expect second order convergence (with decreasing mesh size h) for u(T ) in the H 1 -norm, second order convergence for p 0 (T ) in the L 2 -norm, first order convergence for p j (T ) in the H 1 -norm and second order convergence for p j (T ) in the L 2 -norm (since c j > 0) for j = 1, . . ., A. The numerically computed errors are in agreement with these theoretical results.In particular, we recover the optimal convergence rates of 2 for u h in the H 1 -norm, 2 for p j in the L 2 -norm and 1 for p j in the H 1 -norm.
Additionally, we observe that we recover the optimal convergence rate of 3 for u h (T ) in the L 2 -norm for this test case.Further investigations indicate that this does not hold for general ν: with ν = 0.4, the convergence rate for u h (T ) in the L 2 -norm is reduced to between 2 and 3, cf.Table 3 and dividing each square by a diagonal.We note that the third order convergence rate for u h (T ) in the L 2 -norm observed in Table 3 is reduced to order 2 − 3 in this case with ν = 0.4.

5.2.
Convergence in the vanishing storage coefficient case.We also considered the same test case, total-pressure-based discretization, and set-up as described in Section 5.1, but now with c j = 0 for j = 1, 2. The corresponding errors are presented in Table 4.We note that we observe the same optimal convergence rates as before for this case with c j = 0.

Convergence of the discretization error. Proposition 4.1 indicates superconvergence of the discretization errors e h
u and e h pj .In particular, this result predicts that for the lowest-order Taylor-Hood-type elements, we expect to observe second order convergence for the discretization error of p j in the L 2 (0, T ; H 1 )-norm.To examine this numerically, we consider the same test case, total-pressure-based discretization, and set-up as described in Section 5.1, but now compute the error between the elliptic interpolants and the finite element approximation.The results are given in Table 5 for p 1 .The numerical results were entirely analogous for p 2 and therefore not shown.We indeed observe the second order convergence of e h pj (T ) (for j = 1, 2) in the H  4: Approximation errors and convergence rates for the total pressure-based mixed finite element discretization for the smooth manufactured test case introduced in Example 1.1 but with vanishing storage coefficients (c j = 0 for j = 1, 2).We observe the optimal convergence also for this set of parameter values.The coarsest mesh size H corresponds to a uniform mesh constructed by dividing the unit square into 4 × 4 squares and dividing each square by a diagonal.
6. Simulating fluid flow and displacement in a human brain using a 4-network model.In this section, we consider a variant of the 4-network model presented in [35] defined over a human brain mesh with physiologically inspired parameters and boundary conditions.In particular, we consider the MPET equations (1.1) with A = 4.The original 4 networks of [35] represent (1) interstitial fluid-filled extracellular spaces, (2) arteries, (3) veins and (4) capillaries.In view of recent findings [1] however, we conjecture that it may be more physiologically interesting to interpret the extracellular compartment as a paravascular network.
The computational domain is defined by Version 2 of the Colin 27 Adult Brain Atlas FEM mesh [15], in particular a coarsened version of this mesh with 99 605 cells and 29 037 vertices, and is illustrated in Figure 1 (left).The domain boundary consists of the outer surface of the brain, referred to below as the skull, and of inner convexities, referred to as the ventricles, cf. Figure 1 (right).We selected three points in the domain x 0 = (89.9,108.9, 82.3) (center), x 1 = (102.We consider the following set of boundary conditions for the system for all t ∈ (0, T ).1.57 Table 6: Material parameters used for the multiple network poroelasticity equations (1.1) with A = 4 networks for the numerical experiments in Section 6.We remark that a wide range of parameter values can be found in the literature and the ones used here represents one sample set of representative values.
≈ 133.32 Pa.We assume that the displacement is fixed on the outer boundary and prescribe a total stress on the inner boundary: where n is the outward boundary normal and s is defined as where pj for j = 1, . . ., 4 are given below.We assume that the fluid in network 1 is in direct communication with the surrounding cerebrospinal fluid, and that a cerebrospinal fluid pressure is prescribed.In particular, we assume that the cerebrospinal fluid pressure pulsates around a baseline pressure of 5 (mmHg) with a peak transmantle pressure difference magnitude of δ = 0.012 (mmHg): We assume that a pulsating arterial blood pressure is prescribed at the outer boundary, while on the inner boundaries, we assume no arterial flux: For the venous compartment, we assume that a constant pressure is prescribed at both boundaries: p 3 = 6 ≡ p3 on skull and ventricles.
Finally, for the capillary compartment, we assume no flux at both boundaries: ∇ p 4 • n = 0 on skull and ventricles.
We computed the resulting solutions using the total pressure mixed finite element formulation with the lowest order Taylor-Hood type elements (l = 1 and l j = 1 for j = 1, . . ., 4 in (4.4)), a Crank-Nicolson type discretization in time with time step ∆t = 0.0125 (s) over the time interval (0.0, 3.0) (s).The linear systems of equations were solved using a direct solver (MUMPS).For comparison, we also computed solutions with a standard mixed finite element formulation (as described and used in Example 1.1) and otherwise the same numerical set-up.
The numerical results using the total pressure formulation are presented in Figures 2 and 3.In particular, snapshots of the displacement and network pressures at peak arterial inflow in the 3rd cycle (t = 2.25 (s)) are presented in Figure 2. Plots of the displacement magnitude and network pressures in a set of points versus time are presented in Figure 3.
We also compared the solutions computed using the total pressure and standard mixed finite element formulation.Plots of the displacement magnitude in a set of points over time are presented in Figure 4. We clearly observe that the computed displacements using the two formulations differ.For instance, the displacement magnitude in the point x 0 computed using the standard formulation is less than half the magnitude computed using the total pressure formulation.We also visually compared the pressures computed using the two formulations and found only minimal differences for this test case (data not shown for the standard formulation).

Conclusions.
In this paper, we have presented a new mixed finite element formulation for the quasi-static multiple-network poroelasticity equations.Our formulation introduces a single additional scalar field unknown, the total pressure.We prove, via energy and semi-discrete a priori error estimates, that this formulation is robust in the limits of incompressibility (λ → ∞) and vanishing storage coefficients (c j → 0), in contrast to standard formulations.Finally, numerical experiments support the theoretical results.For the numerical experiments presented here, we have used direct linear solvers.In future work, we will address iterative solvers and preconditioning of the MPET equations.

) for j = 1 ,
. . ., A where g I 0 = λ −1 α • e I p and g I j = −c j ėI pj − α j λ −1 α • ėI p − S j (e I p ).Noting that e h u and e h p satisfy the assumptions of Theorem 3.3 with f = 0, β = −e I p and γ j = −c j ėI pj − S j (e I p ), the semi-discrete discretization error estimate (4.10) follows.Further, by the same techniques as used for the bound (3.14), and assumption A1 combined with (4.14a), we observe that (4.15)

Fig. 1 :
Fig.1: Left: The human brain computational mesh used in Section 6 with 99 605 cells and 29 037 vertices.View from top i.e. along the negative z-axis.The points x 0 (blue), x 1 (orange), x 2 (green) are marked with spheres.Right: The inner (ventricular) boundaries of the computational mesh.View from front i.e. along negative y-axis.

Fig. 2 :Fig. 3 :
Fig. 2: Results of numerical experiment described in Section 6 using the total pressure formulation.Plots show slices of computed quantities at t = 2.25 (s) corresponding to the peak arterial inflow in the 2nd cycle.From left to right and top to bottom: (a) displacement magnitude |u|, (b) extracellular pressure p 1 , (c) arterial blood pressure p 2 , (d) venous blood pressure p 3 and (e) capillary blood pressure p 4 .

Fig. 4 :
Fig.4: Comparison of displacements computed using the standard and total pressure formulation (cf.Section 6).Plots of displacement magnitude |u(x i , t)| versus time t, for a set of points x 0 , x 1 , x 2 (see Figure1for the location and precise coordinates of the points x i ): (a) Total-pressure mixed finite element formulation, (b) Standard mixed finite element formulation (cf.Example 1.1).The computed displacements clearly differ between the two solution methods.

Table 1 :
Approximation errors in the L 2 ( • )-and H 1 ( • H 1)-norms and associated convergence rates for a standard mixed finite element discretization for a smooth manufactured solution test case for a nearly incompressible material (Example 1.1).H corresponds to a uniform mesh constructed by dividing the unit square into 4 × 4 squares and dividing each square by a diagonal.

Table 3 :
. Displacement approximation errors and convergence rates for the total pressure-based mixed finite element discretization for the smooth manufactured test case introduced in Example 1.1 but with ν = 0.4.The coarsest mesh size H corresponds to a uniform mesh constructed by dividing the unit square into 4 × 4 squares

Table
1-norm as indicated by Proposition 4.1.

Table 5 :
Discretization errors and convergence rates for p 1 for the total pressure-based mixed finite element discretization for the smooth manufactured test case for a nearly incompressible material introduced in Example 1.1.We indeed observe the higher (second) order convergence of e h p1 (T ) in the H 1 -norm as indicated by Proposition 4.1.The coarsest mesh size H corresponds to a uniform mesh constructed by dividing the unit square into 4 × 4 squares and dividing each square by a diagonal.
All boundary pressure values are given in mmHg below, noting that 1 mmHg