Mathematical Theory of Nonlinear Single-Phase Poroelasticity

In this paper, we study the equations of nonlinear poroelasticity derived from mixture theory. They describe the quasi-static mechanical behavior of a fluid saturated porous medium. The nonlinearity arises from the compressibility of the fluid and from the dependence of porosity and permeability on the divergence of the displacement. We point some limitations of the model. In our approach, we discretize the quasi-static formulation in time and first consider the corresponding incremental problem. For this, we prove existence of a solution using Brézis’ theory of pseudo-monotone operators. Generalizing Biot’s free energy to the nonlinear setting, we construct a Lyapunov functional, yielding global stability. This allows us to construct bounds that are uniform with respect to the time step. In the case when dissipative interface effects between the fluid and the solid are taken into account, we consider the continuous time case in the limit when the time step tends to zero. This yields existence of a weak free energy solution.


Introduction
The elastic quasi-static deformation of a fluid saturated porous medium received much attention in the civil engineering literature because of its relevance to many problems of practical interest. In the framework of consolidation in soil mechanics, these problems relate to the physical loading of soil layers or the effect of soil subsidence due to groundwater withdrawal for drinking water supply or industrial and agricultural purposes. Examples and underlying theories are given in the well-known works of Coussy (2004), Lewis and Schrefler (1998) and Verruijt (2015). They build on the classical theory of Terzaghi (1951) and the pioneering approach of Biot (1962) and Tolstoy (1992).
In its simplest form, assuming both the fluid and the porous material (grains) to be incompressible and assuming the porous medium to be homogeneous and linearly elastic with small strains, the mathematical formulation reads (see Bear and Bachmat 1990;Verruijt 2015or van Duijn et al. 2019: and where with G E = 2μE + λTr(E)I, for symmetric matrices E.  (4), i.e., Hooke's law, assumes that the skeleton is mechanically isotropic. The linear quasi-static Biot system, as well as its dynamical analog, was also derived by means of a multiscale approach, where the starting point is the linear fluid-structure interaction at the pore level. We refer to the monographs Sanchez-Palencia (1980) and Mei and Vernescu (2010) for derivations using two-scale expansions and to Mikelić and Wheeler (2012) for a rigorous mathematical derivation by means of homogenization. The derivations using multiscale analysis confirm Biot's models in the linear setting. Hence, from different points of view system (1)-(4) is well accepted.
In the engineering literature, one writes α = 1 − K /K g , where K is the drained bulk modulus of the porous skeleton and K g the bulk modulus of the grains. Since it is assumed that K g = +∞, we will set α = 1 in (3).
From a mathematical perspective, Eqs.
(1)-(4) received much attention. Here, we mention the pioneering paper by Auriault and Sanchez-Palencia (1977) and the work of Ženíšek (1984), who were the first to demonstrate existence and uniqueness. More recent studies include Showalter (2000), Owczarek (2010) and Marciniak-Czochra and Mikelić (2015). Later, Cao et al. (2013) considered a nonlinear extension of (1), by replacing the permeability tensor K by the product Kk(div u). The function k(·) is a relative permeability depending on the volumetric strain div u. From (1), we notice that the overall mixture of two incompressible phases is not incompressible itself.
Though system (1)-(4) is linear, its mathematical complexity lies in the fact that it is of quasi-static nature. In particular (2)-(4) allow to control the size of the volumetric strain only through the size of the data. Some authors circumvent this by introducing a time dependence in (2)-(4) as well. For instance, Bociu et al. (2016) replace u in (3) by u + δ∂ t u, where δ ≥ 0 is a visco-elastic parameter. Their study allows δ = 0, hence it includes the true quasi-static case as well. A different regularization was proposed by Murad and Cushman (1996) who replaced (3) and (4) by σ = 2μe(u) + (λdiv u + λ * div ∂ t u − α p)I, with λ * > 0. This form arises in the non-equilibrium theory, where the fluid pressure and the solid pressure differ by λ * div ∂ t u.
In this paper, we propose to study the quasi-static formulation in which we replace Eq. (1) by the nonlinear fluid phase mass balance based on the mixture theory of Bedford and Drumheller (1978) and Bedford and Drumheller (1983), see, e.g., Rutqvist et al. (2001) and Lewis and Schrefler (1998): where j denotes the Darcy mass flux Here, n denotes porosity, ρ = ρ f [kg/m 3 ] fluid density, k relative permeability and Q [kg/m 3 s] sources/sinks.
An explicit expression for (8) is derived from the Lagrangian solid mass balance equation. This is shown in Sect. 2. Through (8), the relative permeability depends on div u.
Since n is the volume fraction of voids in the porous medium, it should satisfy the natural bounds 0 < n < 1.
However, in Sect. 2 we show by means of a counter example that the porosity can attain negative-and thus physically unrealistic-values. Therefore, the bounds in (9) are a major concern in the mathematical model.
To close system (2)-(4), (6)- (7), we introduce a constitutive relation for the fluid density in terms of the pressure. Assuming weak compressibility, we write Further, we propose an explicit expression for the relative permeability in terms of the porosity In (10), ρ 0 and p 0 are reference values for, respectively, density and pressure and β [Pa −1 ] is the fluid compressibility coefficient. The relative permeability in (11) satisfies A well-known example is the Kozeny-Carman formula, see for instance Bear and Bachmat (1990), in a realistic porosity interval, bounded away from n = 0 and n = 1. Thus, taking k such that (12) holds and for 0 < n * < n < n * < 1, for appropriately chosen 0 < n * < n * < 1, gives a relative permeability satisfying (13) in the interval [n * , n * ].
We notice that Eq. (6), coupled with (2) and (4), is nonlinear due to the relation k = k(n) and the products involving time derivatives. Assuming constant fluid phase density in the poroelastic mixture is therefore an important simplification. This is studied in Cao et al. (2013) and Bociu et al. (2016).
In studying system (2)-(4), (6)-(11), a crucial role is played by its free energy. The idea is to generalize Biot's free energy (Biot 1962), which is quadratic in strain and fluid density, to the nonlinear poroelastic setting. This free energy serves as a Lyapunov functional. This approach is linked to general entropy methods for PDEs. For a detailed survey covering various fields of applications, we refer to Evans (2004) and to the recent book by Jüngel (2016). An interesting application of the entropy method is discussed in Mikelić (2010), Cao and Pop (2016) and Milišić (2018), where the authors consider dynamic capillary pressure effects in two-phase porous media flow.
This paper is organized as follows. In Sect. 2, we present details of the model formulation. The starting point is the mass balance for the fluid and the solid phase. The latter implies an explicit expression for (8). Introducing a lower bound for the porosity, we modify the fluid mass balance so that a Lyapunov functional can be constructed for the modified system. This modification is such that the fluid equation reduces to its original form in the physical range of the fluid density ρ and solid volumetric strain E. Section 2 is concluded by a weak formulation of the modified system.
In Sect. 3, we consider, for the relaxation parameter λ * ≥ 0, the incremental version of the modified system. Using Brézis' theory of pseudo-monotone operators, existence is demonstrated. Applying the Lyapunov functional yields global (in time) estimates. Next, in Sect. 4, we use these estimates to solve the time-continuous problem when λ * > 0. In both Sects. 3 and 4, we borrow ideas from Roubiček (2005). Finally, in Sect. 5 we present a discussion and conclusions.

Problem Formulation
In a number of steps, we construct in this section the equations that serve as starting point for the analysis. The general setting of the problem is as follows: Let Ω ⊂ R m (m=2,3) denote a smooth bounded domain, occupied by a linear elastic skeleton. The skeleton material (grains) is assumed incompressible: i.e., the bulk modulus of the grains is infinitely large. The voids in the porous structure are completely filled with a slightly compressible fluid, in the sense that the fluid pressure p and density ρ are related by (10).

Balance Equations
For given ξ ∈ Ω, let x(ξ, t) denote the location of a solid particle at time t > 0, that started at x(ξ, 0) = ξ,. Then the skeleton velocity v s is given by v s = ∂ t x| ξ .
Restricting themselves to small displacements u (within the elastic regime), Rutqvist et al. (2001) and Lewis and Schrefler (1998) argue that in the mass balance equation for the fluid and solid, the material derivative D Dt = ∂ t + v s · ∇ can be replaced by the partial derivative ∂ t . This is made explicit by a scaling argument in van Duijn et al. (2019). The resulting Lagrangian form of the mass balances reads: and where j is mass flux (7). Within the same approximation, one may write div v s = ∂ t div u.
Using this in (15) and (16) gives and Integrating (18) in time from t = 0, say, to t > 0, we have Here, U 0 is the initial displacement and n 0 the initial porosity. With n 0 ∈ (0, 1) in Ω, expression (19) ensures To avoid technical complications, we restrict ourselves to n 0 = constant in Ω. For small displacements u − U 0 , expression (19) is approximated by Remark 1 Frequently, the linear form (21) is used for values of div u in a neighborhood of div U 0 : i.e., in practical circumstances (21) is applied when E * < div (u−U 0 ) < E * , where E * < 0 < E * are appropriately chosen.
Throughout the paper, we redefine where U 0 ∈ H 1 0 (Ω) m ∩ H 2 (Ω) m is the initial displacement. Redefining accordingly we obtain for the fluid pressure p and the skeleton displacement u the system: where Remark 2 Concerning the initial displacement U 0 , we note that only div U 0 , the initial volumetric strain, is used. However, when discussing the free energy, one needs in addition that U 0 is such that the corresponding elastic energy is finite. For simplicity, we suppose U 0 ∈ H 2 (Ω) m .
In the next sections, we will develop the mathematical theory for system (24)-(28). The issue of negative porosity in (27) (or, for that matter, a porosity exceeding one in approximation (28)), is discussed next.

Negative Porosity
We consider a simplified version of the linear problem (1)-(4) and show that div u can attain values for which the porosity from (27)-(28) becomes negative.
For simplicity, we give the construction in R 2 .
Let Ω = (0, L) 2 for some L > 0. We suppose, as in the rest of this paper, that div u| t=0 = 0. Further we set F = 0 in (25). Using (4) in (25) gives Proceeding as in Verruijt (2015), when he discusses the Mandel problem, we take the divergence of (29) to obtain Hence the function The idea is to prescribe boundary conditions for Eqs. (24) and (25) so that H | ∂Ω is given. For instance, if we set along the four edges, see and use Repeating this along the other edges gives where Σ b denotes the given value of σ along the edges. Then, we have Proposition 1 Let E = divu denote the volumetric stress and let n(E) be given by (27). Suppose there exists a constant Σ > 0 such that Σ b ≤ −Σ. Then, for Σ sufficiently large, there exists a T p = T p (Σ) > 0 such that Proof Note that the sign of Σ b implies compression of the medium. Restricting ourselves to the linear case (1) in a homogeneous and isotropic porous medium in which sources/sinks and gravity are absent, we have Since we have for E = div u the problem By the strong maximum principle, E < E in Ω and for t > 0, where E is the solution of problem (34) with E = −Σ/(2μ + λ) on ∂Ω. Writing E as a Fourier series, one observes that the result is immediate.
This example shows that there is a problem with the model. A modification is needed to prevent the porosity (27), or (28), to become negative. Of course, one could argue that this is outside the scope of the model or outside the range of practical applications, since linear elasticity and small strains are supposed. However, since it is not clear how to ensure that indeed small displacements/strains are guaranteed, one needs to impose a porosity modification to prevent negative values.

Modification of Balance Equations
In a number of steps, we modify Eq. (24) so that it becomes well-posed in a mathematical sense and reduces to its original form in the physical range of the unknowns.
First, to satisfy the natural bounds (9), we replace the porosity approximation (28) by a smooth increasing function n : R → R such that Here, E * and E * are practical values chosen such that −n 0 /(1 − n 0 ) < E * < 0 < E * < 1 and δ 0 = n(E * )/2, see Fig. 2for a sketch. This construction ensures that the modified porosity n(E) remains in the physical range (0, 1) and coincides with the linear approximation in the interval (E * , E * ). Realistic porosity measurements are always done away from the bounds n = 0 and n = 1, see, e.g., Bear and Bachmat (1990).
We choose to study Eq. (24) with the fluid density as primary unknown. Hence, we need to express the pressure p in terms of ρ. Using (26), we have explicitly When considering (24), one clearly has in mind that ρ takes values near the reference ρ 0 . However, the mathematical nature of the equations does not guarantee this behavior. Hence, a second modification is needed, now for ρ in the second and third term of the left-hand side of (24). Disregarding gravity, we replace (24) by the modified fluid mass balance equation where n(E) is given by (35) and where ρ * ∈ (0, ρ 0 ) is a small constant. Outside this range we take for d and D extensions that suit the mathematical analysis. We clarify this at a later point in this section.

Remark 3 The composite function
The balance of forces (25) is modified by adding the regularizing term λ * ∂ t E, as in expression (5). This gives where E = div u and where λ * ≥ 0.
We consider system (37), (39) in the set where T > 0 is arbitrarily chosen. To avoid technical complications, we take ∂Ω ∈ C 1 throughout the rest of this paper. As initial conditions, we have where ρ 0 : Ω → (0, +∞) is taken near the reference value ρ 0 . Along the boundary, we prescribe where ν is the outward unit normal at ∂Ω.

Lyapunov Functional
In this section, we derive an expression for the free energy which acts as a Lyapunov functional for system (37), (39). This generalizes the free energy introduced originally by Biot (1962). Let {u, ρ} be a smooth solution of Eqs. (37), (39) that satisfies conditions (40) and (41). Further, let g : R → R be a smooth, strictly increasing and globally Lipschitz function satisfying g(ρ 0 ) = 0.
We first multiply equation (39) by ∂ t u and integrate the result in Ω. This gives Next, we multiply (37) by g(ρ) and integrate the result in Ω. This results in With the first term in (43) can be written as Note that G is a nonnegative, convex function with G(ρ 0 ) = 0. We substitute (45) back into (43). Adding the resulting expression and (42) yields Before considering the general nonlinear case described by this expression, we first show its implication for the simplified linear setting where we have in expression (46) simplifies to Since This gives Using these expressions in (46) yields Hence, acts as a Lyapunov functional for the linear form of system (37), (39). The first term denotes the elastic energy of the skeleton, the second term the compression energy of the fluid, and the third term the work done by the force F. Expression (50) coincides with Biot's original free energy expression from Biot (1962).
Next, we return to the nonlinear case (46). As a first step, we restrict ourselves to the physical range of the porosity. Then, integral (47) becomes Differentiating the expression yields a first-order equation for g. Thus for (51) to vanish, g should satisfy the initial value problem We first consider this problem in the interval |ρ −ρ 0 | < ρ : Then, (53) reduces to Direct integration results in A second integration yields (Fig. 3) Sketch of the free energy βG(ρ/ρ 0 ). The linear case is in blue. The nonlinear case, see (56) and (59) with n 0 = 1/3 and ρ * /ρ 0 = 0.01, is in black has not yet been defined. We do this by first extending g(ρ) for |ρ − ρ 0 | > ρ and then by solving d(ρ) from (52): i.e., Clearly, (55) cannot be used for ρ ≤ 0. Instead, we extend (55) in a linear C 1 -manner yielding Substituting expressions (58) and (59) in (57) yields the desired extension for d(ρ) when |ρ − ρ 0 | > ρ. Thus, with g and G given by (58)and (59) Hence, the triple {g(ρ), G(ρ), d(ρ)} constructed above satisfies (52). For this choice, the integral (51) drops from expression (46). So far, we considered for the porosity the linear approximation n(E) = n 0 + (1 − n 0 )E. To deal with the full cutoff (35), we introduce a second modification. The starting point is (47). This integral vanishes if Keeping g as in (55), (58) and G as in (56), (59), we now modify d(ρ), calling it Using (57) in this expression gives Clearly, for |ρ − ρ 0 | < ρ and E * < E < E * , this expression reduces to Finally, we use in the Darcy mass flux term j from Eq. (37) Thus, in the end we consider the "second" modified fluid mass balance equation System (39), (65) serves as starting point of the analysis. The function D(ρ, E) in (65) generalizes the fluid density. It is chosen so that (66) acts as a Lyapunov functional for the system. The function G : R → R satisfies G(ρ 0 ) = 0, G(ρ) > 0 if ρ = ρ 0 and G is strictly convex, with quadratic behavior for large values of |ρ|. It is explicitly given by (56) and (59).

Summary of Equations and Weak Formulation
The problem describing the nonlinear poroelastic behavior of a fluid saturated porous medium is to find the displacement u : Q T → R m and the fluid density ρ : The coefficients in Eqs. (67)-(68) were introduced in this section. Specifically, n(E) and k(E) satisfy (35) and Remark 3, D(ρ, E), D(ρ) and p(ρ) are given by (62), (64) and (36), and λ * ≥ 0.
We recast this classical formulation in the following weak form.
In Definition 1, we explicitly incorporate energy inequality (72). When dealing with classical solutions, Eqs. (67)-(68) imply the energy balance (see (46), (47) and (66)) However, in the weak formulation (69)-(70) we cannot use Φ = g(ρ) and ξ = ∂ t u, due to lack of smoothness. Therefore, (v) has to be added explicitly. Hence, we consider only those weak solutions satisfying additionally (72). Therefore, they are called weak free energy solutions. In a number of steps, we prove existence of weak solutions when λ * > 0. We achieve this by first considering the incremental formulation. In this approximation, which is clearly relevant when treating the problem numerically, we obtain existence results which hold for all λ * ≥ 0.

Existence of a Solution to the Incremental Problem
In this section, we study the time discretized form of (67), (68).
In doing so we use the function g = g(ρ), defined by (55) and (58) as the primary unknown. This is allowed since g : R → R is smooth and strictly increasing. The switch to g is done for mathematical convenience, because it allows us to obtain Lyapunov functional estimates in a straightforward way. We start with some definitions. Let p(g) := p(ρ(g)) and D(g) := D(ρ(g))ρ (g).
Further, since  (62), Note that the first term in D(g, E) is bounded with respect to E and grows linearly in g for large |g|. The second (pressure) term is bounded with respect to g since Using these definitions in (67) and (68), we find in terms of g in Q T . Next, we turn to the time discretized form of Eqs. (77) and (78). Let τ ∈ (0, 1) denote the time discretization step and N ∈ N a large integer such that N τ = T . At each discrete time t j = jτ , with j = 0, 1, . . . , N , we set Let u j−1 and g j−1 denote, respectively, the displacement and transformed density at t j−1 for some j ∈ {1, 2, . . . , N }: i.e., Then, u and g at time t j are obtained as solutions of the incremental problem (writing The coefficient D τ in Eq. (79) is given by This expression results from D(g, E) in (76), when the derivative n (E) is replaced by the finite difference n( div u) − n( div U) div u − div U . The specific choice of (81) appears convenient in the estimates concerning the time-discrete Lyapunov functional.
Using the weak topology of the space H 1 0 (Ω) m × H 1 (Ω), serious difficulties arise with the coefficients n, D τ and k depending on div u. To remedy this, we introduce a mollifier Υ ε , where ε is a small positive parameter (see, e.g., Roubiček 2005, page 203), and replace div u in the nonlinearities by the convolution div u Υ ε = u ∇Υ ε . Using this substitution one can treat nonlinear coefficients containing div u as lower order terms in the equations. This allows us to use the theory of pseudo-monotone operators. Applying this convolution, the regularized form of Problem (PD) reads: Note that the denominator τ g ε in (82) originates from the time step τ in the discretization and the term g in the denominator of (76). We have the following existence result Proposition 2 Let ε > 0 be a small positive constant. Under the assumptions of Definition 1, Problem (PD) ε admits at least one solution (u ε , g ε ) ∈ V .
Proof We start by introducing a nonlinear operator A, defined on V and with values in its dual V . It results from adding (82) and (83). We write the resulting relation, with (u, g) ∈ V , as where and The idea is to show that A is a perturbed monotone operator: i.e., A is monotone in its principal part containing derivatives of u and g. To be precise, we show that A is pseudo-monotone and coercive. This allows to apply Brézis' theorem to (84) (see Chapter 2 in monographs Lions 1969;Roubiček 2005 or Chapters 26 and 27 in Zeidler (1990)) to conclude existence for Problem (PD) ε . For the comfort of the reader, we recall that an operator A : V → V is pseudomonotone if and only if A is bounded and The boundedness of A is immediate. To show (87), we follow Chapter 2 from Roubiček (2005) or Chapter 17 from Schweizer (2018) and rewrite A in a form having a principal part containing partial derivatives of u (in e(u) and div u) and ∇g, and a lower order part containing u and g. Specifically, we introduce the operator B : We observe that B (u, g), (u, g) = A(u, g). The introduction of B is useful because it reflects the monotonicity of the principal part of A(u, g). This is a direct consequence of B (w, ), (u 1 , g 1 ) − B (w, ), (u 2 , g 2 ) , (u 1 , g 1 ) − (u 2 , g 2 ) ≥ 0, with equality if and only if u 1 = u 2 and g 1 = g 2 . Inequality (89) is checked by a short computation in (88).
To show (87) we consider a sequence {u r , g r } ⊂ V such that (u r , g r ) (u, g) weakly in V and lim sup r →+∞ A(u r , g r ), (u r , g r ) − (u, g) ≤ 0.
As in Roubiček (2005) we set (u δ , Using the monotonicity from (89), we obtain The sequence (u r , g r ) is bounded in V and there exists a subsequence which strongly converges in L 5 (Ω) m and (a.e.) in Ω, to (u, g). Hence, it suffices to pass to the limit along this subsequence. In (91), the terms containing the operator B are fixed with respect to the gradients. Hence, lim r →+∞ B (u r , g r ), (u δ , g δ ) , (u r , g r ) − (u, g) = 0, and lim r →+∞ for any (ξ, ψ) ∈ V . With these results, we are in a position to pass to the limit r → +∞ in inequality (91). It yields By the pseudo-monotonicity hypothesis (90), inequality (94) implies We use this inequality to conclude This completes the proof of the pseudo-monotonicity.
It remains to prove coercivity. We evaluate directly the term Taking ξ = u − U and ψ = g in (85), the cross terms involving the product p(g) div (u − U) cancel and the term p 0 div (u − U)/τ drops out after integration. What remains is The third and fourth terms in the right-hand side need special attention. Since ρ = ρ(g) is a C 1 monotonically increasing function, we have the elementary inequality Using this inequality and the expression for G (see (76)) in these terms gives Applying Korn's inequality, see Theorem 1.33 from Roubiček (2005), and inserting inequality (98) into equality (96) yields where C i , i = 1, . . . , 4 are positive constants. This proves the coercivity.
Having established pseudo-monotonicity and coercivity of the operator A, we are in position to apply Brézis' theorem. This concludes the assertion of the proposition.

Theorem 1 Problem (PD) admits at least one solution (u, g) ∈ V .
Proof For each ε > 0, let (u ε , g ε ) be a solution of Problem (PD ε ) as obtained in Proposition 2. From the coercivity part of the proof of Proposition 2 and Eq. (84), it follows that where C is independent of ε. Estimate (100) yields weak compactness in H 1 . However, this is not enough to prove that u ε ∇Υ ε converges strongly in L 2 and (a.e.) on Ω as ε → 0. The remedy is to consider the momentum Eq. (83), which gives us improved regularity through the elasticity term. Since p(g ε ) is bounded in H 1 (Ω), uniformly with respect to ε, we conclude that where C does not depend on ε. Using estimates (100)-(101), there is a subsequence (u ε , g ε ), denoted by the same subscript, and a pair (u, g) div u ε → div u strongly in L 2 (Ω) and (a.e) on Ω, g ε → g strongly in L 2 (Ω) and (a.e) on Ω, as ε → 0. The convergence properties allow to pass to the limit in system (82)-(83). Hence, the pair (u, g) satisfies the equations of Problem (PD), which proves the theorem.
To complete the study of the incremental problem, we need to estimate the behavior of solutions after at least O(1/τ ) times steps. Here, we use the discrete version of Lyapunov functional (66).
In Problem (PD), where the discrete time step τ enters as parameter, one find after one step (u 1 , g 1 ) from the initial values (div u, ρ)| t=0 = (0, ρ 0 ). The idea is to repeat this procedure for an arbitrary number of steps. If M ∈ N, M ≤ N = T /τ , then (u M , g M ) denotes the time discretized approximation of the original quasi-static equation, at t = t M = Mτ .
The corresponding Lyapunov functional at t = t M reads It satisfies Here, Proof At time t = t j , with j = 1, . . . , N , the equations in Problem (PD) read Note that in Eq. (109) we have used explicitly the form of D τ from (81). Next, we take ξ = (u j − u j−1 )/τ in (108) and ψ = g j in (109). The resulting two equalities are added and summed up with respect to j up from j = 1 to j = M. Using the observations (i) cross terms containing pressure cancel; where (97) is used; one finds inequality (107). The reduced expression for J 0 results from u| t=0 = 0.
Having established existence for the discrete Problem (PD) in Theorem 1 and a Lyapunov estimate in Theorem 2, we are now in a position to obtain estimates that are uniform in the time step τ .

Proposition 3 There exists a constant C > 0 such that
and for all M and τ such that 1 ≤ M ≤ N = T /τ , with τ sufficiently small.
Proof Combining expression (106) for J M and inequality (107) yields for any 1 ≤ M ≤ N By the assumptions on Q and F, the last two terms are uniformly bounded with respect to τ and M. We estimate the left-hand side from below by applying Korn's inequality to the first term and the quadratic growth of G to the second term. Then for δ and τ sufficiently small, we obtain for the combination where C 1 and C 2 do not dependent on τ and M. Next, we apply the discrete Gronwall inequality 1 , see footnote, to find The second estimate follows directly from Theorem 2.
However, to pass to the limit τ → 0 in the nonlinearities, one needs more information on the behavior of the ratios { div (u j − u j−1 )/τ } and {(g j − g j−1 )/τ }. In fact, we must establish relative compactness of the sequences {div u j } and {g j }. We start with a local H 1 -estimate for E j = div u j .
Lemma 1 Let ϕ ∈ C ∞ 0 (Ω) and τ > 0 sufficiently small. Then, there exists a constant C = C(ϕ) such that Proof Let 1 Discrete version of Gronwall's lemma: Let {U n } and {w n } be nonnegative sequences satisfying U n ≤ A + M−1 j=0 U j w j . Then for all n, U n ≤ A exp{ M−1 j=0 w j }.

Inequality (110) implies
Combined with (111) this gives for L j As in the counterexample for negative porosity, we take the divergence of the timediscrete momentum equation. This yields In general, however, there are no boundary conditions for L j available. Here, we must rely on local estimates to obtain (112). Let us first write the equation for ϕ L j ∈ H 2 (Ω) ∩ H 1 0 (Ω): Its weak form reads With C = C(ϕ) denoting a generic constant depending on ϕ, we have for 1 ≤ j ≤ M ≤ N . Combining this inequality with (114) gives Next, we multiply expression (113) by τ ϕ and write it as Taking the H 1 -inner product of this expression with ϕE j gives Using the identity when summing up (119) gives Combining this inequality with (111) and (118), results in the estimate of the lemma.
We conclude this section with an estimate for (ρ(g j )−ρ(g j−1 ))/τ . However, since in Eqs. (67) or (69) the (discrete) time derivative is multiplied by n(E), we look for an estimate for With the results of Proposition 3 and Lemma 1, the space-time compactness of N will imply the same property of g. We summarize our findings in the next proposition Proposition 4 For given τ > 0 and j = 1, . . . , N , let (u τ (t j ), g τ (t j )) ∈ V denote a solution of Problem (PD). Then, we have where and where ϕ ∈ C ∞ 0 (Ω).

Proof
We only need to prove estimate (124). Rewriting Eq. (79), we have Recalling that for m ≤ 3, and the full estimate reads The local estimate for the space derivatives is given by This results in estimate (124).
In this section, we investigate the limit τ 0. Here a crucial role is played by the parameter λ * , which is needed to control the behavior in time of E = div u.
We are now in a position to prove the main existence result for a weak solution of the time-continuous case.

Proof
In the proof, we use approximations (127) and (128), and their convergence properties.
Multiplying Eq. (152) by α ∈ C ∞ 0 (0, T ) and integrating the result over (τ, T ), yields Next, we send τ 0 along the appropriate subsequence to have convergence of the terms containing u, E and F. What remains is the pressure term. We recall that p(g) is the composite function ( p • ρ)(g), where p(ρ) is given by (36)and ρ(g) is defined through (53) and (55). Since g τ → g strongly in L 2 ((0, T ) × ω), see (147), we have similarly This concludes the first part of the proof.

Discussion and Conclusion
In this paper, we study a model that describes the quasi-static mechanical behavior of a fluid saturated porous medium. In it simplest (linear) form, it is described by Eqs. (1)-(4), where (1) results from the fluid phase mass balance in the case that the fluid is incompressible.
We follow Rutqvist et al. (2001) and Lewis and Schrefler (1998) and propose a fluid mass balance that is based on the mixture theory of Bedford and Drumheller (1978) and Bedford and Drumheller (1983). This yields Eq. (6) and the resulting nonlinear system is given by (2)-(4) and (6). Note that the time derivative of the fluid density ρ appears in (6), since the fluid is assumed weakly compressible, see expression (10). Models where the fluid density is constant (see Bociu et al. 2016;Cao et al. 2013) do not contain this source term. Moreover, the porosity n and the deformation of the medium are related through (8). An expression for this relation is derived from the solid phase mass balance. It is given by (19) or, when the deformation is small, by approximation (21).
It is shown by means of a counterexample that the porosity may admit non-physical, i.e., negative, values. This is made precise in Proposition 1. To obtain a well-posed mathematical problem, the porosity is modified according to cutoff (35). This cutoff is chosen such that it reduces to the correct expression in the physical range. Outside this range, it remains positive. Likewise, a cutoff for the density is introduced through expressions (60) and (64).
The momentum balance Eq.
(2)-(4) is modified as well. Following Murad and Cushman (1996), we add the term to the expression for the total stress. This results in expression (5). Murad and Cushman give a thermodynamically based derivation of the equation in which (166) appears as the difference between the fluid and solid pressures. Having λ * > 0, (166) acts as a time regularization of the volumetric stress for our quasi-static problem. An important role in the analysis of the equations is played by the free energy of the system. This free energy acts as a Lyapunov functional. It is given by (66), which generalizes Biot's original expression developed for the linear case (Biot 1962). In the case that the deformation and fluid density are in the physical range, the free energy simplifies to, see also (56), Ge(u) : e(u) − F · u + n(div u) β 0 n 0 (1 − n 0 )ρ 0 (1 − n 0 )ρ − ρ n 0 0 ρ 1−n 0 + n 0 ρ 0 dx. (167) We introduce a weak formulation and prove existence of a solution in a number of steps. Discretizing in time, we first consider the incremental equations. Using Brézis' fundamental theorem for pseudo-monotone operators, see for instance Lions (1969) and Roubiček (2005), we obtain existence for the corresponding incremental problem. The result holds for any λ * ≥ 0. Moreover, using the free energy, estimates that are global in time are derived. These (stability) estimates are crucial when considering the time-continuous, quasi-static, formulation for which we prove existence at the expense of having λ * > 0. The free energy implies global stability of the solution.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.