Port-Hamiltonian formulations of poroelastic network models

We investigate an energy-based formulation of the two-field poroelasticity model and the related multiple-network model as they appear in the geosciences or medical applications. We propose a port-Hamiltonian formulation of the system equations, which is beneficial for preserving important system properties after discretization or model-order reduction. For this, we include the commonly omitted second-order term and consider the corresponding first-order formulation. The port-Hamiltonian formulation of the quasi-static case is then obtained by taking the formal limit. Further, we interpret the poroelastic equations as an interconnection of a network of submodels with internal energies, adding a control-theoretic understanding of the poroelastic equations.


Introduction
The port-Hamiltonian (pH) framework constitutes an energy-based model paradigm that offers a systematic approach for the interactions of (physical) systems with each other and with the environment and thus extends Hamiltonian systems to open physical systems. The pH structure, introduced originally in [32], see [37,44] for an overview, provides a geometric description of the model in terms of a Dirac structure, which directly encodes system properties such as passivity and stability into the system of equations. Structure-preserving methods, such as space-discretization [18,38], model-order reduction [6,15,23,35,47], and time-integration [27,34], ensure the preservation of these properties throughout numerical approximation schemes. Recently, the definition of pH systems was extended to cover implicit systems [7,34,43,45], resulting in so-called port-Hamiltonian differential-algebraic equations (pH-DAEs). Although this system class covers a wide range of equations and simplifies the mathematical theory in many aspects, almost all results are obtained for finite-dimensional DAEs. Extensions of the pH framework to (constrained) infinite-dimensional systems typically do not exist in a general form but rather consider particular applications or special model classes [5,18,24,26,30,36].
In view of extending the range of application classes, this paper is devoted to the generation and analysis of the pH structure for partial differential equations which model the deformation of porous media saturated by an incompressible viscous fluid, namely poroelasticity [8,17,20,39]. These equations represent coupled systems of different types of differential equations and, hence, form a so-called partial differential-algebraic equation (PDAE) in total [3]. The main interest in terms of applications can be found in the field of geomechanics [49], often in connection with heterogeneous structures [2,11,22]. Moreover, also the displacement of a material due to temperature changes can be modeled with a mathematically equivalent model called thermoelasticity, cf. [9,31].
The full strength of the pH approach, however, becomes visible in the network case, including multiple pressure variables. Such models typically appear in medical applications such as cerebral infusion tests [21,41] or the investigation of cerebral edema [46]. In both cases, the brain is modeled as a poroelastic medium saturated by different fluids. The formulation as a pH system then allows an interpretation as a network of submodels with internal energies and interconnections.
This paper is structured as follows. After introducing the classical two-field formulation of poroelasticity in terms of a weak operator formulation in Section 2, we discuss its extension to the network case. Although the system equations are quasi-static, we include the usually omitted second-order term to find a pH formulation. In Section 3 we recall the finite-dimensional pH-DAE framework and show that the poroelastic equations can be cast in a corresponding infinite-dimensional structure. In Section 4 we then derive the pH structure from an interconnection point of view. This means that we already consider the original poroelastic system as a coupled system of pH subsystems. Finally, we summarize the paper in Section 5.

Poroelastic Network Models
In this section, we introduce the model equations of poroelasticity and a corresponding operator formulation, which corresponds to the weak form. In view of the subsequent reformulations, it turns out that the full model is better suited than the often considered quasi-static form. This model is then extended by additional pressure variables leading to the multiple-network case.
2.1. Two-field formulation. The original model of linear poroelasticity in a bounded Lipschitz domain Ω ⊆ R d with d ∈ {2, 3} (and boundary ∂Ω) was introduced in [8], see [39] for a modern formulation covering different modeling components. Considering a time horizon 0 < T < ∞ and setting T := [0, T ], one wants to determine the displacement field u : T × Ω → R and the pressure p : T × Ω → R satisfying the system of coupled partial differential equations with Dirichlet boundary conditions and initial conditions p( · , 0) = p 0 , u( · , 0) = u 0 , as well as ∂ t u( · , 0) =u 0 . Within this system, the stress tensor σ models the linear elastic stress-strain constitutive relation σ(u) = 2µ ε(u) + λ (∇ · u) I, ε(u) = 1 2 ∇u + (∇u) T with the Lamé coefficients µ and λ and the identity tensor I. Further, α denotes the Biot-Willis fluid-solid coupling coefficient, M the Biot modulus, κ the permeability, ρ the density, and ν the fluid viscosity. The right-hand sideĝ represents an injection or production process andf denotes the volume-distributed external forces.
For the weak formulation of the poroelastic model, we homogenize the boundary conditions and assume that u b and p b are zero. Inhomogeneous boundary conditions will be shortly discussed in Remark 3.6. Let us introduce the spaces where L 2 (Ω) is the space of square-Lebesgue integrable functions on Ω and H 1 0 (Ω) is the Hilbert space of functions that have a weak derivative in L 2 (Ω) and satisfy zero boundary conditions on ∂Ω. Hence, V and Q already include the respective homogeneous boundary conditions. With the corresponding dual spaces V * and Q * we have two Gelfand triples at hand, namely V, H V , V * and Q, H Q , Q * , cf. [48,Ch. 23.4]. For both triples we will use the notion · , · for the respective duality pairing. The resulting embedding is denoted by I : Q → Q * , I · , · = ( · , · ) HQ .
We introduce the operators Y : as well as the differential operators A : V → V * , K : Q → Q * , and D : where we use for A the classical double dot notation from continuum mechanics. Note that, due to the integration by parts formula, the operator D can also be considered in the form D : H V → Q * . For the right-hand sides we write f (t) := Ωf (t) · dx ∈ H * V and g(t) := Ωĝ (t) · dx ∈ H * Q . It follows from Korn's inequality [16,Th. 6.3.4] that the operator A is elliptic and, thus, invertible. Further, it is self-adjoint and bounded in V. Similarly, the operator K is selfadjoint, elliptic, and bounded in Q. The operator M mainly contains the multiplication by the constant 1/M > 0 and is thus self-adjoint, elliptic, and bounded in the pivot space H Q , although it may also be interpreted as an operator M : Q → Q * . In the latter setting, however, the operator is only continuous and not elliptic. Similarly, the operator Y represents the multiplication with the constant ρ > 0 and thus is self-adjoint, elliptic, and bounded in H V . Finally, the coupling operator D equals, up to the multiplicative prefactor α, the divergence and is thus bounded.
The weak form of (2.1) can be stated as follows: Determine u : T → V and p : T → Q such that for almost all t ∈ (0, T ) we have in Q * (2.2b) for initial conditions p 0 ∈ H Q , u 0 ∈ V, andu 0 ∈ H V . Note that D * : H Q → V * denotes the dual operator of D, i.e., D * q, v = q, Dv for all v ∈ V and q ∈ Q. Assuming right-hand sides f ∈ L 2 (T; H V ), g ∈ L 2 (T; H Q ), suitable solution spaces (for the weak formulation) read u ∈ L 2 (T; V), p ∈ L 2 (T; Q) withu ∈ L 2 (T; H V ),ü ∈ L 2 (T; V * ), andṗ ∈ L 2 (T; Q * ). Sobolev embeddings and the existence theory of [29,Ch. 3,Sect. 8.4] then imply that p ∈ C(T; H Q ) and u ∈ C(T; V),u ∈ C(T; H V ) such that initial conditions are meaningful.
Within this paper, we will study system (2.2) in the corresponding formulation in terms of operator matrices, i.e., Remark 2.1. In several applications, see for instance [14], the permeability κ, or more precisely, the hydraulic conductivity, depends on the dilatation ∇ · u, i.e., κ = κ(∇ · u). In this case, the term Kp becomes nonlinear. Assuming that there exist positive constants κ − , κ + satisfying 0 < κ − ≤ κ(ξ) ≤ κ + for all ξ ∈ R still renders the operator K(u) selfadjoint and elliptic, which is sufficient for the forthcoming analysis. For instance, this is the case for the Kozney-Carmen type hydraulic conductivity [13,14], see also [4] for a corresponding numerical setup.

2.2.
Quasi-static formulation. In many applications, the coefficient ρ in (2.1a) is considered small enough to be neglected. The formal limit ρ → 0 then leads to the so-called quasi-static formulation, studied, for instance, in [18,39]. In view of the operator formulation (2.2), the quasi-static model is given by Remark 2.2. For the quasi-static formulation (2.4) there is no need to prescribe an initial value foru. In fact, it it not even necessary to prescribe an initial value for u, since u 0 is uniquely determined from (2.4a), see for instance [3] and the forthcoming Section 3.5. The dependence of u 0 on p 0 is a typical consistency condition known for PDAEs.
If f is differentiable in time, then, since the operator A is invertible, we can formally solve (2.4a) for u and insert the result into (2.4b). Then one obtains a parabolic equation for the pressure p given by Note that the differentiability of f may be a critical factor in this approach.
We would like to emphasize that system (2.4) covers general elliptic-parabolic problems [3]. Clearly, this includes (quasi-static) linear poroelasticity and thermoelasticity. In contrast, the full model (2.2) covers a general hyperbolic-parabolic coupling.
2.3. Multi-field network systems. For medical applications, the poroelastic model (2.1) is often extended by additional pressure variables, e.g., to distinguish vessel types in the investigation of cerebral edema, cf. [46]. Again such applications usually come with inhomogeneous boundary conditions, and a homogenization is necessary. For our analysis, however, we restrict ourselves to homogeneous Dirichlet boundary conditions to focus on the structure of the equations.
In the multi-field network case we want to compute the displacement u : T × Ω → R and pressures p i : with initial conditions for u,u, and all p i . Note that the second equation has to be considered for i = 1, . . . , m, leading to m + 1 equations in total. The parameters β ij model exchange rates from one subsystem to the other and are usually small. In several applications, symmetric exchange rates β ij = β ji are considered, see e.g. [42].
To shorten notation, we introduce the matrix B : Thus, B contains all exchange rates of the system in the off-diagonal elements and the row sums are all equal to zero by construction. In this way, the exchange conditions resemble Kirchhoff's laws in electric or energy transport networks [19,43,44].
For the weak formulation of (2.6) we define operators K i : Q → Q * and D i : V → H * Q as in Section 2.1 with the same properties as K and D, respectively. We further introduce the block operators Note that, in order to allow a compact notion of the network case, we write Q m and H m Q for the product spaces Q × · · · × Q and H Q × · · · × H Q , respectively. Further, B ⊗ I mimics the Kronecker product, i.e., for Proof. The claim on the symmetry of K B is obvious, since all operators K i are self-adjoint. Let c denote the (uniform) ellipticity constant for all K i and C β : Here C Q֒→HQ denotes the constant of the continuous embedding Q ֒→ H Q .
Assumption 2.4. Throughout this paper, we assume that the operator K B is elliptic and thus invertible.
Note that Assumption 2.4 does not include the symmetry of B. For the nonsymmetric case we make the following observation.
The proof follows immediately from the identity K B q, q = K Bsym q, q + 1 2 (B ⊗ I − B T ⊗ I)q, q and the observation that Let us emphasize that Lemma 2.3 indicates that small exchange rates in the sense that with c denoting the uniform ellipticity constant of all K i , is a sufficient condition for the ellipticity in Assumption 2.4.
The weak formulation of (2.6) in terms of the introduced operator matrices has the form: Determine u : T → V and p : T → Q m such that for almost all t ∈ T the operator equation is satisfied. Here we have used the notion g := [g 1 , . . . , g m ] T , where the g i are defined in terms ofĝ i as in Section 2.1.
Thus, the network case has the same structure as the two-field poroelastic model, cf. system (2.3). A direct consequence is that a semi-discretization in space of the multi-field network model (2.7) yields an operator equation of the same structure as the two-field model.

2.4.
Quasi-static network case. Similar to the two-field case considered above, we briefly discuss the case where the coefficient ρ is assumed to be small and thus neglected. This then yields the quasi-static network model As in the two-field model, a consistency condition then uniquely determines u 0 .
Assuming once more that f is differentiable in time, we can use the invertibility of A to eliminate the first equation. This then leads to a system of m coupled parabolic equations, namely Mṗ To restore the displacement u, one computes the solution to the elliptic problem Au = D * p + f . Example 2.6 (Poroelastic brain model). The example in [46] considers ρ = 0, m = 4 pressure variables, and no external forces or injections, i.e., f ≡ 0 and g i ≡ 0. This means that the system is driven only by its hydrostatic pressure gradients.

Port-Hamiltonian descriptor system formulation
To render the paper self-contained, we recall the definition and some results for pH-DAEs, also called port-Hamiltonian descriptor systems. For the finite-dimensional setting, this is presented in [34], see also earlier versions in [7,43]. Afterwards, we show that the poroelastic equations can be formulated in a corresponding way, leading to a port-Hamiltonian partial differential-algebraic equation (pH-PDAE).
3.1. Port-Hamiltonian framework. We start with the definition of a pH-DAE in the finite-dimensional framework, which we will mimic in the pH-PDAE setting. We present a special case of the definition given in [7] that is sufficient for our analysis, using a slightly different notation than therein. For a generalization to time-dependent and nonlinear descriptor systems we refer to [34].
called the dissipation matrix, is positive semi-definite.
Note that (3.1) can be equivalently written as We thus obtain the following equivalent formulation of a pH-DAE. Let us emphasize that if no feedthrough term is present, i.e., D = 0 in (3.1) or S = N = 0 in (3.2), then the dissipation matrix (3.3) being semi-definite immediately implies P = 0 such that the pH formulation is given as with skew-symmetric J and symmetric positive semi-definite E and R. For our remaining analysis, the latter formulation is sufficient, serving as the main reason to restrict ourselves to systems of the form (3.4).
The class of pH-DAEs, respectively the extensions presented in [34], have many nice properties, see also [7]. The general definition in [34] extends to weak solutions and infinite dimension, the structure is invariant under state-time diffeomorphisms, the structure is invariant when making the system autonomous, and the class can be described via a Dirac structure. Furthermore, pH-DAEs satisfy a power balance equation and dissipation inequality. [7,34]). Consider a pH-DAE of the form (3.1). Then the power balance equation

Theorem 3.3 (Dissipation inequality for pH-DAEs
holds along any solution z, for any input v. In particular, the dissipation inequality holds.
Another important property, see [34], is that the class of pH-DAEs is invariant under power-conserving interconnection. To see this, consider two pH-DAEs Hamiltonians H i : R n i → R, for i = 1, 2, and assume that the two systems are interconnected via Then, defining the aggregated state, input, and output, respectively, via and aggregated system matrices J := diag(J 1 , J 2 ), R := diag(R 1 , R 2 ), G := diag(G 1 , G 2 ), the coupled system is given by Writing F = F skew + F sym with skew-symmetric part F skew and symmetric part F sym , we immediately observe that (3.7) is again a pH system if, and only if, R − GF sym G T is positive semi-definite. A sufficient condition to retain the pH structure is thus to require that F sym is negative semi-definite.
Finally, a key property for our later analysis is that the pH-DAE structure is invariant under Galerkin projection. This will form the basis for our forthcoming analysis, where we show that the poroelastic (network) equations introduced in Section 2 exhibit such a pH structure. Since we analyze the equations in the operator formulation, which we will then call a pH-PDAE, we cannot follow the exact wording but only the spirit of Definition 3.1, respectively Lemma 3.2. Nevertheless, applying a spatial discretization by a finite element method afterwards will retain the structural properties such that the semidiscretized models are then pH-DAEs. Since the operator matrices are constant in time and there is no feed-through term, we consider the special case that the matrices P , S, and N are zero. As a result, the requirements of Definition 3.1, respectively Lemma 3.2, reduce to J being skew-symmetric and E, R being symmetric as well as positive semidefinite.
3.2. PH formulation of the two-field model. In this section, we discuss a reformulation of the operator form of the poroelastic equations, which exhibits the pH structure in the spirit of Definition 3.1. Consequently, a Galerkin projection of the operator form will then lead to a (finite-dimensional) pH-DAE.
Starting with the second-order formulation (2.3), we first perform a first-order (mixed) formulation by introducing w =u as new variable and adding the equation Aw = Au to the system. Recall that in our case A * = A is elliptic and thus, invertible. Assuming ρ > 0, we obtain the (implicit) first-order system Note that the three operators Y, A, and M appearing on the left-hand side are invertible. A natural Hamiltonian (energy function) associated with this system is given by (3.9) H(w, u, p) := 1 2 Yw, w + Au, u + Mp, p , where 1 2 Yw, w describes the kinetic part of the energy and 1 2 Au, u + 1 2 Mp, p the potential energy. Let us emphasize that although, formally, the two-field model does not feature any control inputs, we will view the source terms f and g as external inputs (respectively forces, cf. [39]) to the system and complement the system with the powerconjugated output, i.e., with We immediately observe that the first-order formulation (3.8) together with the powerconjugated output equation (3.10) emulates the pH-DAE structure with P = 0, S = N = 0, Recall that the operators A, M, K, and Y are positive definite. We expect, similarly to the finite-dimensional case considered in Theorem 3.3, a power balance equation and a dissipation inequality. Let us demonstrate this as an example for the formulation (3.8). In the other formulations that we discuss in the following, the result and proof is analogous.
where the last equality follows from v, y = f, w + g, p .
Note that the pH structure also allows the operators M and Y to become singular, but still semi-definite. In particular, we immediately have the structure also in the (formal) limiting situation ρ = 0. Furthermore, the pH-PDAE structure still holds if K is not selfadjoint, it suffices that the symmetric part of K is positive semi-definite. In addition, we immediately observe that the power balance equation (3.11) is still satisfied for a nonlinear operator K as discussed in Remark 2.1.
The corresponding Hamiltonian reads and is thus equal to the original energy defined in (3.9). Further note that also the output coincides, i.e., y = y. In contrast to the formulation (3.8), however, the pH-PDAE (3.12) only requires the usual regularity assumptions, namelyu ∈ L 2 (T; H V ).
Remark 3.6. If it is not possible or not desirable to homogenize the Dirichlet boundary conditions u b and p b in (2.1c) and (2.1d), respectively, then we can still obtain a pH formulation that explicitly encodes the boundary conditions. In more detail, following [1,40], we can use the trace-operator to write the boundary condition as an operator equation, which is then coupled to the original system via a Lagrange multiplier. To retain the pH structure, we have to choose a negative sign for the Lagrange multiplier and add the derivative of the boundary operator equation for the displacement u. With this choice, the additional terms for the boundary condition enter the skew-symmetric operator matrix J .
3.3. PH formulation of the quasi-static case. First, we observe that the quasi-static case, i.e., ρ = 0, is included in the pH formulation (3.8), by setting Y = 0. Nevertheless, the first-order formulation (3.8) results from introducing the derivative of u as an additional variable. It is worth noting that this strategy also works for the corresponding three-field formulation of poroelasticity, cf. [8]. This formulation considers a mixed (or partitioned) version of the operator K with the fluid flux (or Darcy velocity) as an additional variable. If we directly start with the quasi-static case (2.4) and do not introduce a new variable for the derivative of u, then we can obtain another pH-PDAE formulation. We have already seen that one may eliminate the variable u, which leads to the parabolic equation (2.5). Also this has a pH structure but does not cover the full information due to the missing displacement variable. To reveal the structure of the entire system, we consider a reformulation. Since K is invertible by assumption, we may introduce a new variable q via Assuming sufficient regularity of the given data, it is easy to see that the pairs (u, p) and (u, q) are equivalent in the sense that one can convert them into each other by With this equivalence, we can write the quasi-static formulation (2.4) equivalently in terms of the new variable, as The operator-matrix on the right-hand side is self-adjoint and negative semi-definite revealing that (3.14) is already in pH form. Note, however, that the operator KM −1 K occurs on the right-hand side, which requires strong spatial regularity and is thus not favorable. Instead, we observe that the system matrix in (3.14) is a Schur complement, i.e., we have the identity Hence, if we combine (3.13) and (3.14) and rearrange the equations we obtain the extended system We emphasize that this gives a pH-PDAE in the spirit of Definition 3.1 with Hamiltonian This can be verified immediately by noting that the operator matrices on the left and the second on the right are self-adjoint and positive semi-definite. Further, the first operator matrix on the right is skew-adjoint. Let us emphasize that the formulation (3.15) is not a pH-PDAE in the sense of our definition if K is not self-adjoint. 3.4. PH formulation of the network case. In Section 2.3 we have seen that the multifield network case can be formulated in the same way as the two-field model if the operators are adjusted accordingly. More precisely, we simply need to replace the operators K, M, and D by the operator matrices K B , M, and D, respectively. By Assumption 2.4 we know that K B is elliptic. As a result, the pH formulation of (2.7) is given by As before, we may reinterpret the right-hand sides as external inputs, leading to the power-conjugated output consisting of w and p. The corresponding Hamiltonian has the form H(w, u, p) := 1 2 Yw, w + Au, u + Mp, p Finally, the pH formulation of the quasi-static network case immediately follows by setting Y = 0. Alternatively, we can use the pH formulation for the quasi-static case derived in Section 3.3. Since the operator K B appears on the left-hand side, however, this formulation is only valid for a symmetric network coupling matrix B.

DAE structure and index.
A spatial discretization of the operator equations (3.8), (3.12), or (3.15), e.g., by the finite element method, yields a (finite-dimensional) pH-DAE. However, it is interesting to note that the DAE structure and, in particular, the differentiation index vary in the different formulations. Recall, see e.g. [10,28], that a constant coefficient DAE Eż = Az + k with a regular pair (E, A) (i.e. det(λE − A) is not identically zero) has differentiation index ν = 0 if E is invertible. It has differentiation index one if W T AV is invertible, where V is a matrix that spans the kernel of E and W is a matrix that spans the kernel of E T . Otherwise, it has differentiation index ν ≥ 2.
The spatial discretization of (3.8) with discretization parameter h yields the pH-DAE Here, K A and K K denote the stiffness matrices corresponding to the operators A and K, respectively, and M Y , M M are the mass matrices corresponding to Y and M. Since the finite element method is a Galerkin projection, we immediately observe that these matrices are symmetric. Under reasonable assumptions on the spatial discretization, the matrices are even positive definite (in the case ρ > 0). If ρ > 0, then (3.18) is an implicit equation with nonsingular matrix on the left-hand side. Thus, it is a pH-DAE of differentiation-index ν = 0. If ρ = 0, then the operator Y and the corresponding mass matrix M Y are zero, i.e., the space-discretized system (3.18) reads where we have interpreted the external forcing as control variables and, as before, added the power-conjugated output equation. As before, the mass matrices M u and M p are assumed to be symmetric positive definite. Since the upper-left block in the matrix on the right-hand side of (3.19a) is zero, the pH-DAE cannot be of index ν = 1. Multiplication of the second and third equation with K −1 A and M −1 M , rearranging the equations and variables details that (3.19a) is in Hessenberg form, see [10]. We immediately conclude that (3.19a) has differentiation index ν = 2. The constraint equation K A u h = D T p h + f h is thus complemented with the hidden constraint As a direct consequence, initial values w 0 h and u 0 h for w h and u h , respectively, have to satisfy the consistency conditions showing that it is sufficient to prescribe an initial value p 0 h for p h . We emphasize that this is not specific to the discretized pH-DAE, but also applies to the pH-PDAE (3.8). In this case, the initial values w 0 and u 0 are given implicitly by Remark 3.8. Although the semi-discrete pH-DAE (3.19) has differentiation index ν = 2, it can be regularized via output feedback, and thus has differentiation index ν = 1 in the behavior framework [12]. In more detail, consider a feedback The resulting closed-loop system is given by showing that (3.22) has differentiation index ν = 1, whenever F 11 is nonsingular, and has pH structure, whenever the symmetric part of F 11 is negative semi-definite.
If one does not introduce the new variable w, then it has been shown in [3] that the spatial discretization of the original quasi-static formulation (2.4) results in a system of differentiation index ν = 1. We emphasize that the original quasi-static formulation (2.4) only encodes the consistency condition (3.20) for u but not the consistency condition (3.21) for w (respectivelyu). Remark 3.9. If the inhomogeneity f is sufficiently smooth in time, then there is also another pH-PDAE like formulation. Going back to the original second-order model (2.3) and differentiating the first equation with respect to time yields the third-order (in time) equation This is a formulation where all coefficients are positive semi-definite, except for one which has a positive semi-definite symmetric part. This is the structure of the higher-order dissipative Hamiltonian systems discussed in [33]. In principle, we could work directly with this formulation and do not transform to first order, but we will not discuss this formulation further here.

Interconnection of subsystems
It has been shown in [34], see also Section 3, that the (energy-preserving) interconnection of pH-DAEs is again a pH-DAE. In the following, we show how the presented pH-PDAE formulations can be obtained via the interconnection of subsystems via a suitable output feedback relation.

4.1.
Interconnection in the two-field case. In the two-field case, we construct an interconnection of a hyperbolic (or elliptic) and a parabolic equation. We treat the two cases ρ > 0 and ρ = 0 together. First, consider the system which, by going to the first-order formulation and adding an output equation, can be written as For the system in p we consider the (parabolic) system The Hamiltonian (3.9) of the coupled system is the sum of the Hamiltonians in (4.4).

4.2.
Interconnection for the alternative quasi-static case. If we are solely interested in the quasi-static case, i.e., ρ = 0 and Y = 0, then (4.1) reduces to the elliptic equation In this case, we may choose the output as y u = u. In terms of Definition 3.1, respectively Lemma 3.2, we have R = A, B = Id, and zero operators otherwise. The associated Hamiltonian in this case is constant in time and we may for simplicity choose H u (u) = 0.
To obtain the alternative pH formulation presented in Section 3.3, we consider for p the last two equations in (3.15), given by with associated Hamiltonian H p (p, q) = 1 2 Kq, q . Choosing a suitable output feedback for v u , v p , v p , we recover (3.15) as stated in the following immediate result.
The Hamiltonian (3.16) of the coupled system is the sum of the Hamiltonians of the subsystems.

4.3.
Interconnection in the network case. Clearly, we obtain all the different pH-PDAE formulations similarly as in the two-field case, since the multi-field network structure is precisely the same as that of the two-field system, provided that Assumption 2.4 holds, i.e., we have small exchange rates. To perform the interconnection via output feedback, replace the parabolic equation  If the exchange rates are symmetric, i.e., B = B sym , then we can also use the alternative quasi-static pH formulation derived in Section 3.3. Since the exchange rates, encoded in the matrix B, then appear on the left-hand side in the pH formulation (3.15), the pH formulation cannot be obtained directly via feedback interconnection. Instead, we can first couple the pressure equations (4.7), then introduce the new variable q as in Section 3.3, and afterwards obtain the final form with the coupling from Theorem 4.2. An illustration of the two pH formulations for the quasi-static network case is presented in Figure 4.1.

Summary
We have studied Biot's poroelasticity model from an energy-based perspective and introduced different pH formulations. Our operator matrices mimic the finite-dimensional definition of pH-DAEs such that a spatial discretization via finite elements yields a pH-DAE. To obtain the required properties, we have used the commonly omitted second-order term and performed a first-order reformulation. The quasi-static pH formulation is then obtained by formally setting the second-order term zero, which does not affect the pH structure. The formulation as pH-PDAE is favorable because it automatically induces a dissipation inequality and is invariant under Galerkin projection, making our formulation amendable to structure-preserving discretization and model-order reduction. Besides, the pH formulation offers a system-theoretic interpretation of poroelastic network models. In particular, we have shown that the pH formulation of the network models can be obtained via output feedback interconnection of the different submodels with internal energies. We emphasize that such an interconnection is not directly possible with the quasi-static model typically studied in the literature, since the coupling of the subsystems requires the derivative of the displacement.