GLOBAL EXISTENCE OF WEAK SOLUTIONS TO UNSATURATED POROELASTICITY

We study unsaturated poroelasticity, i.e., coupled hydro-mechanical processes in variably saturated porous media, here modeled by a non-linear extension of Biot’s well-known quasi-static consolidation model. The coupled elliptic-parabolic system of partial differential equations is a simplified version of the general model for multi-phase flow in deformable porous media, obtained under similar assumptions as usually considered for Richards’ equation. In this work, existence of weak solutions is established in several steps involving a numerical approximation of the problem using a physically-motivated regularization and a finite element/finite volume discretization. Eventually, solvability of the original problem is proved by a combination of the Rothe and Galerkin methods, and further compactness arguments. This approach in particular provides the convergence of the numerical discretization to a regularized model for unsaturated poroelasticity. The final existence result holds under non-degeneracy conditions and natural continuity properties for the constitutive relations. The assumptions are demonstrated to be reasonable in view of geotechnical applications. 1991 Mathematics Subject Classification. 35K61, 65M12, 74F10, 76S05. ...


Introduction
Strongly coupled hydro-mechanical processes in porous media are occurring in various applications of societal relevance within, e.g., geotechnical, structural, and biomechanical engineering. Examples for instance are soil subsidence due to groundwater withdrawal, induced seismicity in geothermal reservoirs, swelling and drying shrinkage of concrete, and deformation of soft biological tissue components, to mention a few.
In the field of porous media, such microscopically complex processes are typically modeled by a continuum mechanics approach [6]. The multi-phasic solid-fluid mixture is considered a homogenized continuum, and both geometry, skeleton, and fluid properties are averaged over representative elementary volumes, consisting of a mixture of solid and fluid particles. Ultimately, the interaction of the different microscopic constituents is described by macroscopic effective equations. The simplest macroscopic model accounting for the two-way coupling of single-phase flow and elastic deformation in a porous medium is Biot's linear quasi-static consolidation model. Its phenomenological derivation Keywords and phrases: Poroelasticity; Biot model; variably saturated porous media; Richards' equation; existence of solutions; convergence of mesh-based discretization it is a common discretization in the field of poroelasticity, both in the research literature [14,45] as well as in the industry; and important for this work, (ii) the specific choice of the discretization will allow for straightforward cancelling of the non-linear coupling terms, which will be crucial for deriving a priori estimates.
The analysis here relies on compactness arguments. Due to the nonlinear coupling terms, one needs a sufficient regularity of the solution, and appropriate a priori estimates. To this prior to the discretization, a physicallymotivated regularization is introduced, accounting for the primary and secondary consolidation of variably saturated porous media with compressible grains [17]. From a mathematical point of view, the regularization transforms both balance equations into non-degenerate parabolic ones, and thereby naturally leads to higher regularity in time. Similar ideas have been used previously, e.g., for saturated poroelasticity [5,18,24,31,38], and for unsaturated flow in porous media [1,32,34,36].
In the remaining steps of the proof, a priori estimates are derived and compactness arguments are utilized in order to deduce, first, the convergence of the finite element/finite volume discretization to a (continuous) regularized solution for vanishing discretization parameters, and furthermore, convergence of the regularized solution to a weak solution of the transformed model for vanishing regularization parameters. Thereby, the existence of a weak solution to the transformed model will be proved. In the discussion of the limit of the spatial discretization specific finite volume techniques are utilized, inspired by [21,23,37].
The present analysis requires an overall parabolic character of the coupled problem and natural continuity properties for the non-linearities. When passing the regularization to zero, those are ensured under specific material assumptions. The most important one requires the presence of a compressible fluid or solid grains, and a sufficiently stiff bulk. In the appendix, the assumptions are demonstrated to be satisfied for non-degenerating constitutive relations typically utilized in practical applications. Furthermore, focusing on the non-linear coupled character, some simplifying assumptions are made as the presence of an isotropic material, no gravity, and merely homogeneous essential along with inhomogeneous natural boundary conditions. The rest of the paper is organized as follows. In Sec. 2, the model is introduced as derived in the engineering literature, and then transformed using the Kirchhoff transformation. In Sec. 3, the notion of a weak solution to the transformed problem is introduced, and the main result is stated: existence of a weak solution to the transformed problem under certain model assumptions and non-degeneracy conditions. The idea of the proof, consisting of six steps, is presented. The details of those six steps are the subject of the remaining Sec. 4-Sec. 9. In Sec. 10, a brief numerical test is provided illustrating the convergence of the used numerical scheme. In Appendix A, the feasibility of the required assumptions for the main result are discussed for widely used constitutive models from the literature. In addition, in Appendix B, technical results from the literature used in the proof of the main result are recalled for a comprehensive presentation.

Mathematical model for unsaturated poroelasticity
We consider a continuum mechanics model for unsaturated poroelasticity, a particular simplification of general multi-phase poroelasticity. As Biot's consolidation model, it is based on the fundamental principles of momentum and mass balance combined with constitutive relations. The model accounts for slight compressibility of the fluid as well as changing porosity due to compressible grains and volumetric deformation, as well pore pressure acting on the solid skeleton in the context of unsaturated media. It is valid under the assumptions of infinitesimal strains and the presence of two fluid phases, of which one is an active and the other is a passive phase; the displacement of the passive phase does not impede the advance of the active phase and can be neglected. Effectively, the model non-linearly couples the Richards equation and linear elasticity equations. For a detailed derivation, we refer to the textbooks [17,26].
In the following, we recall the mathematical model employing the mechanical displacement and fluid pressure as primary variables. Additionally, the problem is transformed using the Kirchhoff transformation, a standard tool for the analysis of non-linear diffusion problems, cf., e.g., [1]. The resulting model will be subject of the subsequent analysis. The momentum and mass balance equations as derived in [26] read on Q T : φ∂ t s w (p w ) + φc w s w (p w )∂ t p w + 1 N s w (p w )∂ t p pore (p w ) + αs w (p w )∂ t ∇ · u + ∇ · q = h, (2.1b) where u and p w are the primal variables, and denote the structural displacement and the fluid pressure (of the active phase), respectively. Furthermore, the volumetric flux q is described by the generalized Darcy law q = −κ abs κ rel (s w (p w )) (∇p w − ρ w g) .

(2.2)
Moreover, in (2.1) and (2.2), µ and λ denote Lamé parameters, ε(u) is the linearized strain tensor, α ∈ [0, 1] is the Biot coefficient, p pore is an arbitrary pore pressure (in [26], the derivation is merely formulated using the socalled averaged pore pressure), f is an external volume force density; φ is the porosity, s w is the fluid saturation, c w ∈ [0, ∞) is the constant storage coefficient associated to fluid compressibility, N ∈ (0, ∞] is the constant Biot modulus associated to the compressibility of solid grains, q is the volumetric flux, and h is an external source term; κ abs is a spatially varying absolute permeability, κ rel is the relative permeability, ρ w is the fluid density, and g is the gravitational acceleration. Typical ranges of values in geotechnical applications are µ, λ ∼ 10 7 ..10 11 Pa, φ ∼ 0.01..0.7, c w ∼ 10 −11 ..10 −7 Pa −1 , N ∼ 10 10 ..10 12 Pa, α ∼ 0.05..1, κ abs ∼ 10 −21 ..10 −9 m 2 ; see [17] for concrete examples. To briefly explain the structure of (2.1), we note that the total stress in (2.1a) is defined via the Biot effective stress [17]; the first four terms in (2.1b) correspond to the weighted change in fluid mass ρ −1 w ∂ t (φs w ρ w ). They are obtained after applying the product rule, employing ρ −1 w ∂ t ρ w = c w ∂ t p w for slightly compressible fluids, as well as the relation ∂ t φ = α∂ t ∇ · u + 1 N ∂ t p pore for varying porosity [17,26]. Constitutive laws are assumed to be given for the fluid saturation s w , the relative permeability κ rel , and the pore pressure p pore . The medium is assumed homogeneous, aside from a heterogeneous matrix permeability. Consequently, the constitutive relationships remain the same everywhere. Commonly accepted in practice are the van Genuchten-Mualem relations [25,30] with a non-decreasing saturation and a possibly degenerate relative permeability κ rel (s w ) := κ vG (s w,res + s ε ), s w < s w,res + s ε , κ vG (s w ), s w ≥ s w,res + s ε , (2.3d) where m vG ∈ (0, 1), n vG = (1 − m vG ) −1 , and α vG > 0 are model parameters, s w,res ∈ [0, 1) is the residual saturation, s w,eff is the effective saturation, and s ε ∈ [0, 1 − s w,res ) is a regularization parameter, ensuring non-degeneracy when chosen positive. In addition, we consider the widely-used equivalent pore pressure [17] p pore (p w ) := pw 0 s w (p) dp. (2.4) In the fully saturated regime, when p w ≥ 0, the equivalent pore pressure reduces to the fluid pressure, which is consistent with poroelasticity models for fully saturated media. In the unsaturated regime, the pore pressure is equal to the volume averaged fluid pressure with a correction accounting for interfacial energy. As the fluid pressure equals the negative capillary pressure, it becomes negative in the unsaturated regime. We stress that although the subsequent analysis has these particular relations as examples, it does not depend on these particular choices for s w , p pore and κ rel .
Under the hypothesis of small perturbations of the porosity, which are often applied along with the assumptions of linear elasticity [17], we can assume that the porosity φ, acting as weight, is constant in time, equal to some reference porosity field φ 0 . With this, we also note that in the fully saturated regime, i.e., when s w ≡ 1 and, equivalently, p w ≥ 0, the model equations reduce to the classical quasi-static Biot equations.
From now on, we consider a compact form of (2.1)-(2.2). Specifically, we seek (u, In order to close the system (2.5), we consider non-overlapping partitions Furthermore, the pressure satisfies the initial condition Finally, for the initial displacement, u(0), we only consider compatible states, satisfying This also implicitly defines an initial condition for the volumetric deformation ∇ · u(0). Putting the focus on the non-linear and coupled character of the balance equations, in the subsequent, mathematical analysis, we consider a simplified setting, resulting in particular in a simpler notation. We neglect gravity and consider merely homogeneous essential boundary conditions, i.e., from now on we set g ≡ 0, u D ≡ 0 and p w,D ≡ 0.

The mathematical model under the Kirchhoff transformation
The Kirchhoff transformation defines a new pressure-like variable χ(p w ) = pw 0 κ rel (s w (p)) dp. (2.10) Assuming the constitutive laws satisfy κ rel (s w (p)) > 0, for all p ∈ R, the transformation (2.10) can be inverted. Given constitutive laws for b, s w , κ rel , p pore , as e.g. the equivalent pore pressure (2.4) and the van Genuchten-Mualem relations (2.3), from now on, we consider the Kirchhoff transformation χ as unknown instead of p w , and redefine all functions accordinglŷ Then under the assumption of a homogeneous relative permeability and saturation, the non-linear Biot equations (2.5)-(2.9) reduce to finding (u, χ), satisfying subject to the adapted boundary conditions and the initial conditions with u(0) satisfying −∇ · [2µε(u(0)) + λ∇ · u(0)I − αp pore (p w,0 )I] = f (0). (2.15) We highlight that after applying the Kirchhoff transformation, the transformed system (2.12) remains nonlinear and coupled. The main nonlinearities areb,p pore ,ŝ w . Yet, the diffusion term is now linear.

Main result -existence of a weak solution for the unsaturated poroelasticity model
The main result of this work is the existence of a weak solution for the unsaturated poroelasticity model under the Kirchhoff transformation, cf. Sec. 2.2. It should be noted, that the proof utilizes the convergence of a numerical approximation towards a weak solution and thereby also suggests a numerical scheme. The scheme involves the finite element and the finite volume method.
In the following, we state the main result. This includes the notion of a weak solution, required assumptions and the idea of the proof. The details of the proof are the subject of the remainder of this paper.

Definition of a weak solution
We use the standard notation for L p , Sobolev (H k and W k,p ) and Bochner spaces, together with their inherent norms and scalar products. Let ·, · denote the duality pairing between a dual and its primal space; for L 2 (Ω) this becomes the standard scalar product for scalars, vectors and tensors. Let denote the function spaces corresponding to mechanical displacement and fluid pressure, respectively, incorporating essential boundary conditions. We abbreviate the bilinear form associated to linear elasticity and define · V := a(·, ·) 1/2 , which induces a norm on V due to Korn's inequality. Moreover, we combine the external body and surface sources as elements in V and Q , the duals of V and Q, respectively. Let f ext = (f , σ N ) and h ext = (h, w N ) be defined by Remark 3.2 (Discussion of (W1)). For the constitutive relations (2.3) and the equivalent pore pressure (2.4), (W1) is satisfied if, e.g., χ ∈ L 2 (Q T ), as discussed in more detail in Appendix A.
Remark 3.3. We note that the weak formulation of the initial conditions (W3) of the volumetric deformation allows for a stronger formulation. See Lemma 9.6 for more information. In fact, more regularity is asked for in (W3) than required forŝ w (χ)∂ t ∇ · u ∈ L 2 (0, T ; Q ) to be well-defined; however, the higher regularity will be required for interpreting the initial data properly.
The assumptions on the external load and source terms read: and analogously norms corresponding to Bochner spaces as f ext H 1 (0,T ;V ) , etc.
The assumptions on the initial data read: whereB andB are energies defined aŝ Note, that (A1) implies the existence of a convex C 1 -potentialψ : R → R, such thatb =ψ . In this context,B can be related to the Legendre transform of ξ →ψ(ξ) −ψ(0) composed withb, [1] see also Lemma B.12 -similarly forB. Additionally, the following non-degeneracy conditions are assumed: (ND1) There exists a constant C ND,1 > 0 such that (ND2) There exists a constant C ND,2 > 0 such that (ND3) The bulk modulus K dr = 2µ d + λ is sufficiently large and satisfies Physical interpretation of the non-degeneracy conditions. The condition (ND1) essentially means that the equivalent pore pressure essentially behaves like the transformed pressure. The condition (ND2) states that the pore pressurep pore essentially behaves as the transformed pressure χ; in Appendix A, (ND2) is showed to be satisfied for the equivalent pore pressure involving regularized hydraulic properties. Finally, the condition (ND3) essentially requires the mechanical system to be sufficiently stiff in relation to the system's effective compressibility governed by the compressibilities of the fluid and solid grains, as well as the hydraulic properties. In Appendix A, the conditional relation (ND3) is discussed for the constitutive laws presented in Sec. 2. To conclude, (ND3) can be expected to be satisfied in several practical situations.

Existence of solutions for the unsaturated poroelasticity model
This section is presenting the main result together with the main steps of the proof. We observe that uniqueness is not addressed here. The fully coupled, and nonlinear character of the problem makes it in particular difficult to use any monotonicity arguments. This aspect is left open, for further research.
The main idea of the proof of Theorem 3.4 is to use the Galerkin method in combination with compactness arguments. The main difficulty here is the control over the non-linear coupling terms. For this a regularization approach is used. After all, the proof consists of six steps. In the following, we present the idea of each step. Details are subject of the remainder of the article and will be presented in the six, subsequent sections.
Step 1: Physically meaningful regularization. Applying the Galerkin method along with compactness arguments for the original problem (3.2) is challenging due to the coupling terms. A simple way to control the term ∂ t ∇ · u is to add a suitable regularization term in the mechanics equation (3.2a). As the coupling terms also involve non-linearities in the Kirchhoff pressure, ultimately strong compactness is required. Therefore, we add a coercive term in the flow equation, which allows for controlling the term ∂ t χ. In this way, one can control the coupling terms, and eventually leading to convergence.
From a physical point of view, the regularized model accounts for secondary consolidation and compressible solid grains. In mathematical terms, it reads as follows. For given regularization parameters ζ, η > 0, find (u ζη , χ ζη ) to be the solution to the variational equations , whereb η is a strictly increasing regularization ofb (see (A1) η below for further properties). The next two steps prove that the regularized problem has a weak solution in an analogous sense to Definition 3.1.
Step 2: Discretization in space and time. To obtain a fully discrete counterpart of (3.4), we employ the implicit Euler scheme for the discretization in time, a conforming Galerkin finite element discretization of the elasticity equation, and a finite volume discretization of the fluid flow equation, based on two-point flux approximations [21,22]. We stress that such a combination of discretization methods is common in the context of poromechanics [14,45], in particular in the engineering community as well as the industry. Let N ∈ N and {t n ; n = 0, ..., N } be an equidistant partition of the interval [0, T ] with constant time steps τ = T /N . Furthermore, given an admissible mesh T = {K} K , cf. Definition 5.1, let V h ⊂ V denote a conforming, discrete function space for displacements. Let Q h ⊂ Q denote the discrete function space of piecewise constant functions. Then the discretization for time step n ≥ 1 reads: Problem P hn : Given the solution at the previous time step (u n−1 This specific choice for the discretization will turn out to be crucial due to two reasons: (i) the piecewise constant approximation of the pressure allows for choosing non-linear test functions and thereby cancelling the coupling terms in the analysis; (ii) the two-point flux approximation encoded by the discrete gradients ∇ h retain the local character of the differential operator ∇. These together allow for simultaneously cancelling the coupling terms and utilizing the coercivity of the diffusion term, which will be crucial for proving existence of discrete solution via a corollary of Brouwer's fixed point theorem, and for deriving stability estimates.
Step 3: Existence of a weak solution to the regularized model. Based on the discrete values {(u n h , χ n h )} n , we define suitable interpolations in time, (u hτ , χ hτ ), yielding approximations of (u ζη , χ ζη ). We remark that various interpolations are in fact introduced in the course of step 3 and 4. To avoid an excess in notations and for the ease of the presentation, we use the same notation, (u hτ , χ hτ ), for all interpolations throughout this section.
• Weak convergence (up to a subsequence) of the discrete gradients ∇ h χ hτ towards ∇χ ζη is not an obvious consequence of uniform stability. For this, we apply techniques from the finite volume literature [21,37].
Motivated by that, we first derive stability estimates that are uniform wrt. the discretization parameters for some constant C ζη > 0 independent of h, τ -as already indicated in Step 2, the specific spatial discretization is beneficial for obtaining this result. Therefore, one obtains weak convergence for subsequences (denoted the same as before) for h, τ → 0 Moreover, by employing finite volume techniques the following convergence of the discrete diffusion term can be showed for arbitrary discrete test functions q h , which strongly converge towards continuous functions q. Here, ·, · κ abs denotes suitably defined κ abs -weighted scalar products. Finally, the limit, (u ζη , χ ζη ), can be identified as weak solution of the regularized problem (3.4).
Step 4: Increased regularity for the weak solution of the regularized model. When discussing the limit ζ → 0 in step 5, it will be beneficial to have access to the derivative in time of the mechanics equation (3.4a). Under the additional non-degeneracy condition (ND2), stating thatp pore is Lipschitz continuous, an increased regularity can be showed for the weak solution of the regularized model, (u ζη , χ ζη ). For instance, for all The proof follows the same line of argumentation as step 3. First a fully discrete counterpart of (3.6) is constructed by considering differences of (3.5a) between subsequent time steps In addition, suitable interpolationsû t,hτ andp pore,hτ of the discrete values {τ −1 (u n h −u n−1 h )} n and {p pore (χ n h )} n , respectively, define approximations of ∂ t u ζη andp pore (χ ζη ). The uniform stability estimate up to subsequences, for h, τ → 0. Finally, one can identify (3.6) in the limit.
Step 5: Vanishing regularization in the mechanics equation. For each ζ, η > 0, there exists a solution (u ζη , χ ζη ) to (3.4). For the discussion of the limit ζ → 0, we employ compactness arguments similar to step 3. We derive the uniform stability estimates and In order to derive (3.7), we take inspiration from the analysis of the linear Biot equations in [28] and utilize the differentiated-in-time momentum equation (3.6). We use v = ∂ t u ζη as test function in (3.6) (essentially generating ∂ t u ζη H 1 (0,T ;V ) ), and q = ∂ t χ ζη in the fluid flow equation (3.4b), tested with (allowing for simple discussion of the transformed diffusion term). This generates mostly positive terms, which in principle would lead to (3.7). However, unlike in the linear Biot case, the coupling terms do not cancel, and leave behind a nonpositive term. The main idea to recover (3.7) is then to compensate the non-positive term, under a condition on the data and constitutive laws. For this, we apply the mean inequality v V ≥ K dr ∇ · v L 2 , the binomial identity (App. B.2), together with the non-degeneracy condition (ND3) to obtain This intermediate calculation allows to drop the coupling terms in the analysis and retrieve the uniform bound (3.7) at the cost of (ND3). This calculation carefully demonstrates the interpretation of uniform compressibility provided by (ND3), cf. the physical interpretation discussed in Sec. 3.2.
With this, letting ζ → 0, one obtains for subsequences (denoted the same as before) Finally, it is straightforward to see that the limit (u η , χ η ) is weak solution of (3.4) for ζ = 0.
We underline, that for showing (3.9), the time-continuous character of the variational problem is required. It is not obvious how to use a similar strategy on time-discrete level. Therefore, step 5 has been performed separately from step 3 and 4.
Step 6: Vanishing regularization in the flow equation. In the presence of fluid or solid grain compressibility in the original formulation, i.e., c w > 0 or 1 N > 0, respectively, this final step is obsolete. Otherwise, we consider the limit process η → 0 for the sequence of solutions {(u η , χ η )} η , derived in step 5. The overall idea is the same as in step 5, namely to obtain estimates that are uniform wrt. η and to use compactness arguments. Referring to (3.7), the following estimate is uniform in η (3.10) For estimating ∂ t χ η , we first show that the time derivative of the mechanics equation (3.5a) is well-defined for Since ∂ t χ η ∂ tppore (χ η ) , the uniform stability for ∂ t χ η follows by an inf-sup argument, (3.11), and the stability bound (3.10). Due to the lack of a suitable bound on ∂ tt u ζη in step 5, this approach only works for ζ = 0. Standard compactness arguments allow for extracting subsequences (again denoted as before) such that for η → 0 it holds that Ultimately, (u, χ) can be identified as a weak solution to the unsaturated poroelasticity model in the sense of Definition 3.1. This finishes the proof of Theorem 3.4.

4.
Step 1: Physical regularization -secondary consolidation and enhanced grain compressibility We introduce a physical-based regularization of the weak formulation (3.2) by enhancing the mechanics and the flow equations. Specifically, we let ζ > 0 and η > 0 be two regularization parameters. We include secondary consolidation, which effectively adds a linear viscoelastic contribution in the mechanics equation of the form ζa(∂ t u, v). Additionally, we assume non-vanishing solid grain compressibility by defining the regularizationb η ofb asb i.e.,b η has the same structure asb, but with 1 N + η replacing 1 N refering to the physical example (2.6). Refering to Sec. 3.2, the functionb η still satisfies (A1). Additionally, now the following growth condition holds cf. also Appendix A. In the subsequent discussion, a growth condition forb η (orb) of type (A1) η will be required in order to utilize strong compactness arguments for the pressure variable. Note, that if min c w , 1 N > 0 in (2.6), the growth condition (A1) η is fulfilled even for η = 0, and the regularization of the flow equation actually is not necessary, cf.
For a non-degenerate initial condition χ 0 , the additional terms inB η andB η can be essentially bounded by η χ 0 2 L 2 (Ω) , which itself is bounded by (A8). Finally, we introduce the doubly regularized unsaturated poroelasticity model by defining the notion of a weak solution.
Furthermore, we call (u ζη , χ ζη ) a weak solution for the doubly regularized unsaturated poroelasticity model with increased regularity if it additionally satisfies (W1) ζη -(W4) ζη and: We will later separately consider ζ → 0 and η → 0. Therefore, we give the definition of a weak solution for the simply regularized unsaturated poroelasticity model, obtained for η > 0 and ζ = 0. To distinguish between the equations satisfied by the weak solution of a doubly regularized model and the one of the simply regularized one, where ζ = 0, we use the notation (W1) η -(W4) η . Proof. The assertion follows from steps 2-3.
Remark 4.4. The proof of Lemma 4.3 presented in Sec. 6 does in fact allow for replacing the regularizing term ζa(∂ t u, v) in (4.4a) with a reduced regularization merely applied to the temporal derivative of the volumetric deformation, i.e., ζ ∂ t ∇ · u, ∇ · v . The statement of Lemma 4.3 remains true with an adapted version of (W3) ζη addressing only the initial condition for the volumetric deformation. Such regularized model is of wider interest in the literature [17,18,24,31,38], in particular in the context of biomedical applications [5].  Proof. The assertion follows from step 5.

5.
Step 2: Implicit Euler non-linear FEM-TPFA discretization The next two sections, identified with steps 2 and 3, are providing the proof of Lemma 4.3. To this aim, we employ a discretization in space and time. We apply the implicit Euler time stepping method, combined with a conforming Galerkin finite element method for the mechanics equation (4.4a) and a cell-centered finite volume method utilizing a two point flux approximation (TPFA) for the flow equation (4.4b). In this section, we establish the existence of a fully discrete solution. We start with introducing the notations used in the discretization.

Finite volume and finite element notation
We use standard notations in the finite volume literature, cf., e.g., [21,37]. In particular, we introduce notation for elements, faces, their measures, transmissibilities etc. We assume that the domain Ω is polygonal such that it can be discretized by an admissible mesh, as introduced in [22].
Definition 5.1 (Admissible mesh T ). Let T be a regular mesh of Ω with mesh size h, consisting of simplices in 2D or 3D, or convex quadrilaterals in 2D and convex hexahedrals in 3D. Furthermore, we introduce the following terminology: • K ∈ T denotes a single element.
• N (K) := L ∈ T | L = K,L ∩K = ∅ denotes the set of neighboring elements of K ∈ T .
• E denotes the set of all faces, i.e., boundaries of all elements; let E K denote the faces of a single element K ∈ T ; let E ext denote the faces lying on the boundary ∂Ω. • K|L ∈ E denotes the face between two neighboring elements K, L ∈ T .
• {x K } K∈T is such that for all K ∈ T , L ∈ N (K) the connecting line between x K and x L is perpendicular to K|L . • d K,σ denotes the distance between center of K and σ ∈ E K ; • τ σ = |σ|/d σ denotes the transmissibility through σ ∈ E. The additional regularity property is assumed. There exists a constant C > 0 such that In addition, we introduce a dual grid T with diamonds as elements. It will be used for the approximation of heterogeneous permeability fields. Additionally, it will be utilized within the proof, to define suitable projection operators.
Definition 5.2 (Dual grid to T ). Let T be an admissible mesh, cf. Definition 5.1. For each face K|L ∈ E, K ∈ T , L ∈ N (K), define a prism P K|L ⊂ Ω with x K , x L and the vertices of K|L as vertices. For all σ ∈ E ext ∩ E K , K ∈ T define P σ ⊂ Ω to be the prism with x K and the vertices of σ as vertices. By construction, T := {P σ } σ∈E defines a partition of Ω. Fig. 1 displays a two-dimensional admissible mesh and its auxiliary dual grid. The final discrete scheme is written in variational form. Given an admissible mesh T , we introduce the discrete function spaces and implicitly their bases providing spaces for the discrete displacement and pressure, respectively. For the analysis below, we assume that the discrete function spaces satisfy the following conditions: (D1) Q h is the space of all piecewise constant functions (P 0 ) on T and the basis {q h,j } j is equal to the indicator functions of all single elements.
In the analysis, (D1) will allow for intuitively handling non-linearities in the pressure variable. Assumption (D2) will allow for using standard inf-sup arguments. In two dimensions, e.g., piecewise quadratic elements can be used for V h [7]; alternatively, both in two and three dimensions, piecewise linear elements enhanced by face bubbles result in a fairly cheap choice, in particular when utilizing localization techniques [35]. In order to define the discrete scheme, we use a discrete H 1 0 (Ω) inner product and the corresponding norm, as introduced in [21], incorporating averaging of the mobility field.
(Ω) inner product and norm on Q h ). We define the inner product for any χ h , q h ∈ Q h , where the the weight ω evaluated at faces is approximated as weighted average incorporating the neighboring elements {ω} σ := 1 |Pσ| Pσ ω(x) dx, σ ∈ E, utilizing the dual mesh T . In addition, we define the induced norm · 1,T ,ω := ∇ h ·, ∇ h · 1/2 ω , as well as the special case for Remark 5.4 (Two-point flux approximation). Let K ∈ T be a single element. When choosing q h = 1 K as the indicator function of the element K, defined as 1 K (x) = 1 if x ∈ K and 1 K (x) = 0 otherwise, simple calculations yield Hence, a standard two point flux approximation with averaged transmissibilities t K|L = 2|σ| {ω}σ dσ for neighboring cells K, L with common face σ = K|L; and t K,σ = for boundary faces σ.

Approximation of source terms and initial conditions
Let 0 = t 0 < t 1 < ... < t N = T define a partition of the time interval (0, T ) with constant time step size τ = t n − t n−1 , n ∈ N. We interpolate the source terms at discrete time steps. Let Discrete initial conditions are chosen to imitate the compatibility assumption (3.1). Let χ 0 h ∈ Q h be defined by the piecewise constant projection of χ 0 , i.e., on K ∈ T , we define As χ 0 ∈ L 2 (Ω), cf. (A8) η , it follows by classical approximation theory for h → 0 and it holds that χ 0 h 1,T ,κ abs ≤ C χ 0 1 for some constant C > 0, cf., e.g., [22]. Furthermore, sincep pore ∈ C(R), . Furthermore, we define the initial approximate displacement u 0 h ∈ V h to satisfy the compatibility condition Using standard finite element techniques and the convergence of χ 0 h it follows as h → 0 with u 0 defined in (W3). All in all, due to the convergence, (A8) η also applies on discrete level.

Approximation of the evolutionary problem
The discretization of (4.4) is defined by the Galerkin method combined with the standard implicit Euler time discretization: for n ≥ 1, given (u n−1 for all K ∈ T . and with the transmissibilities t K|L and t K,σ as introduced in Remark 5.4. For the remaining discussion the variational formulation (5.5b) is more convenient and is therefore used. Proof. The proof is by induction. We present only the general step, since the proof for n = 1 is similar. We employ a corollary of Brouwer's fixed point theorem, cf. Lemma B.4, to show the existence of a solution of a non-linear algebraic system, which is equivalent to (5.5).
Introduction of a pressure-reduced algebraic problem. We introduce an isomorphism between the discrete function space corresponding to the fluid pressure χ and a suitable coefficient vector space Due to (A4), χ h is well-defined. Similarly, let For given β ∈ R dQ , define α = α(β) ∈ R dV to be the unique solution to Finally, we define F : R dQ → R dQ component-wise by We note, the existence of a discrete solution of Eq. (5.5) is equivalent to the existence of β ∈ R dQ , satisfying F (β) = 0. To prove the existence of a zero of F , we employ Lemma B.4; we consider the expression where we dropped the explicit dependence of α on β and used the identity dQ j=1 β j q h,j =p pore (χ h (β)) s w (χ h (β)) .
We discuss the terms T 1 , T 2 , T 3 , T 4 separately.
Discussion of T 2 . The coupling term T 2 can be reformulated and estimated by employing . Under the use of the binomial identity (App. B.2), the Cauchy-Schwarz inequality and Young's inequality, the coupling term T 2 can be bounded as Discussion of T 3 . By the mean value theorem and (A4), the diffusion term T 3 can be estimated from below Here, the specific definition of the two-point flux approximation is crucial.
Discussion of T 4 . Employing the definition of h ext = (h, w N ), the non-degeneracy condition (ND1), a discrete trace inequality (introducing C tr ), cf. Lemma B.2, together with a discrete Poincaré inequality (introducing C Ω,DP ), cf. Lemma B.1, we obtain h n ext ,p pore(χh (β)) for a constant C (C ND,1 , C tr , C Ω,DP ) > 0 Hence, by (A6) and Young's inequality, for the term T 4 it holds that Combination of all results. By inserting the estimates for T 1 , T 2 , T 3 , and T 4 , (5.8) becomes Finally, since · 1,T ,κ abs defines a norm on Q h and (5.6) holds by induction for n − 1 if n ≥ 2 or from (A8) η for n = 1, by Lemma B.4 (a corollary of Brouwer's fixed point theorem), there exists a β ∈ R dQ such that F (β) = 0. This implies the existence of a discrete solution. The bound (5.6) for n follows immediately from (5.9) due to (A7).

6.
Step 3: Limit h, τ → 0 In the following, we show that the fully-discrete FEM-TPFA discretization, introduced in the previous section, converges (up to subsequence) to a weak solution of the doubly regularized unsaturated poroelasticity model, i.e., we prove Lemma 4.3. For this, we employ standard compactness arguments. Throughout the entire section, we assume (A0)-(A8) and (ND1) hold true. 6.1. Stability estimates for the fully-discrete approximation Proof. The proof follows essentially the same steps as the proof of Lemma 5.6. Therefore, we are quick on similar steps. We consider the reduced displacement-pressure formulation (5.5). We choose v h = u n h − u n−1 h and q h =p pore(χ n h ) sw(χ n h ) as test functions and sum the two equations; note that the second is well-defined asŝ w (χ) > 0 for all χ ∈ R, by (A2). We obtain .
On the left hand side, we employ the binomial identity(App. B.2), Lemma B.12 for the energy termB η defined in (4.3), and the uniform increase ofp porê sw , cf., (A4). It holds that .
Later on, we sum over the time steps 1 to N and utilize a telescope sum. Let us first separately discuss the sums of the two terms on the right hand side. For the first of them, we employ summation by parts, cf. Lemma B.6, as well as the Cauchy-Schwarz inequality and Young's inequality. We obtain The second term is estimated as in the discussion of T 4 within the proof of Lemma 5.6. We obtain Summing (6.1) over time steps 1 to N , inserting (6.2) and (6.3), and rearranging terms, yields Finally, the last term on the right hand side can be controlled after applying a discrete Grönwall inequality, cf. Lemma B.7, using that 2τ < 1 4 . The thesis follows from the assumptions on the regularity of the source terms (A7) (together with a Sobolev embedding) and initial data (A8) η,h .
where b χ,m is from the growing condition (A1 ), C (1) is the stability constant from Lemma 6.1, and C 0 is the bound in (A8 ) h .
Proof. We test (5.5b) with q h = χ n h − χ n−1 h . By using the binomial identity (App. B.2), for the diffusion term, we obtain Dividing by τ and summing over time steps 1 to N , yields Employing the growth condition (A1) η for the first term on the left hand side, and the Cauchy-Schwarz inequality and Young's inequality for the last two terms on the right hand side, yields Finally, the first two terms on the right hand side are bounded by (A8 ) h and (A7), whereas the last term can be bounded by Lemma 6.1. On the left hand side, we employ (A6).
where C 0 and C (1) are the constants from (A8 ) h and Lemma 6.1, respectively.
Proof. Testing (5.5b) with q h = χ n h and employing Lemma B.12, yields for all n For the first term on the right hand side, we employ a similar bound as in the discussion of T 4 within the proof of Lemma 5.6. For the second term, we employ the Cauchy-Schwarz inequality, a discrete Poincaré inequality (introducing C Ω,DP ), and (A6). We obtain Finally, summing over time steps 1 to N and employing Lemma 6.1 and (A7) proves the assertion.
Proof. We utilize a standard inf-sup argument (introducing C Ω,is ), cf. Lemma B.11. Due to (D2), there exists a v h ∈ V h such that Employing the mechanics equation (4.4a), we obtain Finally, the assertion follows from Lemma 6.1, assuming wlog. ζ is bounded from above.
Lemma 6.5 (Stability for the temporal change ofb). There exists a constant C where C (1) is the stability constant from Lemma 6.1.
be an arbitrary sequence of test functions. Employ q n h as test function in (5.5b). Summing over time steps 1 to N and applying the Cauchy-Schwarz inequality, yields For the last term, we employed a discrete trace inequality (introducing C tr ), cf. Lemma B.2, and a discrete Poincaré inequality (introducing C Ω,DP ), cf. Lemma B.1. Finally, after utilizing a discrete Poincaré inequality for the first term on the right hand side, (A6) on the second term, the assertion follows by Lemma 6.1 with

Stability estimates for interpolants in time
Utilizing the discrete-in-time approximations (u n h , χ n h ) n , defined by (5.5), we define continuous-in-time approximations on (0, T ] by piecewise constant interpolation u hτ (t) := u n h , t ∈ (t n−1 , t n ], χ hτ (t) := χ n h , t ∈ (t n−1 , t n ], and by piecewise linear interpolation We deduce stability for the interpolants from the stability of the fully discrete approximation. Lemma 6.6 (Stability estimate for time interpolants of the mechanical displacement). For all h, τ > 0 and τ ∈ [0, τ ) it holds that where C (1) is the stability constant from Lemma 6.1.

Proof. Application of the definition of the interpolation operators and simple calculations yield
Hence, the assertion follows directly from Lemma 6.1.
Analogously, we conclude stability for the interpolants of the Kirchhoff pressure.

Lemma 6.7 (Stability estimate for time interpolants of the Kirchhoff pressure). For all
ζητ , (6.6c) where C (1) and C (2) ζη are the stability constants from Lemma 6.1 and Lemma 6.2, respectively, and C Ω,DP is the discrete Poincaré constant, cf. Lemma B.1.
Proof. The proof is analogous to the proof Lemma 6.6. For (6.6c) and (6.6d), a discrete Poincaré inequality, cf. Lemma B.1, has to be applied before utilizing the stability bound on Similarly, by definition of the piecewise constant interpolation, we deduce stability for some of the nonlinearities used in the model. Lemma 6.8 (Stability estimates for non-linearities evaluated in interpolants). It holds that and C (4) are the bounds from Lemma 6.3 and Lemma 6.4, resp. Lemma 6.9 (Stability estimate for the temporal change ofb). Let It holds that is the bound from Lemma 6.5, and C Ω,P is a Poincaré constant.
Proof. Let q ∈ L 2 (0, T ; Q). We define a piecewise constant interpolation in both space and time, and only time. Then by Lemma 6.5 it holds that Exploiting the stability of discrete gradients, cf. Lemma B.3, and a (continuous) Poincaré inequality (combined introducing C Ω,P ), the triangle inequality and the Cauchy-Schwarz inequality, we obtain which concludes the proof.

Relative (weak) compactness for the limit h, τ → 0
We utilize the stability results from the previous section to conclude relative compactness. We deduce limits for the interpolants which eventually converge towards a weak solution of the doubly regularized unsaturated poroelasticity model, i.e., it fulfills (W1) ζη -(W4) ζη . Lemma 6.10 (Convergence of the mechanical displacement). There exist subsequences of {ū hτ } h,τ and {û hτ } h,τ , denoted by the same subscript, and u ζη ∈ L ∞ (0, Proof. The assertion essentially follows from the stability result in Lemma 6.6. In particular, (6.8a) follows together with the Eberlein-Šmulian theorem, cf. Lemma B.8. For (6.8b), we employ a relaxed Aubin-Lions-Simon type compactness result for Bochner spaces, cf. Lemma B.9. Furthermore, by the Eberlein-Šmulian theorem, there exists aû ∈ L 2 (0, T, V ) such that up to a subsequencê We can identifyû = u ζη as follows. Employing the triangle inequality and Lemma 6.6, yields which converges to zero for h, τ → 0. This proves (6.8c) and (6.8d).
In order to discuss the limit of the pressure, we utilize techniques employed in the finite volume literature [21,37]. For this, we define a piecewise constant discrete gradient ofχ hτ utilizing the dual grid T , cf. Definition 5.2, where n K|L denotes the outward normal on K|L ∈ E, pointing towards L; and n σ,K denotes the outward normal on σ ∈ E ext ∩ E K , pointing towards K.
For the proof of (6.9b), we refer to Proposition 8.1 in [37], which essentially proves the same assertion under the same assumption on uniform stability for N n=1 τ χ n h 2 1,T , here following from Lemma 6.7. The proof of (6.9c) is standard and follows mainly from the stability results in Lemma 6.7 and the Eberlein-Smulian theorem, cf. Lemma B.8. This concludes the proof.
The main purpose of the double regularization has been the aim to get control over the non-linear coupling terms, and eventually establish convergence. Lemma 6.12 (Convergence of the coupling terms). There exists a subsequence of {χ hτ } h,τ , denoted by the same subscript, such thatp pore (χ hτ ) p pore (χ ζη ) weakly in L 2 (Q T ), (6.10a) Proof. By the Eberlein-Šmulian theorem, cf. Lemma B.8, and Lemma 6.8, there exist a subsequence of {χ hτ } h,τ , denoted by the same index, and ap ∈ L 2 (Q T ) such that p pore (χ hτ ) p weakly in L 2 (Q T ).
Proof. Let v ∈ H 1 (0, T ; V ) with v(T ) = 0. We obtain, using the same calculations as in the proof of Lemma 6.13, with the piecewise linear function defined on t ∈ (t n−1 , t n ], n ∈ {1, ..., Based on the stability and convergence results for {ū hτ } hτ and {u 0 h } h , cf. (5.4), we get (for a subsequence, denoted by the same index) for h, τ → 0 and thereby (W3) ζη . 6.4. Identifying a weak solution for h, τ → 0 Finally, we show the limit (u ζη , χ ζη ), introduced above, is a weak solution of the doubly regularized unsaturated poroelasticity model, cf. Definition 4.1.
Proof. The limit (u ζη , χ ζη ) satisfies (W1) ζη -(W3) ζη by Lemma 6.10, Lemma 6.11 Lemma 6.12, and Lemma 6.14. It remains to show (W4) ζη , i.e., that (u ζη , χ ζη ) satisfies the balance equations (4.4). We first consider sufficiently smooth test functions and then use a density argument. Let (v, q) ∈ L 2 (0, T ; V ∩C ∞ (Ω) d )×L 2 (0, T ; Q∩C ∞ (Ω)). For a given mesh T , we define spatial projection and interpolation operators, respectively, by 14) Similarly, letf ext,τ (t) := f n ext , t ∈ (t n−1 , t n ], h ext,τ (t) := h n ext , t ∈ (t n−1 , t n ]. Combining classical results, based on the assumed regularity (A7), for h, τ → 0 it holds that v hτ → v strongly in L 2 (0, T ; V ), We choose v h = v n h and q h = q n h as test functions in (5.5), multiply both equations with τ and sum over all time steps 1 to N ; we obtain For most terms we can apply the fact that the product of weakly and strongly convergent sequences converge to the product of their limits. The only term needing special attention is the diffusion term in the flow equation. For this, we follow a similar argument as in the proof of Theorem 3.4 in [37] -here for a heterogeneous permeability field. The technical calculations are the same. Hence, we only summarize the idea. We exploit the definition of the piecewise constant gradient ∇χ hτ . Similarly, we define a piecewise constant gradient of the smooth test function and the permeability field ∇q hτ (x, t) := ∇q n (x σ ), (x, t) ∈ P σ × (t n−1 , t n ], σ ∈ E, where x σ ∈ P σ is chosen by the mean value theorem satisfying (here only for an internal edge) 1 d K|L (q n (x K ) − q n (x L )) = ∇q n (x K|L )·n L|K , i.e., connecting the two-point flux approximation with the actual gradient of q. A calculation, just as in [37], yields Due to sufficient regularity, it holds that {κ abs } T → κ abs strongly in L ∞ (Q T ) and ∇q hτ → ∇q strongly in L 2 (Q T ). Since by Lemma 6.11 it also holds that ∇χ hτ ∇χ ζη weakly in L 2 (Q T ), we finally conclude that All in all, together with the convergence properties of the test functionsv hτ ,q hτ , the source termsf ext,τ , h ext,τ , and the interpolants for the fully discrete approximations (cf. Lemma 6.10, Lemma 6.11, Lemma 6.12 and Lemma 6.13), we conclude that (6.18) converges to (4.4), evaluated in (u ζη , χ ζη ) and tested with (v, q) ∈ L 2 (0, T ; V ∩ C ∞ (Ω) d ) × L 2 (0, T ; Q ∩ C ∞ (Ω)). Finally, a density argument yields the final result.

Step 4: Increased regularity in a non-degenerate case
In the following, further stability estimates for the fully-discrete problem are derived, allowing for showing that the limit (u ζη , χ ζη ) introduced in the previous section also satisfies (W5) ζη and (W6) ζη , i.e., we prove Lemma 4.5. For this, non-degeneracy assumptions are required. For compact presentation throughout the entire section, we assume (A0)-(A8) and (ND1)-(ND2) hold true, and we define u −1 h := u 0 h . 7.1. Improved stability estimates for fully-discrete approximation Lemma 7.1 (Improved stability estimate for the structural velocity). There exists a constant C ζη is the stability constant from Lemma 6.2, C ND,2 corresponds to the non-degeneracy condition (ND2), and b χ,m corresponds to the growth condition (A1 ).
Proof. First we observe, that the compatibility condition for the initial conditions (5.3) is equivalent to the mechanics equation (5.5a) for n = 0, since u 0 h − u −1 h = 0. This allows for considering the difference of the mechanics equation (5.5a) at time steps n and n − 1, n ≥ 1, ) and using the binomial identity (App. B.2), we obtain Summing over n ∈ {1, ..., N }, yields after applying the Cauchy-Schwarz inequality and Young's inequality for the right hand side terms Due to (ND2),p pore =p pore (χ) is Lipschitz continuous. Therefore, by Lemma 6.2 it holds that which together with (7.1) concludes the proof.
Lemma 7.2 (Consequence for the structural acceleration). There exists a constant C where C (6) ζη is the stability constant from Lemma 7.1. Proof. Let {v n h } n ⊂ V h \ {0} be an arbitrary sequence of test functions. Consider the difference of (5.5a) at n and n − 1, n ≥ 1, tested with v n h ; it holds that Summing over n ∈ {1, ..., N }, applying the Cauchy-Schwarz inequality and Lemma 7.1, yields ζη .

Improved stability estimates for interpolants in time
We define piecewise linear interpolations of the discrete structural velocity and the pore pressure. For t ∈ (t n−1 , t n ], n ≥ 1,  4) and (7.2a). For all h, τ > 0 andτ ∈ (0, τ ), it holds that where C (6) ζη and C (7) ζη are the stability constants from Lemma 7.1 and Lemma 7.2, respectively, and C Ω,PK is the product of the Poincaré and the Korn constants.
Proof. By construction, it holds that Hence, (7.3a) follows directly from Lemma 7.1. The time-translation property (7.3b) follows from the fact that ∂ tûhτ is piecewise constant. Analogous to the proof of Lemma 6.6, one can show Finally, after using a Poincaré inequality and Korn's inequality, (7.3b) follows from Lemma 7.1. In order to show (7.3c), we expand the integral over the time interval. By definition ofû t,hτ , it holds that û t,hτ Hence, (7.3c) follows by Lemma 7.1. In order to show (7.3d), we again expand the integral over the time interval.
Employing the definitions ofû t,hτ andû hτ , simple calculation yield Hence, after employing a Poincaré inequality and Korn's inequality, (7.3d) follows from Lemma 7.2. Finally, (7.3e) follows directly from Lemma 7.2, since Lemma 7.4 (Stability result for the interpolation of the pore pressure). Forp pore,hτ defined in (7.2b). It holds that where C (6) ζη is the stability constant from Lemma 7.1. Proof. By construction, it holds that where the second result follows by expanding time integration. Hence, the assertion follows directly from the stability result for the discrete time derivative of the pore pressure, cf. Lemma 7.1.

7.3.
More relative (weak) compactness for h, τ → 0 The previous stability results allow for analyzing the limit in relation to (u ζη , χ ζη ).

Stability estimates independent of ζ
The key ingredients for the subsequent discussion are stability estimates, which are independent of ζ. As a result of the stability results derived in Sec. 6.2, there exists a constant C = C C (1) , C (4) > 0 (independent of ζ > 0 and η > 0), such that where Lemma 6.6 and Lemma 6.10 yield stability for the displacement, and Lemma 6.8 and Lemma 6.12 yield stability for the pore pressure. Further stability bounds can be obtained by exploiting the continuous nature of the balance equations and the time derivative of the mechanics equation. The following result is crucial.
We discuss the individual terms separately. For the first two terms on the left hand side of (8.2), we employ the fundamental theorem of calculus where we used that ∂ t u ζη (0) = 0. This follows from the temporal derivative of the mechanics equation (4.5) and the compatibility condition for the initial conditions (3.1).
For the remaining terms on the left hand side of (8.2), we apply the binomial identity (App. B.2), employ that b η ≥b ≥ 0 (cf. (A1)), condition (ND3), and For the first term on the right hand side of (8.2), we apply the Cauchy-Schwarz and Young's inequality For the second term on the right hand side of (8.2), we apply integration by parts, a Cauchy-Schwarz inequality and Young's inequality, a Poincaré inequality and a Sobolev embedding, involving the constants C Ω,P and C T,Sob , as well as (A6). This gives Altogether, (8.2) becomes A Grönwall inequality, and the general choice ofT , yields the assertion under the given assumptions.
The previous stability estimate allows for deriving revising stability estimates in the previous section.
Lemma 8.2 (Stability for the energyB η ). There exists a constant C (9) > 0 (independent of ζ, η), such that where C (8) is the stability constant from Lemma 8.1, and C 0 is the stability constant from (A8 ).
Proof. Testing the flow equation (4.4b) with q = χ ζη , yields For the first term on the left hand side, we apply a property ofB η , cf. Lemma B.12, .
Proof. The proof is similar to the proof of Lemma 6.5. However, this time, we exploit We will require to show strong convergence of the Kirchhoff pressure. Having that in mind, we conclude with a stability estimate for ∂ t χ ζη . We note, this is the only stability estimate in this section, requiring the regularizing growth condition (A1) η . Lemma 8.4 (Stability estimate for the temporal change of the Kirchhoff pressure). There exists a constant C where C (8) is the stability constant from Lemma 8.1, b χ,m is from (A1 ), and C 0 is from (A8 ).
Proof. We repeat parts of the proof of Lemma 8.1. We test the flow equation (4.4b) with q = ∂ t χ ζη and apply (A1) η and the Cauchy-Schwarz inequality; we obtain After rearranging terms, applying the regularity of the data, and applying Lemma 8.1, the assertion follows.

Relative (weak) compactness for ζ → 0
We utilize the stability results from the previous section to conclude relative compactness.
Lemma 8.7 (Initial conditions for the fluid flow). Up to subsequences it holds for ζ → 0 where ∂ tbη (χ η ) is understood in the sense of (W2) η .

9.
Step 6: Limit η → 0 in the incompressible case In this section, we show the main result, Theorem 3.4, for the more demanding case: the presence of an incompressible fluid and incompressible solid grains. Otherwise, by Remark 8.10 the main result of this paper follows already. In the incompressible case, b as in (2.6) is monotone but withb = 0 on a part of the domain with non-zero measure. Under the use of regularization with η > 0, it holds that b χ,m = η. In the following, we prove that the limit of {(u η , χ η )} η exists for η → 0, and that it is a weak solution of (2.12)-(2.15) according to Definition 3.1. Throughout the entire section, we assume (A0)-(A8) and (ND1)-(ND3) hold true.
numerical convergence, consistent with the above analysis. Further study of accuracy properties, as well as existence and uniqueness, of the discretization of the unregularized model are planned in the future.
Appendix A. Practical set of constitutive laws In this section, we consider the constitutive laws for s w , κ rel and p pore , as presented in Sec. 2.1, and show that in this case that the existence theory presented in this work is applicable. In particular, we consider  where · 1,T denotes the discrete H 1 0 (Ω) norm, introduced in Definition 5.3. Lemma B.3 (Stability of discrete gradients [21]). Let T be an admissible mesh of some domain Ω, cf. Definition 5.1, and u ∈ H 1 0 (Ω). Define a piecewise constant functionũ bỹ u(x) := 1 |K| K u(x) dx, x ∈ K ∈ T .
Then there exists a constant C > 0 (independent of h for regular meshes) such that ũ 1,T ≤ C u H 1 (Ω) .
Lemma B.4 (Corollary of Brouwer's fixed point theorem [15]). Let ·, · denote the standard R d scalar product and let F : R d → R d be a continuous function, satisfying for all x ∈ R d with x, x ≥ M for some fixed M ∈ R + . Then there exists a x ∈ R d with x , x ≤ M and F (x ) = 0. b n (a n+1 − a n ).
Lemma B.7 (Discrete Grönwall inequality [16]). Let (a n ) n ⊂ R + , (λ n ) n ⊂ R + , B ≥ 0. Assume for all n ∈ N it holds that a n ≤ B + n−1 k=0 λ k a k . Then it follows a n ≤ B n−1 k=0 (1 + λ k ). In particular, if λ k = λT N for all k ∈ N for some λ, T ∈ R + and N ∈ N, it holds that a N ≤ B exp(λT ).
Lemma B.8 (Eberlein-Šmulian theorem [15]). Assume that B is a reflexive Banach space and let {x n } n ⊂ B be a bounded sequence in B. Then there exists a subsequence {x n k } k that converges weakly in B.
Lemma B.9 (Relaxed Aubin-Lions lemma [41]). Let {f n } n ⊂ L p (0, T ; B), 1 ≤ p < ∞, B a Banach space. {f n } n is relatively compact in L p (0, T ; B) if the following two are fulfilled: • {f n } n is uniformly bounded in L p (0, T ; X), for X ⊂ B with compact embedding. • Lemma B.10 (Riesz-Frechet-Kolmogorov compactness criterion [11]). Let F be a bounded set in L p (R N ) with 1 ≤ p < ∞, N ∈ N. Assume that Then the closure of F | Ω := {f : Ω → R | f ∈ F } is compact for any measurable set Ω ⊂ R N with finite measure.
Lemma B.11 (Standard inf-sup argument [7]). Let V and Q be Hilbert spaces, and let B be a linear continuous operator from V to Q . Denote by B t the transposed operator of B. Then, the following two statements are equivalent: • B t is bounding, i.e., there exists a γ > 0 such that B t q V ≥ γ q Q for all q ∈ Q.