A posteriori error estimation and adaptivity for multiple-network poroelasticity

The multiple-network poroelasticity (MPET) equations describe deformation and pressures in an elastic medium permeated by interacting fluid networks. In this paper, we (i) place these equations in the theoretical context of coupled elliptic-parabolic problems, (ii) use this context to derive residual-based a posteriori error estimates and indicators for fully discrete MPET solutions and (iii) evaluate the performance of these error estimators in adaptive algorithms for a set of test cases: ranging from synthetic scenarios to physiologically realistic simulations of brain mechanics.


Introduction
At the macroscale, the brain and other biological tissues can often be viewed as a poroelastic medium: an elastic structure permeated by one or more fluid networks.Such structures can be modeled via Biot's equations in the case of a single fluid network [8,9,35] or by their generalization to the equations of multiple-network poroelasticity (MPET) which describe the case of two or more interacting fluid networks [2,4,5,6,18,20,36,37,39].However, the computational expense associated with the numerical solution of these equations, over complex domains such as the human brain, is substantial.A natural question is therefore whether numerical error estimation and adaptivity can yield more accurate simulations of the MPET equations within a limited set of computational resources.
The quasi-static MPET equations read: given a domain Ω, a finite final time T > 0 and a set of J fluid networks, find the displacement field u : [0, T ] × Ω → R d and pressure fields p 1 , p 2 . . .p J : [0, T ] × Ω → R such that ∂ t (s j p j + α j div u) − div (κ j ∇ p j ) + T j = g j . (1.1b) The quantity σ (u) = 2µε (u) + λ tr (ε (u)) I in (1.1a) is the elastic stress tensor and involves the displacement u, the linearized strain tensor ε (u) = 1 2 ∇ u + ∇ u T , the d×d identity matrix I and the material Lamé coefficients µ and λ.Each one of the J fluid networks is associated with a Biot-Willis coefficient α j , a storage coefficient s j , and a hydraulic conductivity κ j .An interpretation of the Biot-Willis and storage coefficients, in the MPET context, appears in [4,Section §3].We use transfer terms T j in (1.1b) of the form T ji , T ji = γ ji (p j − p i ). (1. 2) The coefficients γ ji regulate the interplay between network i and network j and T j is the total transfer out of network j (into the other networks).The transfer term vanishes when J = 1 and (1.1) coincides with Biot's equations for a single fluid in a poroelastic medium.We also note that the fluid (Darcy) velocity v j in network j is defined by Over the last decade, several authors have studied a posteriori error estimation and adaptivity related to (1.1) in the case of J = 1; that is, for Biot's equations of poroelasticity.Depending on the application of interest, different formulations of Biot's equations have been used which introduce additional solution fields such as the Darcy velocity, the total pressure or the effective stress.In each case, a posteriori methods have been developed to facilitate adaptive refinement strategies.In [30], the authors consider the standard two-field formulation of Biot's equation in two spatial dimensions, develop an a posteriori error analysis based on H(div) reconstructions of the flux and effective stress and apply the resulting estimators to construct a time-space adaptive algorithm.A posteriori estimators have also been used to provide error estimates for the popular fixed-stress iterative solution scheme applied to the two-field formulation [22].Formulations with additional fields have also been considered for Biot's equations.The total pressure formulation [29] is a locking-free, three-field formulation, ideal for a nearly-incompressible poroelastic material.A priori estimates, and an adaptive refinement strategy, for this formulation are constructed in [21] for quadrilateral and simplicial meshes.Residual-based a posteriori error estimates have also been advanced [24] for a lowest-order discretization of the standard Darcy-flux three-field formulation which, as shown in [31], robustly preserves convergence in the presence of vanishingly small hydraulic conductivity.Finally, a four-field formulation, with symmetric stress and a Darcy velocity, of Biot's equations has been used to develop [1] a posteriori error estimates, and an adaptive refinement, based on post-processed pressure and displacement fields.
The a posteriori landscape for the more general MPET system (1.1) is considerably sparse.A posteriori error estimates for the two-field formulation of the Barenblatt-Biot equations (corresponding to the J = 2 case of (1.1)) have indeed been obtained by Nordbotten et al. [27].In general, though, there has been little work on the development of a posteriori error estimators for (1.1), for formulations with any number of fields, in the case of more than one fluid network (i.e.J > 1).However, the recent work of [15] developed an abstract framework for a posteriori error estimators for a general class of coupled elliptic-parabolic problems.
In this manuscript, our focus is three-fold.First, we rigorously place the MPET equations in the context of coupled elliptic-parabolic problems.In particular, we consider extended spaces, bilinear forms and augmentation with a semi-inner product arising from the additional transfer terms.Second, we use this context to derive specific a posteriori error estimates and error indicators for the space-time finite element discretizations of the multiple-network poroelasticity equations in general.In biomedical applications, two-field variational formulations are often used to numerically approximate the multiple-network poroelasticity equations [36,37], and we therefore focus on such here.Third, we formulate a physiological modelling and simulation-targeted adaptive strategy and evaluate this strategy on a series of test cases including a clinically-motivated simulation of brain mechanics.

Notation and preliminaries
This section provides a brief account of the notation and relevant results employed throughout the remainder of the manuscript.
2.1.Domain, boundary and meshes.It is assumed that the poroelastic domain Ω ⊂ R d with d ∈ {1, 2, 3} is a bounded, convex domain with ∂Ω Lipschitz continuous.We consider a family of mesh discretization {T h } h>0 of Ω into simplices; triangles when d = 2 and tetrahedra when d = 3.Here, h > 0 is a characteristic mesh size such as the maximum diameter over all simplices.Furthermore, we assume that each mesh T h in the family is quasi-uniform.

Material parameters.
For simplicity, we assume that all material parameters are constant (in space), and that the following (standard) bounds are satisfied: µ > 0, 2µ + λ > 0, κ j > 0, α j ∈ (0, 1], s j > 0 for j = 1, . . ., J, and γ ji = γ ij ≥ 0 for i, j = 1, . . ., J with γ ii = 0.The analysis can be extended to the case where the parameters, above, vary with sufficient regularity in space and time provided the above bounds hold uniformly.For each parameter, ξ i or ξ ij above, the minimum and maximum notation will be used throughout the manuscript; the notational extension to the case of smoothly varying parameters, on a bounded domain with compact closure, is clear.
When the context is evident in praxis the domain, Ω, is suppressed in the above expressions.Given w a positive constant, positive scalar field, or positive-definite tensor field, the symbolics w , refer to a w-weighted inner product and norm, respectively.
The Sobolev space H 1 (Ω), often abbreviated as simply H 1 , consists of those functions f ∈ L 2 whereby ∂ xj f exists, in the sense of distributions, for every j = 1, 2, . . ., d and ∂ xj f ∈ L 2 .The associated norm is given by the expression The subset H 1 0 ⊂ H 1 signifies functions with zero trace on the boundary; that is, those functions f ∈ H 1 such that f (x) = 0 for almost every x ∈ ∂Ω.In addition, given a Hilbert space X with inner product The additional decoration of the inner product, for the case of a Hilbert space X, will be omitted when the context is clear.For X any Banach space, the notation X * denotes the dual space,and we also write x * , x X ×X for the duality pairing.Accordingly, the operator norm on Unlike the inner product case, the decoration of the duality pairing bracket notation will always be made explicit and never omitted.
We also recall the canonical definition [16] of some useful time-dependent spaces whose codomain is also a given Hilbert space X.With X selected we consider a strongly measurable function We now discuss those strongly measurable functions, f : [0, T ] → X, which possess weakly differentiability in time.The space H 1 (0, T ; X) denotes the collection of functions, f ∈ L 2 (0, T ; X), such that ∂ t f exists, in the weak sense, and also resides in L 2 (0, T ; X).This is similar to the usual definition of H 1 (Ω), given above, and the norm corresponding to this space is also similar; it is given by Likewise, f ∈ C k (0, T ; X) implies that f and its first k weak derivatives in time, ∂ j t f for j = 1, 2, . . ., k, all reside in C(0, T ; X).

2.4.
Mesh elements and discrete operators.For a fixed h, the mesh T h is composed of simplices, denoted T ∈ T h , and faces (edges in 2D) e ∈ ∂T .Let Γ denote the complete set of faces of simplices T ∈ T h ; then Γ can be written as the disjoint union where e ∈ Γ int if e is an interior edge and e ∈ Γ bd if e is a boundary edge.
Let f be a scalar or vector valued function and suppose e is an interior edge e ∈ T + ∩ T − where T + and and T − are two simplices with an arbitrary but fixed choice of labeling for the pairing.We denote by n e the outward facing normal associated to T + .We use an explicit jump operator defined, for e ∈ Γ int , by where f + denotes f restricted to e ∈ T + and f − denotes f restricted to e ∈ T − .For an edge e ∈ Γ bd we have that there exists only one T + = T ∈ T h such that e ∈ ∂T and in this case we define [f ] = f + .
2.5.Boundary and initial conditions.We assume homogeneous boundary conditions for the displacement and pressures; though, as in [15], these conditions can easily be generalized [32].

Coupled elliptic-parabolic problems as a setting for poroelasticity
To consider the a posteriori error analysis of the generalized poroelasticity equations (1.1), we follow the general framework for a posteriori error analysis for coupled elliptic-parabolic problems presented by Ern and Meunier [15].In Section 3.1 below, we briefly overview this general framework and its application to Biot's equations.Next, we show that (1.1) can be addressed using this general framework, also for the case where J > 1 under appropriate assumptions on the transfer terms T m→n , in Section 3.2.Based on the general framework, Ern and Meunier derive and analyze several a posteriori error estimators.These estimators, and their corresponding extensions to the generalized poroelasticity equations, will be discussed in Section 5.

3.1.
The coupled elliptic-parabolic problem framework.The setting introduced by Ern and Meunier [15] for coupled elliptic-parabolic problems provides a natural setting also for generalized poroelasticity.The general coupled elliptic-parabolic problem reads as: find (u, p) ∈ H 1 (0, T ; V a )× H 1 (0, T ; V d ) that satisfy (for almost every t ∈ [0, T ]): The data, f and g in (1.1), are general and assumed to satisfy f ∈ H 1 (0, T ; V * a ) and g ∈ H 1 (0, T ; V * d ).The initial pressure is assumed to satisfy p(0) ∈ V d .Moreover, it is assumed that (1) V a and V d are Hilbert spaces.
(2) a : Example 3.1.Biot's equations of poroelasticity (i.e.(1.1) for J = 1 fluid networks) fit the coupled elliptic-parabolic framework with and with the standard (vector) H 1 0 -inner product and norm on V a and V d , and L 2 -inner product and norm on L a and L d .It is readily verifiable that the general conditions described above are satisfied under these choices of spaces, norms, inner products and forms [15].The existence and uniqueness of solutions to Biot's equations of poroelasticity (J = 1 in (1.1) and (3.1) with the above bilinear forms) is now a classical result [32].

3.2.
Generalized poroelasticity as a coupled elliptic-parabolic problem.In this section, we derive a variational formulation of the generalized poroelasticity equations (1.1) for the case of several fluid networks (i.e.J ≥ 1) and show how this formulation fits the general framework presented above.The extension from Biot's equations to generalized poroelasticity is natural in the sense it coincides with the original application of the general framework to Biot's equations when J = 1.Suppose that the total number of networks J is arbitrary but fixed.We define the spaces We consider data such that f ∈ H 1 (0, T ; L a ) and (g 1 , g 2 , . . ., g J ) ∈ H 1 (0, T ; L d ) with given initial network pressures determined by p(0) ∈ V d .A standard multiplication, integration and integration by parts yield the following variational formulation of (1.1): find u ∈ H 1 (0, T ; V a ) and p = (p 1 , . . ., p J ) ∈ H 1 (0, T ; V d ) such that for a.e.t ∈ (0, T ]: J j=1 ∂ t s j p j , q j + ∂ t α j div u, q j + κ j ∇ p j , ∇ q j + T j , q j = J j=1 g j , q j .(3.3b) for v ∈ V a , q = (q 1 , . . ., q J ) ∈ V d .As noted in [15], (3.3a) holds up to time t = 0 so that u 0 is determined by the initial data p(0) and initial right-hand side f (0).By labeling the forms c(p, q) = J j=1 s j p j , q j , (3.4c) we observe that the weak formulation (3.3) of (1.1) takes the form (3.1) where T j , in (3.4d), is given by (1.2).
We now show that the associated assumptions on these forms and spaces hold, beginning with properties of the form d in Lemma 3.2 below.Lemma 3.2.The form d given by (3.4d) defines an inner product over [H 1 0 (Ω)] J with associated norm where | • | T is defined by (3.7), which is such that with inequality constant depending on J, κ max , γ max and Ω.
Proof.By definition (1.2) and the assumption of symmetric transfer γ ji = γ ij ≥ 0, we have Given p, q ∈ V d , the bilinear form defined by (3.7), that is p, q T = J j=1 T j , q j , is clearly symmetric and satisfies the requirements of a (real) semi-inner product on L d × L d in the sense of [11].It follows that defines a semi-norm on L d × L d and that the corresponding Cauchy-Schwarz inequality holds.Using the triangle inequality, the definition (3.7), the bounds for γ ji and the Poincaré inequality, we have that with constant depending on γ max , J, and the domain via the Poincaré constant.Under the assumption that κ min > 0, we observe that as a result d defines an inner product and norm on [H 1 0 ] J × [H 1 0 ] J .Similarly, (3.6) holds with with constant depending on κ max in addition to γ max , J, and the domain Ω.Lemma 3.2 will be used in the subsequent sections.We next show that the choices of spaces (3.2) and forms (3.4) satisfy the abstract assumptions of the framework as overviewed in Section 3.1, and summarize this result in Lemma 3.3.
Lemma 3.3.The problem (3.3), arising from the equations of generalized poroelasticity (1.1) with material parameters as in Section 2.2, posed on the spaces (3.2) with bilinear forms defined via (3.2) is a coupled elliptic-parabolic problem and satisfy the assumptions set forth in [15].
Proof.We consider each assumption in order.These are standard results, but explicitly included here for the sake of future reference.
(1) V a and V d defined by (3.2) are clearly Hilbert spaces with natural Sobolev norms • H 1 0 .(2) a is symmetric, coercive on V a by Korn's inequality and the lower bounds on µ, 2µ + dλ, and continuous (with continuity constant depending on µ max and λ max ).d is clearly symmetric by the transfer symmetry assumption and (3.7), coercive by | • | T ≥ e0 and the assumption that κ min > 0: by applying Cauchy-Schwarz and Hölder's inequality, with constant depending on α min > 0 and the coercivity constants of a and c.
In light of Lemma 3.3, the generalized poroelasticity system (3.3) is of coupled elliptic-parabolic type and takes the form of (3.1) with bilinear forms defined by (3.4).
Corollary 3.4.The following energy estimates hold for almost every t ∈ [0, T ] Proof.The proof follows directly from Lemma 3.3, the corresponding energy estimates for coupled elliptic-parabolic systems [15, Prop.2.1] and the definition of the norms arising from the forms (3.4).Moreover, the proportionality constant in the estimates is independent of all material parameters and the number of networks.
Remark 3.2.The elliptic-parabolic MPET energy estimates, of Corollary 3.4, are similar to those of the total pressure formulation [23,Theorem 3.3] when the second Lamé coefficient, λ, is held constant in the latter.The primary difference is that [23] separates the estimates of u from that of the solid pressure, λ div u, by including the latter term into a 'total pressure' variable.This allows for ||u|| 1 to be estimated directly, in [23], regardless of the value of λ used in the definition of ||u|| a .
Remark 3.3.The conditions of Section 3.1, i.e. conditions (2)-( 5), can place restrictions on the generalized poroelastic setting.As an example, the assumption of a vanishing storage coefficient has appeared in the literature as a modeling simplification [25,40].However, the coercivity requirement of condition (4) precludes the use of a vanishing specific storage coefficient, s j in (3.4c), for any network number j = 1, 2, . . ., J. Care should be taken to ensure that any modeling simplifications produce forms that satisfy the conditions of Section 3.1 in order for the results of Corollary 3.4 and Section 4 to hold.
4. Discretization and a priori error estimates 4.1.An Euler-Galerkin discrete scheme.We now turn to an Euler-Galerkin discretization of (1.1) in the context of such discretizations of coupled elliptic-parabolic problems in general [15].We employ an implicit Euler discretization in time and conforming finite elements in space.
We consider a family of simplicial meshes {T h } h>0 with h a characteristic mesh size such as the maximal element diameter For functions and fields, we use the superscript n to refer values at time point t n .We also utilize the discrete time differential notation δ t where With this notation, the discrete problem is to seek such that for all time steps t n with n ∈ {1, 2, . . ., N }: where the spaces and forms are defined by (3.2) and (3.4).The right-hand sides, above, express the inner product of the discrete approximations f n h ∈ L a,h , to f and g n h ∈ L d,h , to g, at time t n .By Lemma 3.3 and [15, Lemma 2.1], the discrete system (4.2) is well-posed.

4.2.
A priori error estimates.Now, let V a,h and V d,h be spatial discretizations arising from continuous Lagrange elements of order k a and k d , respectively, where k d = k a − 1; this relation on relative degree results directly from the framework hypotheses [15, Section 2].Let P k (T ) denote polynomials of order k on a simplex T ∈ T h .We consider the continuous Lagrange polynomials of order k a and k a−1 defined by as the discrete spaces for the displacement and network pressures, respectively.When J = 1 this choice coincides with the previous [15] discretization considered for Biot's equations.
The general framework stipulates that three hypotheses [15, Section 2.5], restated here for completeness, should be satisfied for the discretization. ) Hypothesis 4.2.There exists a real number δ such that for every r ∈ L d , the unique solution We now state the primary result of this section.Proof.Choose s a = k a and s d = k a−1 .Then the conditions of Hypothesis 4.1 follow, as in [15], from choosing J where H k (T h ) denotes the broken Sobolev space of order k on the mesh T h .The estimate (4.5) follows, without extension, directly from classical results in approximation theory [14]; precisely as discussed in [15].Similarly, the estimate (4.6) follows from the properties of d, standard interpolation estimates [14], and the product structure of V d and V d,h .For hypothesis 4.2, we use the elliptic regularity, c.f. standard well posedness and interior regularity arguments in [16,Chp. 6], of the solution to the coupled linear diffusion-reaction equation of finding φ ∈ V d such that and Γ is a matrix composed of the transfer coefficients γ ij : Γ ii = j γ ji , and It follows from Lemma 3.2, and Γ symmetric and positive-semi definite, that the solution φ ∈ V d to the dual problem lies in H where δ = 1; exactly as in [15].Finally, with δ = 1 and the choices s a = k a and s d = k a − 1, hypothesis 4.3 also holds.
A priori estimates for the Euler-Galerkin discretization (4.2) of the generalized poroelasticity equations (1.1) then follow directly from [15,Theorem 3.1].These estimates will be used in the a posteriori analysis and are restated from [15], subject to the extended spaces and forms of (3.2) and (3.4).

Corollary 4.5 (A priori estimates for generalized poroelasticity). Let
where W a and W d are given above with V a and L d as in (3.2).It is also assumed that the initial data satisfies Setting, for simplicity, s = k a then s = k d + 1, by the selection of the discrete spaces, and we have that for each n ∈ {1, 2, . . ., N } and (4.10)

A posteriori error estimation for generalized poroelasticity
We now turn to discuss the implications to a posteriori error estimates for generalized poroelasticity as viewed through the lens of the coupled elliptic-parabolic problem framework.Our focus is to derive, apply and evaluate residual-based error estimators and indicators in the context of generalized poroelasticity.We will therefore present an explicit account of abstractly defined quantities presented in [15,Sec. 4.1], including e.g. the Galerkin residuals, applied in our context.
Similarly (p 1hτ , p 2,hτ , . . ., p J,hτ ) = p hτ is defined by p hτ (t n ) = p n h and extended linearly in time.As a result, ∂ t u hτ , ∂ t p hτ are defined for almost every t ∈ (0, T ).Define the corresponding continuous, piecewise linear in time variants of the data, f hτ and g hτ , by the same approach; i.e. f hτ (t n ) = f n h and g hτ (t n ) = g n h .Before rephrasing (4.2) using the time-interpolated variables we define piecewise constant functions in time for the pressure and right-hand side data.These are defined as π 0 p hτ = p n h and π 0 g hτ = g n h on I n = (t n−1 , t n ).Using the above notation, the discrete scheme for almost every t ∈ (0, T ) becomes (5.1b) Remark 5.1.Using the linear time interpolations defined above, such as u hτ or p hτ , we have the following identity so that the left-hand sides of and (4.2b) and (5.1b) are identical.However, as noted in [15], the interpolants of the data, f hτ and g hτ , are continuous and facilitate the definition of the continuoustime residuals (5.2) and (5.3).

Galerkin residuals.
The Galerkin residuals [15, Section 4.1] are functions of time whose co-domain lies in the dual of either V a or V d .More specifically, the residuals are continuous, piecewise-affine functions In our context, of generalized poroelasticity, the Galerkin residual G a is, given any v ∈ V a , defined by the relation Similarly, G d is, given any q = (q 1 , q 2 , . . ., q J ) ∈ V d , defined by the relation Again, we note that (5.2) and (5.3) generalize the corresponding [15, Sec.4.1] residuals for the case of single-network poroelasticity studied therein.
5.3.Data, space and time estimators.The general coupled elliptic-parabolic problem framework gives a posteriori error estimates, and in particular so-called data, space and time estimators for the discrete solutions.In our context of generalized poroelasticity, these can be expressed explicitly as follows.We have terms for the data f and g given by and the framework data, space and time estimators are defined, respectively, as where we recall that the norms are now defined according to the extended generalized poroelasticity spaces (3.2) and forms (3.4).The following a posteriori error estimate for the general MPET equations holds: Proposition 5.1.For every time t n , with n ∈ {1, 2, . . ., N }, the following inequality holds Proof.The proof follows from [15,Thm. 4.1] and the arguments of Section 3.2.
Remark 5.2.Note [15, eqn.(4.10)] that (5.6) is equivalent to The above expression will be used in Section 5.6 and follows, via the definition of p hτ and π 0 , from the calculation 5.4.Element and edge residuals.In this section we state the definition of the element and edge residuals (c.f.[15, Sec.4.1]) adapted to generalized poroelasticity.We then define from these residuals a set of a posteriori error indicators.These indicators can be used to bound the Galerkin residuals defined in Section 5.2.The a posteriori error indicators defined in this section will be used to carry out adaptive refinement for the numerical studies in Section 6.

5.4.1.
Element and edge residuals for the momentum equation.The residuals associated with the displacement are derived from the Galerkin residual (5.2).We give them explicitly here for the sake of clarity and to facilitate implementation.For v ∈ V a and at time t n , with n ∈ {1, 2, . . ., N }, we have where the notation f, g K = K f g dx denotes local integration over a simplex K ∈ T h and we have used that u n hτ = u n h and p n j,hτ = p n j,h for every n ∈ {1, 2, . . ., n}.Integrating the above by parts over each K ∈ T h gives where Γ int denotes the set of interior edges, f, g e denotes integration over the edge e and where where the additional subscript denotes the restriction K.To define the term J n uh above we use the standard notation, of (2.1), and define where e is an edge and n e is the fixed choice of outward facing normal to that edge.The corresponding time-shifted local residual and jump operators are then (5.12) To close, we note that the conditions in [15] on the jump operator, J uh,e above, are general and other choices satisfying the abstract requirements can be used if desired.

5.4.2.
Element and edge residuals for the mass conservation equation.The residuals associated with the network pressures are derived from the Galerkin residual G d (5.3).Integrating the diffusion terms by parts, over T ∈ T h , gives (5.13) In the context of the extended multiple-network poroelasticity framework the strong form of the mass conservation residual, R n ph,K ∈ L d of (5.13), has a j th component, for j ∈ {1, 2, . . ., J}, with recalling that T j is given by (1.2), and T j,h is its discrete analogue.In (5.14), we have also used that ∂ t p n j,hτ = δ t p n j,h = τ −1 m (p n j,h − p n−1 j,h ), ∂ t u n h,τ = δ t u n h , p n j,hτ = p n j,h and u n j,hτ = u n j,h .The corresponding jump term J n ph,e for e ∈ Γ int has j th component and we once more remark that other jump operators satisfying the abstract conditions in [15] can also be considered.We also have the analogous time-shifted versions of the above, δ t R n ph,K and δ t J n ph,K , just as in (5.12).

5.5.
Error indicators in space and time.We now define the element-wise error indicators; these indicators will inform the construction of the a posteriori error indicators of Section 5.6.In turn, these indicators will form the foundation of the adaptive refinement strategy of section 7. Specifically, we define element-wise indicators denoted η n u,K and η n p,K such that the following equalities hold for all v ∈ V a and q ∈ V d (5.16) First, we define the following local error indicator associated with the momentum equation: Similar to (5.12) we will use the time-shifted version of the local spatial error indicator for the momentum equation.This expression is given by where the right-hand is analogous to that of (5.17 In Section 7 we will use the above expressions to define the a posteriori error indicators informing a simple adaptive refinement strategy for the numerical simulations of Section 6.

5.6.
A posteriori error estimators.We close this section by defining the final a posteriori error estimators: (5.20) The summed term ||p n h − p n−1 h || 2 d , in η 4 above, can be expanded using the definition of (3.4d) as Finally, a bound on the MPET discretization errors in terms of the a posteriori error estimators follows: Proposition 5.2.For each time t n , n ∈ {0, 1, . . ., N }, the following inequality for the discretization error holds where E h0 (u 0 , p 0 ) and E h (f, g) are determined by the fidelity in the approximation of the initial data and source terms, respectively, as only the relevant left-hand side quantities for our computations have been restated.It is also interesting to ask whether the framework of Ern and Meunier [15] can be extended to yield a posteriori error estimators for higher-order time discretizations of elliptic-parabolic systems (e.g.(3.1)).One might ponder, for instance, the use of the generalized θ scheme ) for which (4.2) is θ = 1 and θ = 1/2 yields the trapezoidal time integration method.Adapting [15] to this context could be approached by generalizing Lemma 2.1 and Theorem 3.1 alongside extending the discrete scheme interpolation, Galerkin residuals, element and jump residuals, Theorem 4.1 and Proposition 4.1-4.3 of [15, Sec.4.1].Though higher-order time discretization schemes are of practical importance, the analytic extension of [15] to this context is a topic for future work.

Numerical convergence and accuracy of error estimators
To examine the accuracy of the computed error estimators and resulting error estimate, we first study an idealized test case with a manufactured smooth solution over uniform meshes.We will consider adaptive algorithms and meshes in the subsequent sections.All numerical experiments were implemented using the FEniCS Project finite element software [3].
We approximate the solutions using Taylor-Hood type elements relative to given families of meshes; i.e. continuous piecewise quadratic vector fields for the displacement and continuous piecewise linears for each pressure.The exact solutions were approximated using continuous piecewise cubic finite element spaces in the numerical computations.
6.1.Convergence and accuracy under uniform refinement.We first consider the convergence of the numerical solutions, their approximation errors and error estimators η 1 , η 2 , η 3 , η 4 under uniform refinement in space and time.We define the meshes by dividing the domain into N × N squares and dividing each subsquare by the diagonal.The errors and convergence rates for the displacement and pressure approximations, measured in natural Bochner norms, are listed in Table 1.We observe that both the spatial and the temporal discretization contributes to the errors, and that all variables converges at at least first order in space and time -as expected with the implicit Euler scheme.For coarse meshes, we observe that the displacement converges at the optimal second order under mesh refinement (Table 1a).
We next consider the convergence and accuracy of the error estimators η 1 , η 2 , η 3 , η 4 for the same set of discretizations (Table 2).We observe that each error estimator converge at at least first order in space-time, with η 2 and η 3 converging at second order in space and η 4 converging at first order in time 1 .
We also define two efficiency indices Ĩeff and I eff with respect to the Bochner and energy norms, respectively, for the evaluation of the approximation error: where Note that we use both Bochner-and energy norms to investigate the practical quality and efficiency of the approximations and estimators as well as in terms of the energy/parameter-weighted norms appearing in the theoretical bound (Proposition 5.2).For this test case, we find Bochner efficiency indices between 1.8 and 5.7, with little variation in this efficiency index between time steps for coarse meshes, and efficiency indices closer to 1 for finer meshes.
The energy-norm efficiency indices I eff are above 1 for all material variations considered.Variations in the specific storage coefficients, Biot-Willis coefficients and transfer coefficients have minimal effect on both the energy-and Bochner norm efficiency indices: the efficiency indices are ∼ 4 for all variations in each of these parameters.For variations in the permeabilities κ i , we observe some reduction in the energy-norm efficiency indices as the permeability is reduced (from 4.1 to 2.6), but that the index value stabilizes around 2.5 for the smaller permeabilities.We observe similar behaviour for the Bochner-norm efficiency index, but with index values of ∼ 0.8 for the smaller permeabilities, and thus efficiency indices below 1.0.For the elastic parameters, the results are quite different.Both the energy-norm and Bochner efficiency indices increase substantially with increasing Lamé parameters µ and λ, though with Bochner efficiency indices increasing more.Remark 6.1.The boundedness of efficiency indices, e.g.I eff and Ĩeff , is canonically provided by the reverse inequality of Proposition 5.2, whereby the the error indicators are bounded above in terms of a constant times the norm of the discretization error.That is, one establishes If the constant of proportionality, in the above inequality, does not involve specific material parameters, then the efficiency indices are robust with respect to variations in those parameters.However, as in the case of the use of higher-order time discretizations, the general Euler-Galerkin elliptic-parabolic framework of Ern and Meunier [15] does not provide this bound and, as a result, its extension to the equations of generalized poroelasticity (MPET), presented herein, is limited in this same regard.The computational experiments of this section suggest that such an estimate will entail a constant of proportionality that scales strongly with both µ and λ.

Adaptive strategy: algorithmic considerations and numerical evaluation
We now turn to consider and evaluate two components of an overall adaptive strategy: (i) temporal adaptivity (only) and (ii) temporal and spatial adaptivity.Our choices for the adaptive strategy can be viewed in light of the observations on the convergence of η 1 , η 2 , η 3 , η 4 for the previous test case, as well as the following characteristics of MPET problems arising in e.g.biological applications: • Mathematical models of living tissue are often associated with a wide range of uncertainty e.g. in terms of modelling assumptions, material parameters, and data fidelity.Simulations are therefore often not constrained by a precise numerical error tolerance, but rather by the limited availability of computational resources.• Living tissue often feature heterogeneous material parameters, but typically with small jumps, and in particular smoother variations than e.g. in the geosciences.The corresponding MPET solutions are often relatively smooth.• Even for problems with a small number of networks J such as single or two-network settings, the linear systems to be solved at each time step are relatively large already for moderately coarse meshes.In light of these points, our target is an adaptive algorithm robustly reducing the error(s) given limited computational resources.We therefore consider an error balancing strategy in which we adaptively refine time steps such that the estimated temporal and spatial contributions to the error is balanced, then refine the spatial mesh to reduce the overall error, and repeat.This approach to spatial adaptivity seeks to balance the computational gains associated with adaptively refined meshes and the computational and implementational overhead costs associated with more sophisticated time-space adaptive methods, see e.g.[1,7].While a full space-time adaptive algorithm could yield time-varying meshes of lower computational cost, the computational costs associated with finite element matrix assembly over separate meshes, and interpolation of discrete fields between different meshes are often substantial.Without time-varying meshes, the blocks of the MPET linear operator can be reused (using the time steps τ n as weights) which may reduce assembly time, and potentially linear system solver times.We note though that there is ample room for more sophisticated time step control methods than we consider here, see e.g.[34] and related works.7.1.Time adaptivity.We consider the time-adaptive scheme listed in Algorithm 1. Overall, for a given mesh T h , we step forward in time, evaluate (an approximation to) the error estimators at the current time step, compare the spatial and temporal contributions to the error estimators, and coarsen (or refine) the time step if the spatial (or temporal) error dominates.
Evaluation of the time-adaptive algorithm on a smooth numerical test case.We evaluate Algorithm 1 using the numerical test case with smooth solutions defined over 3 networks as introduced in Section 6, the default material parameters, and different uniform meshes (defined by 2 × N × N triangles as before).We also verified the adaptive solver by comparing the solutions at each time step and error estimators resulting from rejected coarsening and refinement (resulting in τ n = 0.2 for each n) with the solutions and error estimators computed with a uniform time step (τ = 0.2).We let T = 1.0, and considered an initial time step of τ 0 = 0.2, adaptive weight α η = 0, a coarsening/refinement factor β = 2, and time step bounds τ max = T and τ min = 0.0.
The discrete times t n resulting from the adaptive algorithm, error estimators η n h and η τ n are shown in Figure 2 for different uniform mesh resolutions.For N = 8 (Figure 2a), we observe that the adaptive algorithm estimates the initial time step of 0.2 to be unnecessarily small in light of the dominating spatial error, and coarsens the time step to 0.4 before quickly reaching Algorithm 1 Time-adaptive algorithm 1: Define adaptive parameters α η ∈ [0, 1) and β ≥ 1, τ max > 0 and τ min ≥ 0. 2: Assume that a mesh T h and an initial time step size τ 0 is given.Set t 0 and set the time step iterator n = 0. 3: while t n < T do 4: while True do 5: Set n = n + 1.

6:
Set t * = t n−1 + τ n , and solve (4.2) over T h for (u * h , p * h ) with time step size τ n 7: Compute error estimator approximations at the current time step and set , coarsen the next τ n+1 = βτ n , and break loop. 11: Discard the solution and refine the time step: set 13: else 14: end while 17: end while the end of time (T ).The u − u hτ L ∞ , p − p hτ L ∞ , and p − p hτ L 2 errors are 4.61 × 10 −3 , 3.86 × 10 −2 and 6.83 × 10 −1 .For comparison, with a uniform time step τ = 0.2, the u − u hτ L ∞ , p − p hτ L ∞ , p − p hτ L 2 , and p − π 0 p hτ L 2 (as listed in Table 1) are 4.71 × 10 −3 , 4.38 × 10 −2 , 4.96 × 10 −1 , and 1.33 respectively, and thus the errors with the adaptively defined coarser time step are very comparable -as targeted by our error balancing principle.The picture changes for N = 16 (Figure 2b), in this case the temporal error initially dominates the spatial error, and the time step is reduced substantially initially before a subsequent increase and plateau at 0.1 − 0.2.The value of the adaptive error estimator η 4 is lower than for the uniform solution (1.23 vs 2.17), but the exact errors are comparable between the uniform and adaptive scheme in this case.By setting τ min = τ 0 /4, the unnecessarily high initial time step refinement is limited (Figure 2c), and again comparable errors as for the uniform time step are observed.For higher spatial resolution and thus lower spatial errors (N = 32), similar observations hold (Figure 2d), but now the adaptive solutions approximately halve the exact errors compared to the uniform τ 0 = 0.2 case (as expected).We conclude that the time adaptive scheme efficiently balances the temporal and spatial error, but does little for reducing the overall error -as the spatial error dominates this case.For N = 64 and the same configurations, the adaptive time step reduces to the minimal threshold τ 0 /4 = 0.05 and remains there until end of time T , with the expected quartering of the exact errors compared to the τ 0 = 0.2 case (and the first order accuracy of the temporal discretization scheme).7.2.Spatial adaptivity.For the spatial adaptivity, we use adaptive mesh (h-)refinement based on local error indicators {η K } K∈T h derived from the global error estimators (5.20).In light of the theoretical and empirical observation that η 4 primarily contributes to the temporal error, we will rely on local contributions to η 1 , η 2 and η 3 only for the local error indicators.Specifically, we will let where , η   for each K ∈ T h .The complete space-time adaptive algorithm is given in Algorithm 2. We here choose to use Dörfler marking [13] or a maximal marking strategy in which the γ M of the total number of cells with the largest error indicators are marked for refinement, but other marking strategies could of course also be used.
Remark 7.1.In Algorithm 2, we suggest adapting the mesh in each outer iteration via only (local) mesh refinements.One could equally well consider a combination of local mesh refinement and coarsening, and/or other adaptive mesh techniques such as r-refinement.Indeed, this could be particularly relevant in connection with complex geometries, for which the initial mesh may be overly fine (in terms of approximation power) in geometrically involved local regions.
Evaluation of the space-time-adaptive algorithm on a smooth numerical test case.We evaluate Algorithm 2 using the numerical test case with smooth solutions defined over 3 networks as introduced in Section 6 with the default material parameters.As this is a smooth test case in a regular domain, we expect only moderate efficiency improvements (if any) from adaptive mesh refinement, and therefore primarily evaluate the accuracy of the error estimators on adaptively refined meshes and the balance between temporal and spatial adaptivity.
We set T = 1.0, τ 0 = 0.5, and begin with a 2×4×4 mesh of the unit square as T 0 h .We set = 0, but instead prescribe a resource tolerance L. We first set a fine initial time step τ 0 = T /64, and let β = 2.0, α η = 0.3, τ min = τ 0 /16, and τ max = τ 0 in Algorithm 1.We note that a Dörfler marking fraction γ M of 1.0 yields a series of uniformly refined meshes.For marking fractions between Algorithm 2 Space-time adaptive algorithm 1: Assume that an error tolerance and/or a resource limit L and an initial mesh T = T 0 h are given.Set a marking fraction parameter γ M ∈ (0, 1].2: while True do 3: Set the parameters (τ 0 , α η , β, τ max , τ min ) required by Algorithm 1.

4:
Solve (4.2) over T via the time-adaptive scheme defined by Algorithm 1.

11:
Refine T (locally) based on the markers {y K }. end if 15: end while 0.0 and 1.0, we obtain adaptively refined meshes, yet for this test case, the time step remains uniform throughout the adaptive loop.The resulting errors and error estimates at each adaptive refinement iteration are shown in Figure 3a.We observe that the errors decay as expected, and that the error estimates provide upper bounds for the errors E at each refinement level for all marking fractions tested.
We next let τ 0 = T /4 and γ M = 0.3 (and all other parameters as before), and consider the results of the adaptive algorithm (Figure 3b-3d).We find that the adaptive algorithm keeps the initial time step and refines the mesh only for the first 4 iterations, which substantially reduces the L ∞ H 1 0 displacement approximation error and moderately reduces the L ∞ H 1 0 pressure approximation error.For the next iterations, both the mesh and the time step is refined.The pressure errors seem to plateau before continuing to reduce given sufficient mesh refinement, while the displacement errors steadily decrease.

Adaptive brain modelling and simulation
We turn to consider a physiologically and computationally realistic scenario for simulating the poroelastic response of the human brain.Human brains form highly non-trivial, non-convex domains characterized by narrow gyri and deep sulci, and as such represent a challenge for mesh generation algorithms.Therefore, brain meshes are typically constructed to accurately represent the surface geometry, without particular concern for numerical approximation properties.We therefore ask whether the adaptive algorithm presented here can effectively and without further human intervention improve the numerical approximation of key physiological quantities of interest starting from a moderately coarse initial mesh and initial time step.

Parameter
Value Note E (Young's modulus) Table 3. Material parameters corresponding to a human brain at body temperature.The hydraulic conductances κ are defined in terms of the permeabilities and the fluid viscosities κ j = k j /µ j , µ 2 = µ 1 .Values marked by the * are estimates, yielding physiologically reasonable brain displacements, fluid pressures, and fluid velocities.
Specifically, we let Ω be defined by a subject-specific left brain hemisphere mesh (Figure 4a) generated from MRI-data via FreeSurfer [17] and SVMTK as described e.g. in [26].The domain boundary is partitioned in two main parts: the semi-inner boundary enclosing the left lateral ventricle ∂Ω v and the remaining boundary ∂Ω s (Figure 4b).
Over this domain, we consider the MPET equations (3.3) with J = 3 fluid networks representing an arteriole/capillary network (j = 1), a low-pressure venous network (j = 2), and a perivascular space network (j = 3).We assume that the two first networks are filled with blood, while the third network is filled with cerebrospinal fluid (CSF).8.1.Pulsatility driven by fluid influx.We consider a scenario in which fluid influx is represented by a pulsatile uniform source in the arteriole/capillary network (j = 1): while we set g 2 = g 3 = f = 0. From the arteriole/capillary network, fluid can transfer either into the venous network or into the perivascular network with rates γ 12 , γ 13 > 0, while γ 23 = 0.All material parameters are given in Table 3

Table 4
We assume that the venous network is connected to a low (zero) pressure compartment and set: We assume that the perivascular space is in direct contact with its environment, and set: Last, we model p csf via a simple Windkessel model at the boundary: with compliance C and a resistance R (see Table 3), and where Q is the outflow: Q = ∂Ω u • n ds.After an explicit time discretization, we define at each time step and use (8.7) in (8.2) and (8.5).Finally, we let all fields start at zero.We let T = 2.0 corresponding to two cardiac cycles, and an initial time step of τ 0 = 0.1.The given fluid influx induces pulsatile tissue displacements and pressures in the different networks with varying temporal and spatial patterns (Figure 4c, Figure 5).The brain hemisphere expands and contracts with peak changes in volume dV = Ω div u dx, of up to 1200 mm 3 .The largest displacements occur around the lateral ventricle with peak displacement magnitudes of ≈0.5 mm.The arteriole/capillary pressure varies in space and time with a peak pressure max p 1 of up to 1200 Pa, a pressure pulse amplitude ∆p 1 of ≈ 560 Pa and a pressure difference in space of ≈ 400 Pa.The venous pressure field show similar patterns, though with lower temporal variations and higher spatial variability inducing higher venous blood velocities of above 2.0 mm/s (Figure 5).The perivascular pressure shows a steady increase of up to ≈ 200 Pa at T = 2.0, but only moderate pulsatility and lower fluid velocities than both the arteriole/capillary and venous networks.The local error indicators {η K } K as defined by (7.1) show substantial local error contributions with substantial spatial variation (Figure 6): the values range from the order of 10 3 to 10 9 on the initial mesh T 0 h .This large variation in error indicator magnitude makes the choice of marking strategy important: the Dörfler marking strategy would lead to the marking of perhaps only a handful of cells in this case as the local error indicators for a few cells would easily add up to a significant percentage of the total error.Therefore, we instead choose to employ a maximal marking strategy with a marking fraction γ M = 0.03 for this test scenario.
The adaptive algorithm yields locally refined meshes with around 67 000 cells after one refinement and 198 000 cells after two (Table 4a).The error estimates cf.(5.20) decrease with the adaptive refinement (Table 4b).The contribution from η 2 and η 3 dominates the error estimate, and both of these as well as the total error estimate η seem to halfen for each adaptive refinement level.We also note that in this simulation scenario, for all time steps n and refinement levels, the We also inspect the computed quantities of physiological interest (Table 4a).Using the uniform refinement as an intermediate reference value, we observe that the adaptive algorithm seems to produce more accurate estimates of these quantities of interest even after a single adaptive refinement, and that the quantities of interest after two refinements are more accurate than those of a uniform refinement.The adaptive procedure is therefore able to drive more accurate computation of quantities of interest at a lower or comparable cost as uniform refinement.
Finally, we observe that a single uniform refinement yields a mesh with around 167 000 cells (Table 4a).Thus, even with a small marking fraction of 3%, the mesh growth in each adaptive iteration is substantial.In the Plaza algorithm [28] and other similar conforming mesh refinement algorithms, both cells marked for refinement as well as neighboring cells will be refined to avoid mesh artefacts such as hanging nodes.Therefore, the domain geometry and initial mesh connectivity may strongly influence the adaptive mesh growth, and mesh growth may be more rapid than anticipated, especially in 3D.A more targeted adaptivity and more gradual growth could possibly be achieved with even lower marking fractions, though in the current case the propagation of cell refinement to neighboring cells seems to dominate.In any case, allowing for meshes with hanging nodes could be an effective albeit more disruptive strategy for reducing the computational complexity.
8.2.Pulsatility driven by boundary pressure.We also consider an alternative, more localized, scenario in which, instead of considering a Windkessel model for the CSF pressure and the directly coupled PVS pressure p 3 , we directly prescribe a variation in the boundary PVS pressure.and thus boundary CSF pressure variations of up to ±2 × 133 Pa (corresponding to approx.±2 mmHg) in each cycle.In this scenario, we set the bulk fluid influx to zero (g 1 = g 2 = g 3 = 0.0).We consider otherwise the same experiments as in Section 8.1 and the same adaptive parameters.Also for this case, we observe that the adaptive refinement -even with a small marking fraction and maximal marking yields non-localized marking patterns and a relatively rapid growth in the number of mesh cells.Three adaptive refinements yields meshes with 20 911, 68 608 and 197 975 cells and no refinement of the time steps; numbers which are comparable with the previous case.These results corroborate the observation that the adaptive algorithm drives distributed mesh refinement, and that the spatial errors overall dominate.Moreover, further studies may consider finer initial meshes and further refinements.A prerequisite for this would be robust parallel adaptive refinement algorithms including robust transfer of fields within the mesh hierarchy, and considered the topic of later work.

2. 3 .
Norms and function spaces.Let f , g denote real-valued functions with domain Ω ⊂ R d .If there exists a generic constant C with f ≤ Cg then we write f g.The notation f, g signifies the usual Lebesgue inner product defined by f, g = Ω f g dx, and ||f || = f, f 1/2 is the corresponding norm on the Hilbert space of square-integrable functions

0 ,
and continuous by Lemma 3.2.(3)The embedding of (V a , • a ) into L a follows from Poincare's inequality and the coercivity of a over V a and similarly forV d → L d .(4) c is symmetric by definition, continuous over L d with continuity constant depending on c max , and coercive with coercivity constant depending on c min > 0. (5) The form b given by (3.4b) is clearly bilinear and continuous on V a × L d as

Hypothesis 4 . 1 .
There exists positive real numbers, denoted s a and s d , and subspaces, W a ⊂ V a and W d ⊂ V d equipped with norms || • || Wa and || • || W d , such that the following estimates hold independently of h ∀v

Lemma 4 . 4 .
The discrete two-field variational formulation of the MPET equations (4.2) with the choice of discrete spaces V a,h (4.3) and V d,h (4.4) satisfy the elliptic-parabolic framework hypotheses 4.1-4.3,above.

5. 1 .
Time interpolation.We now recall additional notation for time interpolation (from [15, Sec.4.1]), and rewrite (4.2).Let u hτ denote the continuous and piecewise linear function in time, ) by taking the time-shift of the expressions appearing inside the norm.With the local indicators in hand we immediately have the global indicators and their time-shifted version given by

Figure 2 .
Figure 2. Evaluation of adaptive time stepping for a smooth numerical test case, given uniform meshes and different adaptive parameter configurations.All plots show the approximated error estimators η n h and η n τ at each time step versus adaptive times t n .
Discrete times (and time steps) generated by the space-time adaptive algorithm (same parameters as in Figure3b).

Figure 3 .
Figure 3. Evaluation of the space-time adaptive algorithm on a smooth test case.

( a )
Initial mesh of a brain hemisphere (sagittal view, along positive x-axis) with 20911 cells and 6325 vertices, and a volume of 4.37 × 10 5 mm 3 .(b) Illustration of the semi-inner ventricular boundary (in white), coronal and sagittal clips (view from the yand negative x-axes, respectively).

Figure 4 .
Figure 4.The human brain as a poroelastic medium: meshes, boundaries, and snapshots of solution fields.

Figure 5 .
Figure 5. Left to right, top to bottom: Volume change dV , peak pressure p i and average velocity v i for i = 1, 2, 3 and integrated transfer rates T 12 and T 13 over time for a series of adaptively refined meshes (a1, a2, a3).The opacity indicates the adaptive level: the more opaque, the finer the mesh.

Figure 6 .
Figure 6.Error indicators {η K } K for three levels of adaptive refinement T 0 h , T 1 h , T 2 h the brain simulation scenario.Refinement levels from top to bottom (a1, a2, a3), sagittal views from right and left on the left and right.
are symmetric, coercive, and continuous bilinear forms, thus inducing associated inner-products and norms (denoted by • a and • d ) on their respective spaces.(3) There exist Hilbert spaces L a and L d with V a ⊂ L a and V d ⊂ L d , where the inclusion is dense and such that ||f || La ||f || a for f ∈ V a , and ||g|| L d ||g|| d for g ∈ V d .(4) c : L d × L d → R is symmetric, coercive and continuous; thereby defining an equivalent norm • c , on L d .(5) There exists a continuous bilinear form b 2 e , where the norms || • || K and || • || e represent the usual L 2 , or or d-dimensional L 2 , norm over a simplex, K, and edge, e, respectively.Likewise, the local error indicators associated with the mass conservation equations are 0,T ;Va) .Proof.The above follows from the results of [15, Thm.4.1, Prop.4.1, Thm 4.2] applied in the context of generalized poroelasticity in light of the results of Section 3.

Table 2 .
Error estimators η 1 , η 2 , η 3 , η 4 and their rates of convergence, and Bochner efficiency indices Ĩeff for the smooth 3-network test case under uniform refinement in space (horizontal) and time (horizontal) T = 0.4, τ 0 = T /2.Rate (τ ) is the rate for the finest mesh, under time step refinement.Rate (h) is the rate for the finest time step, under mesh refinement.The diagonal rate (in bold) is the final space-time (diagonal) rate.
and break loop.
j α j p j I) • n = −p csf n on the inner boundary ∂Ω v .(8.2b) for a spatially-constant p csf to be defined below.For the arteriole space, we assume no boundary flux:∇ p 1 • n = 0 on ∂Ω.(8.3)#cells h