<figure><img src="image_1.png" /></figure>

HAL
open science

# Analysis of a linearized poromechanics model for incompressible and nearly incompressible materials

Mathieu Barré, Céline Grandmont, Philippe Moireau

<figure><img src="image_2.png" /></figure>

## To cite this version:

Mathieu Barré, Céline Grandmont, Philippe Moireau. Analysis of a linearized poromechanics model for incompressible and nearly incompressible materials. Evolution Equations and Control Theory, In press, 12 (3), 10.3934/eect.2022053. hal-03501526v2

HAL Id: hal-03501526

https://inria.hal.science/hal-03501526v2

Submitted on 13 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Analysis of a linearized poromechanics model for incompressible and nearly incompressible materials

Mathieu Barré¹,², Céline Grandmont³,⁴,⁵, Philippe Moireau¹,²

¹ Inria, 1 Rue Honoré d'Estienne d'Orves, 91120 Palaiseau, France

² LMS, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91120 Palaiseau, France

³ Inria, 2 Rue Simone Iff, 75012 Paris, France

⁴ LJLL, Sorbonne Université, CNRS, 4 Place Jussieu, 75005 Paris, France

⁵ Département de Mathématique, Université Libre de Bruxelles, CP 214, Boulevard du Triomphe, 1050 Bruxelles, Belgium

January 11, 2023

In this work, we thoroughly analyze the linearized version of a poromechanics model developed to simulate soft tissues perfusion. This is a fully unsteady model in which the fluid and solid equations are strongly coupled through the interstitial pressure. As such, it generalizes Darcy, Brinkman and Biot equations of poroelasticity. The mathematical and numerical analysis of this model was initially performed for a compressible porous material. Here, we focus on the nearly incompressible case with a semigroup approach, which also allows us to prove the existence of weak solutions. We show the existence and uniqueness of strong and weak solutions in the incompressible limit case, for which a divergence constraint on the mixture velocity appears. Due to the special form of the coupling, the underlying problem is not coercive. Nevertheless, by using the notion of *T*-coercivity, we obtain stability estimates and well-posedness results. Our study also provides guidelines to propose a stable and robust approximation of the problem with mixed finite elements. In particular, we recover an inf-sup condition that is independent of the porosity. Finally, we numerically investigate the elliptic regularity of the associated steady-state problem and illustrate the sensitivity of the solution with respect to the different model parameters.

## Introduction

Poromechanical models aim at describing the mechanical behavior of saturated porous media with the interaction of a fluid flow within a deformable porous structure through the definition of a multi-phase continuum framework [51, 55]. The initial introduction of such models concerns geophysics [16, 102], but these models have been recently used for biomechanical applications, in particular to represent perfused living tissues. If the heart perfusion remains a leading example of application [109, 66, 79, 45, 75, 44], poroelastic models have also been considered to simulate lipid and drug transport in blood vessel walls [67, 39, 9, 53, 40], water transport and drug delivery in the brain [13, 99, 103, 104, 48, 64, 70], ocular diseases such as glaucoma [42, 87], fibrosis diagnosis in the lungs [15, 61], or also tissue growth [4, 89, 56].

In these biomedical applications, physical phenomena such as fluid inertia and solid quasi-incompressibility, generally neglected in civil engineering, may play an important role. Therefore, the original poroelasticity model derived by Terzaghi [101] and Biot [16] must be revised to include inertial effects. Note that in the many applications of poroelasticity, unsteady behavior for the fluid and the solid is typically included when studying wave propagation in porous media, see [93] and references therein. It is also an important topic for the simulation of fluid-porous structure interaction (FPSI) occurring in living tissues [95, 9, 34, 33, 35, 43, 5, 3, 1, 40, 20]. In addition to inertial effects, perfused organs such as the heart or the
lungs are subject to finite strains, so their modeling must also account for these non-linear effects and, in particular, consider porosity – which represents the fraction of fluid in the porous material – as a primary variable. Such modeling extensions were proposed within the framework of Biot theory were the solid skeleton plays a special role [30, 73, 62], or in the context of mixture theory treating equivalently all components of the mixture [31, 108, 84]. All these models suppose – explicitly or implicitly – that the frictional effects within the fluid can be neglected due to its viscosity, and rarely take into consideration the influence of solid viscosity. Recently, authors in [46] have revisited the framework of Biot theory at finite strain to derive general formulations adapted to soft tissues perfusion, including inertial and viscous effects both for the fluid and the solid.

Their formulation is compatible with thermodynamical principles. In particular, the solution of the linearized version of the fully coupled model proposed in [46] satisfies energy estimates, opening the way to prove well-posedness. In [38, 10] the case where the structure is compressible is considered for a linearized system closed to the one considered here. Still, the general resulting formulation can exhibit – when solid viscosity is neglected – a hyperbolic-parabolic coupling between the structure and the fluid, with – when the skeleton is incompressible – an additional incompressibility constraint involving a mixture velocity, and therefore leads to challenging questions of analysis.

From a mathematical point of view, there is a large literature related to the existence and uniqueness of solutions for linear Biot’s consolidation models, namely systems of the form

$$ \rho \partial_{tt} u_s - (\lambda + \mu) \nabla (\operatorname{div} u_s) - \mu \Delta u_s + \alpha \nabla p = f, \quad (1a) $$

$$ \partial_t(c_0 p + \alpha \operatorname{div} u_s) - \operatorname{div}(k_f \nabla p) = g, \quad (1b) $$

where the two unknowns are the displacement of the structure $u_s$ and the interstitial pressure $p$, which corresponds to the fluid pressure in the pores. For the unsteady system ($\rho > 0$), the existence of strong solutions was first derived in [52] using Laplace transform and then completed by [60], and the existence of weak solutions was obtained in [12] with a Galerkin method and a regularization technique. The quasi-static case ($\rho = 0$) was first studied in [6] where it was recovered using homogenization techniques, leading to the existence of strong solutions. Existence of weak solutions was shown in [111] using a Galerkin approach, which was recently refined to get a more regular solution [80]. In [94], existence of strong but also weak solutions is established by means of a semigroup approach. This article also handles secondary consolidation phenomena occuring in clays [78], modeled by the presence of an extra term $-\nabla(\lambda^* \partial_t(\operatorname{div} u_s))$ in (1a). Non-linear extensions of (1) were also analyzed [97, 96, 41, 22, 29, 24, 23]. Yet, in the previous models, fluid inertial effects are neglected and, apart from [94, 23], little attention is paid to the incompressible case $c_0 = 0$. Moreover, fluid inertial effects are included in porous wave propagation models [18], whose well-posedness was studied in [92] and [59] using respectively Galerkin and semigroup approaches. However, the fluid viscosity is still not considered and the existence of solutions is carried out only for a compressible fluid, while the fluids present in biomedical applications (blood, lymph, cerebrospinal fluid) are mostly incompressible. Finally, [38] and [10] take into account inertial and viscous fluid effects as their formulation are derived from the linearization of [46] and show respectively the existence of a strong solution when solid viscosity is included, and the existence of a weak solution in absence of solid viscosity, both for a compressible solid. The existence result for incompressible or nearly-incompressible materials was not covered by their results.

In the present work, we study the well-posedness for a linearized system, obtained by linearizing the fully coupled system introduced in [46], by unifying semigroup and variational approaches. The considered model takes into account both fluid and structure inertia, the fluid viscosity, possible damping in the structure, a friction force between both phases, and the interstitial pressure. The elastic or viscoelastic skeleton can be compressible or incompressible, so that we consider four different cases. Our results include the compressible fully viscous case originally studied in [38] and generalize, by relaxing the condition on the fluid mass source term, the results on the compressible elastic case obtained in [10]. Note moreover that the linearized system here considered differs slightly from the one studied in [38, 10], since it incorporates the Biot-Willis coefficient that models pressure-deformation coupling, hence relating the proposed model to the forementioned Biot-type systems. In addition to the compressible case, we fully analyze the incompressible limit case, which corresponds to the physiological regime when considering living tissues. Our approach exploits the notion of
T-coercivity [50, 47] to prove, when no damping is added to the structure, the surjectivity of the underlying operator that involves the resolution of a non-coercive problem. Furthermore, we also take advantage of the recent parallel between inf-sup conditions and T-coercivity [11] to prove the fundamental inf-sup condition associated with the mixture velocity constraint that we have to deal with in the incompressible case. It appears that the inf-sup condition is ultimately independent of the porosity. This result, already conjectured in [38] and partially justified in [10], is crucial to be able to use generic finite-element discretization. It would also be essential when considering the discretization of the non-linear model from [46], in which the porosity is an unknown of the system.

The paper is organized as follows. Section 1 presents the poromechanics model under study, its connection with standard Biot models, and general preliminaries such as the formal derivation of energy estimates on the system. Further details concerning the full non-linear model introduced in [46] and its linearization are given in Appendix A. In Section 2, we unify the semigroup and variational approaches used in [38] and [10] by proving the existence and uniqueness of strong and weak solution for a compressible porous material. We highlight the role of solid viscosity on the model by pointing out the differences that appear in the variational formulation when this coefficient vanishes. Section 3 is devoted to the incompressible regime and more specifically to the saddle-point structure of the problem arising in this case, with a particular attention dedicated to the existence and regularity of pressure. Next, in Section 4, we establish a link between the results of Sections 2 and 3 by passing to the incompressible limit obtained when the bulk modulus of the structure skeleton goes to infinity. Finally, these theoretical results are complemented with numerical experiments exploring the regularity of solutions and the domain of the underlying semigroup operator.

# 1 Problem setting

The model, close to the one we consider here, was introduced in [38] and further explored in [10, 27]. This model comes from the linearization of the poromechanical model developed in [46]. For the sake of completeness, we refer the reader to Appendix A for a brief presentation of the non-linear model proposed in [46] and details about the linearization process in which we introduce the Biot-Willis coefficient that was not taken into account in [38, 10, 27]. We will explain in Section 1.1 that the resulting linearized model is a variant of the well-known Biot systems [16, 17, 19]. Its peculiarity compared to Biot-type models is that it incorporates inertial and viscous effects for both the fluid and the solid, and satisfies an energy balance which is formally derived in Section 1.2 and further rigorously justified.

We consider a porous medium in a bounded domain $\Omega \subset \mathbb{R}^d$ ($d=2,3$) with Lipschitz boundary. In each point of the domain $\Omega$, we consider a mixture of fluid and structure and we denote by $\phi$ the porosity. For all $x \in \Omega$, $\phi(x) \in (0,1)$ represents the fraction of fluid whereas $1-\phi(x)$ represents the fraction of elastic medium. The fluid phase is assumed to be an homogeneous, viscous, Newtonian and incompressible fluid. We denote by $v_f$ its velocity, $\rho_f$ its density and $\mu_f$ its viscosity. Since the fluid is incompressible (resp. homogeneous), $\rho_f$ is independent of time (resp. of space). We also assume that the structure is elastic and, to simplify, that its macroscopic behavior law is linear and isotropic and, thus, characterized by two Lamé constants $\lambda$ and $\mu$. They stand for the elastic parameters characterizing the macroscopic behavior of the elastic part of the mixture (i.e. the homogenized behavior of a perforated elastic medium with no fluid). We denote by $u_s$ the structure displacement and by $v_s = \partial_t u_s$ the structure velocity. The density of the structure is denoted by $\rho_s$ and its viscosity by $\eta$. The fluid and the structure are coupled through a friction force that depends linearly on the relative velocity $v_f - v_s$ and reads $\phi^2 k_f^{-1} (v_f - v_s)$, where $k_f$ is the hydraulic conductivity tensor. In addition, they are coupled through the interstitial pressure $p$, that is further linked to the incompressibility of the whole fluid-structure mixture. Finally, $\alpha(x) \in (\phi(x), 1)$ is the Biot-Willis coefficient, which takes into account the pressure-deformation coupling. This coefficient depends on space for a compressible material but tends to 1 in the incompressible limit as the skeleton elastic bulk modulus, denoted $\kappa$, tends to $+\infty$, see Remark A.1.
The fully coupled model then reads

$$
\left\{
\begin{aligned}
& \rho_s(1-\phi) \partial_{tt} u_s - \operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(\partial_t u_s)) \\
& \quad - \phi^2 k_f^{-1} (v_f - \partial_t u_s) + (\alpha - \phi) \nabla p = \rho_s(1-\phi) f, && \text{in } \Omega \times (0, T), \quad (2a) \\
& \rho_f \phi \partial_t v_f - \operatorname{div}(\phi \sigma_f(v_f)) + \phi^2 k_f^{-1} (v_f - \partial_t u_s) - \theta v_f + \phi \nabla p = \rho_f \phi f, && \text{in } \Omega \times (0, T), \quad (2b) \\
& \frac{\alpha - \phi}{\kappa} \partial_t p + \operatorname{div}((\alpha - \phi) \partial_t u_s + \phi v_f) = \frac{\theta}{\rho_f}, && \text{in } \Omega \times (0, T), \quad (2c)
\end{aligned}
\right.
$$

where the structure stress tensor is given by Hooke's law

$$
\sigma_s(u) = \lambda \operatorname{Tr}(\varepsilon(u)) \mathcal{I} + 2\mu \varepsilon(u),
$$

with $\varepsilon(u) = \frac{1}{2}(\nabla u + \nabla u^T)$, the structure additional viscosity is given by

$$
\sigma_s^{\text{vis}}(v) = 2\eta\varepsilon(v),
$$

and the fluid stress tensor reads

$$
\sigma_f(v) = \lambda_f \operatorname{Tr}(\varepsilon(v)) \mathcal{I} + 2\mu_f \varepsilon(v).
$$

In the above system, the data are the applied exterior force *f* and the additional fluid mass input *θ*. The coupled system (2) describes the mixture of an elastic, possibly viscous medium and an incompressible Newtonian flow. The first equation (2a) represents the momentum conservation law of the elastic phase including inertial effects, macroscopic elastic behavior, possible viscous damping, friction force between the fluid and the structure, and the gradient of the interstitial pressure. The second equation (2b) stands for the momentum conservation law of the fluid phase including inertial effects, macroscopic viscous effects, friction force and the gradient of the interstitial pressure. The third equation (2c) traduces the total mass conservation dynamic, it involves the parameter *κ*, that represents the bulk modulus of the elastic medium constituting the porous matrix, and the Biot-Willis parameter *α*. When *κ* < +∞ it corresponds to a compressible skeleton, whereas when *κ* = +∞ (that implies *α* = 1) we have an incompressible elastic skeleton. Since the fluid is assumed to be incompressible we deal in the limit case *κ* = +∞ with an incompressible porous medium. This latter case is crucial when considering living tissues since they are nearly incompressible. Note that when *κ* = +∞ there is no dynamic of the pressure since the term ∂<sub>t</sub>p in (2c) vanishes, but the pressure is the Lagrange multiplier associated with the mixture constraint ρ<sub>f</sub> div ((1 − φ) ∂<sub>t</sub>u<sub>s</sub> + φv<sub>f</sub>) = θ involving the mixture velocity v<sub>m</sub> = (1 − φ) ∂<sub>t</sub>u<sub>s</sub> + φv<sub>f</sub>. Further details on the derivation of this linearized coupled system are gathered in the Appendix A.

For a presentation of typical boundary conditions for such systems, we refer to [37, 38, 90, 40] and
their analysis will imply further development. In the present work, we limit our analysis to the case of
homogeneous Dirichlet boundary conditions for the structure and for the fluid:

$$
\left\{
\begin{array}{ll}
u_s = 0, & \text{on } \partial\Omega \times (0, T), \\
v_f = 0, & \text{on } \partial\Omega \times (0, T).
\end{array}
\right.
\qquad (3a)
$$

This coupled problem has to be completed with initial data:

$$
\left\{
\begin{array}{ll}
u_s(0) = u_{s0}, & \text{in } \Omega, \\
\partial_t u_s(0) = v_{s0}, & \text{in } \Omega, \\
v_f(0) = v_{f0}, & \text{in } \Omega,
\end{array}
\right.
\qquad (4a)
$$

and in the case $\kappa < +\infty$

$$
p(0) = p_0, \quad \text{in } \Omega. \tag{5}
$$

Before detailing the well-posedness analysis of the considered coupled system, let us first emphasize its
links to other systems modeling porous media.
## 1.1 Related poromechanics models

As shown in [85], Darcy, Brinkman and Biot equations can be derived within the framework of mixture theory under specific assumptions. The system (2), which arises from Biot theory, can be seen as a combination of these models. Indeed, (2) is close to the fully dynamic Biot system introduced in [18] for the study of acoustic waves in saturated porous media, but also includes a viscous fluid term as in Brinkman equation.

More precisely, denoting by $u_f$ the displacement of fluid particles within the porous medium and by $w = \phi(u_f - u_s)$ the relative displacement of the fluid phase with respect to the solid one, the model from [18] reads

$$ \rho \partial_{tt} u_s + \rho_f \partial_{tt} w - \operatorname{div} (\sigma_s(u_s)) + \alpha \nabla p = g, \quad (6a) $$

$$ \rho_f \partial_{tt} u_s + a \rho_f \partial_{tt} \left( \frac{w}{\phi} \right) + k_f^{-1} q + \nabla p = h, \quad (6b) $$

$$ c_0 p + \alpha \operatorname{div} u_s + \operatorname{div} w = k, \quad (6c) $$

where $\rho = \rho_s(1 - \phi) + \rho_f\phi$ corresponds to the density of the mixture, $a \ge 1$ is a coefficient describing tortuosity effects, and

$$ c_0 = \frac{\phi}{\kappa_f} + \frac{\alpha - \phi}{\kappa} $$

is the storage coefficient, with $\kappa_f$ the fluid bulk modulus.

In our case, the fluid is assumed to be incompressible and thus $\kappa_f = +\infty$, so that $c_0 = \frac{\alpha - \phi}{\kappa}$. To link (2) and (6), let us assume that we have no additional fluid mass input, namely $\theta = 0$, and that we can neglect viscous effects, which amounts to take $\eta = \mu_f = \lambda_f = 0$. Introducing the new unknown

$$ q = \phi(v_f - \partial_t u_s) = \partial_t w, $$

which corresponds to the filtration velocity, (2) becomes

$$ \rho_s(1 - \phi) \partial_{tt} u_s - \operatorname{div} (\sigma_s(u_s)) - \phi k_f^{-1} q + (\alpha - \phi) \nabla p = \rho_s(1 - \phi) f, \quad (7a) $$

$$ \rho_f \phi \partial_{tt} u_s + \rho_f \partial_t q + \phi k_f^{-1} q + \phi \nabla p = \rho_f \phi f, \quad (7b) $$

$$ c_0 \partial_t p + \operatorname{div} (\alpha \partial_t u_s + q) = \rho_f^{-1} \theta. \quad (7c) $$

Replacing (7a) by (7a) + (7b) and dividing (7b) by $\phi$, we get

$$ \rho \partial_{tt} u_s + \rho_f \partial_t q - \operatorname{div} (\sigma_s(u_s)) + \alpha \nabla p = \rho f, \quad (8a) $$

$$ \rho_f \partial_{tt} u_s + \rho_f \partial_t \left( \frac{q}{\phi} \right) + k_f^{-1} q + \nabla p = \rho_f f, \quad (8b) $$

$$ \partial_t (c_0 p + \alpha \operatorname{div} u_s) + \operatorname{div} q = \rho_f^{-1} \theta, \quad (8c) $$

which, provided that $a = 1$, corresponds exactly to (6) since $q = \partial_t w$ and (8c) $= \partial_t (6c)$. Note that if $c_0 > 0$, equation (6c) can be used to eliminate the pressure unknown as done in [112, 92, 59], but it is no longer the case if we consider (8c). The assumption $a = 1$ indicates that (2) does not take into account tortuosity effects since they are not compatible with the first principle of continuum mechanics introduced in [46], see [62, Section 5.3.4] for a discussion on the thermodynamical compatibility of these effects and [73] for a fully unsteady poromechanical model in which they are included.

If the fluid and solid inertial effects are also neglected, (8a) reduces to (1a) and (8b) implies that $q = -k_f \nabla p + \rho_f k_f f$. Substituting this result in (8c), we recover the quasi-static Biot's consolidation model, namely (1) with $\rho = 0$. Therefore, the model studied in this paper is connected to Darcy, Brinkman and Biot equations, but the presence of inertial and viscous terms both for the fluid and the solid requires a separate study. In particular, because of these extra terms, the functional setting adapted to the problem differs from the one developed for Biot models. This functional setting is guided by the energy balance presented below.
**Remark 1.1.** Darcy, Brinkman and Biot models have been justified a posteriori using homogeneization techniques, see for instance [7, 65, 76, 86] and references therein. The justification of (2) by homogeneization is an open problem.

## 1.2 Energy estimate

Existence of solutions of such a coupled system, in the compressible case $\kappa < +\infty$, has been partially obtained in [38, 10]. More precisely, the case $\kappa < +\infty$, $\eta > 0$ and $\theta = 0$ has already been studied in [38], where existence of strong solutions thanks to the semigroup formalism has been derived. The case $\kappa < +\infty$, $\eta = 0$ is treated in [10], where existence of variational solutions is obtained under a smallness assumption on $\theta$. Here, we consider all the different cases $\kappa \le +\infty$, $\eta \ge 0$, and any given $\theta$ sufficiently smooth. We prove existence of unique strong and mild solutions – in a sense to be made precise later – using semigroup theory, from which we deduce existence of a unique variational solution. We eventually show that one can pass to the limit in the weak formulation as $\kappa$ goes to infinity.

Before going through the proves, let us first derive formally some energy bounds satisfied by any smooth enough solutions of the coupled problem. We first derive them in the case $\kappa < +\infty$ and then in the limit case $\kappa = +\infty$. Let us multiply (2a) by the structure velocity $\partial_t u_s$, integrate over $\Omega$ and integrate by parts in space. No boundary terms appear thanks to the homogeneous Dirichlet boundary conditions (3a) and we obtain

$$ \begin{aligned} & \frac{\rho_s}{2} \frac{d}{dt} \int_{\Omega} (1-\phi) |\partial_t u_s|^2 dx + \frac{1}{2} \frac{d}{dt} \int_{\Omega} \sigma_s(u_s) : \varepsilon(u_s) dx + 2\eta \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(\partial_t u_s) dx \\ & \qquad - \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot \partial_t u_s dx - \int_{\Omega} p \operatorname{div}((\alpha-\phi) \partial_t u_s) dx = \int_{\Omega} \rho_s (1-\phi) f \cdot \partial_t u_s dx. \end{aligned} $$

Let us also multiply (2b) by the fluid velocity $v_f$, integrate over $\Omega$ and integrate by parts in space. No boundary terms appear thanks to the homogeneous Dirichlet boundary conditions (3b) and we get

$$ \begin{aligned} & \frac{\rho_f}{2} \frac{d}{dt} \int_{\Omega} \phi |v_f|^2 dx + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) dx \\ & \qquad + \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot v_f dx - \int_{\Omega} \theta |v_f|^2 dx - \int_{\Omega} p \operatorname{div}(\phi v_f) dx = \int_{\Omega} \rho_f \phi f \cdot v_f dx. \end{aligned} $$

The last equation (2c) is multiplied by $p$ and integrated over $\Omega$, which leads to

$$ \frac{1}{2} \frac{d}{dt} \int_{\Omega} \frac{\alpha - \phi}{\kappa} |p|^2 dx + \int_{\Omega} \operatorname{div}((\alpha - \phi) \partial_t u_s + \phi v_f) p dx = \int_{\Omega} \frac{\theta}{\rho_f} p dx. \quad (9) $$

Adding these three contributions, we see that the terms involving the divergence of the mixture velocity $v_{m,\alpha} = (\alpha - \phi)\partial_t u_s + \phi v_f$ cancel, and we have the following energy equality

$$ \begin{aligned} & \frac{\rho_s}{2} \frac{d}{dt} \int_{\Omega} (1-\phi) |\partial_t u_s|^2 dx + \frac{1}{2} \frac{d}{dt} \int_{\Omega} \sigma_s(u_s) : \varepsilon(u_s) dx + 2\eta \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(\partial_t u_s) dx + \frac{\rho_f}{2} \frac{d}{dt} \int_{\Omega} \phi |v_f|^2 dx \\ & \qquad + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) dx + \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (v_f - \partial_t u_s) dx + \frac{1}{2} \frac{d}{dt} \int_{\Omega} \frac{\alpha - \phi}{\kappa} |p|^2 dx \\ & = \int_{\Omega} \rho_s (1-\phi) f \cdot \partial_t u_s dx + \int_{\Omega} \rho_f \phi f \cdot v_f dx + \int_{\Omega} \frac{\theta}{\rho_f} p dx + \int_{\Omega} \theta |v_f|^2 dx. \quad (10) \end{aligned} $$

Consequently, in order to obtain an energy estimate, we impose the following assumptions on the data:

(h1) The constants $\rho_s, \rho_f, \mu_f, \lambda, \mu$ are assumed to be strictly positive, whereas $\eta \ge 0$;

(h2) The porosity $\phi \in H^{d/2+r}(\Omega)$ with $r > 0$, and is such that there exists $(\phi_{\min}, \phi_{\max})$ satisfying

$$ 0 < \phi_{\min} \le \phi(x) \le \phi_{\max} < 1, \quad \forall x \in \Omega; $$
(h3) The friction tensor $k_f$ is invertible and there exists $k_0 > 0$ such that

$$k_f^{-1}v \cdot v \geq k_0|v|^2, \quad \forall v \in \mathbb{R}^d;$$

(h4) $f \in L^2((0,T) \times \Omega)$;

(h5) $\theta \in C^0([0,T] \times \Omega)$;

(h6) The (non-homogeneous) Biot-Willis coefficient $\alpha \in H^{d/2+r}(\Omega)$ with $r > 0$, and is such that there exists $((\alpha - \phi)_{\min}, (\alpha - \phi)_{\max})$ satisfying

$$0 < (\alpha - \phi)_{\min} \leq \alpha(x) - \phi(x) \leq (\alpha - \phi)_{\max} < 1, \quad \forall x \in \Omega;$$

**Remark 1.2.** The hypotheses (h2) and (h6) imply that the porosity $\phi$ and the Biot-Willis coefficient $\alpha$ belong to a multiplier space of $H^1(\Omega)$. These assumptions are needed to define the term $\text{div}((\alpha - \phi)\partial_t u_s + \phi v_f)$ in (9). Indeed, if $\alpha, \phi \in H^{d/2+r}(\Omega)$ with $r > 0$, then for any $(w_s, w_f) \in [H_0^1(\Omega)]^d$ we have $\text{div}((\alpha - \phi)w_s + \phi w_f) \in L^2(\Omega)$.

Under these assumptions, using Young inequality to bound the right-hand side of (10) by

$$\frac{1}{2} \int_{\Omega} \rho_s(1-\phi) |f|^2 dx + \frac{1}{2} \int_{\Omega} \rho_s(1-\phi) |\partial_t u_s|^2 dx + \frac{1}{2} \int_{\Omega} \rho_f \phi |f|^2 dx + \frac{1}{2} \int_{\Omega} \rho_f \phi |v_f|^2 dx \\ + \frac{1}{2} \int_{\Omega} \frac{\kappa}{\rho_f^2(\alpha-\phi)_{\min}} |\theta|^2 dx + \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |p|^2 dx + \frac{2\|\theta\|_{C^0([0,T]\times\Omega)}}{\rho_f \phi_{\min}} \cdot \frac{1}{2} \int_{\Omega} \rho_f \phi |v_f|^2 dx,$$

integrating in time from 0 to $t$ and applying Grönwall Lemma, we obtain the following energy bound

$$\begin{aligned} & \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |\partial_t u_s(t)|^2 dx + \frac{1}{2} \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(u_s(t)) dx + 2\eta \int_0^t \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(\partial_t u_s) dx ds + \frac{\rho_f}{2} \int_{\Omega} \phi |v_f(t)|^2 dx \\ & + \int_0^t \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) dx ds + \int_0^t \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (v_f - \partial_t u_s) dx ds + \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |p(t)|^2 dx \\ & \leq \exp \left( \max \left( 1, \frac{2\|\theta\|_{C^0([0,T]\times\Omega)}}{\rho_f \phi_{\min}} \right) t \right) \left( \left( \frac{\rho_s}{2}(1-\phi_{\min}) + \frac{\rho_f}{2} \phi_{\max} \right) \int_0^t \int_{\Omega} |f|^2 dx ds \right. \\ & \qquad \left. + \frac{\kappa}{2\rho_f^2(\alpha-\phi)_{\min}} \int_0^t \int_{\Omega} |\theta|^2 dx ds + \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |v_{s0}|^2 dx \right. \\ & \qquad \left. + \frac{1}{2} \int_{\Omega} \sigma_s(u_{s0}) : \varepsilon(u_{s0}) dx + \frac{\rho_f}{2} \int_{\Omega} \phi |v_{f0}|^2 dx + \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |p_0|^2 dx \right). \end{aligned} \quad (11)$$

Note that the friction contribution induces dissipation in the system since

$$\int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (v_f - \partial_t u_s) dx \geq 0,$$

in virtue of (h3).

Moreover, Korn inequality [57, 49] implies that the fluid and structure dissipative terms are coercive in $[H_0^1(\Omega)]^d$. Namely, there exists $C_d > 0$ such that

$$\forall v \in [H_0^1(\Omega)]^d, \quad \int_{\Omega} \varepsilon(v) : \varepsilon(v) dx \geq C_d \|v\|^2_{[H_0^1(\Omega)]^d}, \quad (12)$$
which implies that the bilinear elastic form is coercive in $[H_0^1(\Omega)]^d$ and verifies

$$ \forall v \in [H_0^1(\Omega)]^d, \quad \int_{\Omega} \sigma_s(v) : \varepsilon(v) \, dx \geq 2\mu C_d \|v\|^2_{[H_0^1(\Omega)]^d}. \quad (13) $$

For the fluid part, thanks to assumption (h2), one also has

$$ \forall v \in [H_0^1(\Omega)]^d, \quad \int_{\Omega} \phi \sigma_f(v) : \varepsilon(v) \, dx \geq 2\mu_f \phi_{\min} C_d \|v\|^2_{[H_0^1(\Omega)]^d}. \quad (14) $$

Consequently, assuming that $(u_{s0}, v_{s0}, v_{f0}, p_0) \in [H_0^1(\Omega)]^d \times [L^2(\Omega)]^d \times [L^2(\Omega)]^d \times L^2(\Omega)$, it follows that $u_s \in L^\infty(0, T; [H_0^1(\Omega)]^d)$, $\partial_t u_s \in L^\infty(0, T; [L^2(\Omega)]^d)$, $v_f \in L^\infty(0, T; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$, $p \in L^\infty(0, T; L^2(\Omega))$, and that $\partial_t u_s \in L^2(0, T; [H_0^1(\Omega)]^d)$ if $\eta > 0$.

Note that the energy bound (11) depends on the bulk modulus $\kappa$. Nonetheless, if $\theta$ is regular enough, we can recover an energy estimate independent of $\kappa$ by coming back to the case where the right-hand side of (2c) is equal to zero, as we are now going to perform it in the incompressible case.

Let us now focus on the case $\kappa = +\infty$ for which $\alpha = 1$. The equation (2c) reduces to

$$ \operatorname{div}((1 - \phi)v_s + \phi v_f) = \frac{\theta}{\rho_f}, \quad \text{in } \Omega. \quad (15) $$

Without loss of generality we can assume that the right-hand side of (15) is equal to zero. Indeed, provided that $\int_{\Omega} \theta \, dx = 0$, there exists $v_\theta$ such that $\operatorname{div} v_\theta = \frac{\theta}{\rho_f}$. Considering the system satisfied by $v_s - v_\theta$ and $v_f - v_\theta$, namely defining the new displacement

$$ u_{s0} + \int_{0}^{t} (v_s - v_{\theta}) \, ds = u_s - \int_{0}^{t} v_{\theta} \, ds, $$

we end up with a system for which the constraint reads $\operatorname{div}((1 - \phi)v_s + \phi v_f) = 0$. To obtain the energy estimates, we proceed as for the case $\kappa < +\infty$ by multiplying (2a) by the structure velocity $\partial_t u_s$, and (2b) by the fluid velocity $v_f$. After integration over the domain and integration by parts, adding these two contributions and taking into account the mixture incompressibility constraint $\operatorname{div}((1 - \phi)v_s + \phi v_f) = 0$ yields

$$ \begin{aligned} & \frac{\rho_s}{2} \frac{d}{dt} \int_{\Omega} (1-\phi) |\partial_t u_s|^2 \, dx + \frac{1}{2} \frac{d}{dt} \int_{\Omega} \sigma_s(u_s) : \varepsilon(u_s) \, dx + 2\eta \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(\partial_t u_s) \, dx \\ & \qquad + \frac{\rho_f}{2} \frac{d}{dt} \int_{\Omega} \phi |v_f|^2 \, dx + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) \, dx + \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (v_f - \partial_t u_s) \, dx \\ & = \int_{\Omega} \rho_s (1-\phi) f \cdot \partial_t u_s \, dx + \int_{\Omega} \rho_f \phi f \cdot v_f \, dx + \int_{\Omega} \theta |v_f|^2 \, dx. \quad (16) \end{aligned} $$

Grönwall Lemma then implies

$$ \begin{aligned} & \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |\partial_t u_s(t)|^2 \, dx + \frac{1}{2} \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(u_s(t)) \, dx + 2\eta \int_{0}^{t} \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(\partial_t u_s) \, dx \, ds \\ & \qquad + \frac{\rho_f}{2} \int_{\Omega} \phi |v_f(t)|^2 \, dx + \int_{0}^{t} \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) \, dx \, ds + \int_{0}^{t} \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (v_f - \partial_t u_s) \, dx \, ds \\ & \leq \exp\left(\max\left(1, \frac{2\|\theta\|_{C^0([0,T]\times\Omega)}}{\rho_f \phi_{\min}}\right)\right) \left(\left(\frac{\rho_s}{2}(1-\phi_{\min}) + \frac{\rho_f}{2}\phi_{\max}\right) \int_{0}^{t} \int_{\Omega} |f|^2 \, dx \, ds \right. \\ & \qquad \left. + \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |v_{s0}|^2 \, dx + \frac{1}{2} \int_{\Omega} \sigma_s(u_{s0}) : \varepsilon(u_{s0}) \, dx + \frac{\rho_f}{2} \int_{\Omega} \phi |v_{f0}|^2 \, dx\right). \quad (17) \end{aligned} $$
Thanks to Korn inequality (12), coercivities (13), (14), assumptions (h1) – (h5) and assuming that $(u_{s0}, v_{s0}, v_{f0}) \in [H_0^1(\Omega)]^d \times [L^2(\Omega)]^d \times [L^2(\Omega)]^d$, we have $u_s \in L^\infty(0, T; [H_0^1(\Omega)]^d)$, $\partial_t u_s \in L^\infty(0, T; [L^2(\Omega)]^d)$, $v_f \in L^\infty(0, T; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$, and if $\eta > 0$, $\partial_t u_s \in L^2(0, T; [H_0^1(\Omega)]^d)$. Here the energy bounds does not give bounds on the pressure, which is the main difference between the cases $\kappa < +\infty$ and $\kappa = +\infty$.

We then propose the following milestones for our analysis. We start by considering the compressible case for which $\kappa < +\infty$. In this case the pressure $p$ has its own dynamic. Then the incompressible case, namely $\kappa = +\infty$, is treated and we have to deal with a divergence-free constraint on the mixture velocity. Each case is split into two cases: the viscous one (namely $\eta > 0$) for which we have a parabolic-parabolic coupling between the solid and fluid equations, and the inviscid one (namely $\eta = 0$) for which we have a hyperbolic-parabolic coupling. For each four cases we prove existence of strong, mild and variational solutions and give the link between the three types of solutions. In particular, existence of strong and mild solutions relies on the study of the first order system of the form $\dot{z} + Az = g$ associated with (2) and the underlying unbounded operator $A$ using semigroup theory. By strong solution, we mean that the solution is regular in time and that the equations are satisfied almost everywhere in the sense that all the components of $\dot{z}$ and $Az$ are defined almost everywhere, whereas mild solutions are solutions satisfying the Duhamel formula. Note that in the case $\eta = 0$ in order to prove that the operator is maximal accretive we need to take care of the non coercivity of the associated bilinear form. This issue is solved thanks to the notion of $T$-coercivity introduced in [50, 47]. Next the variational solutions are obtained by an approximation strategy as the limit of a sequence of strong solutions. Note that, as we will see, the definition of the variational formulations is different when considering $\eta > 0$ or $\eta = 0$. The main difference comes from the fact that, in the latter case, the structure velocity is not in $[H_0^1(\Omega)]^d$ in space but only in $[L^2(\Omega)]^d$. We end up with the study of the incompressible limit, which allows to pass to the limit in the weak formulation for $\kappa < +\infty$, to recover the weak formulation associated with $\kappa = +\infty$. The theoretical results are further completed by numerical illustrations to investigate the regularity of the solutions.

## 2 Existence of solutions for a compressible skeleton $\kappa < +\infty$.

In this section, we study the poromechanical problem for a compressible skeleton, that corresponds to $\kappa < +\infty$. First, we write (2) as a first-order evolution equation and we define the associated unbounded operator. Then, by investigating the properties of this operator, we use a semigroup approach to show existence and uniqueness of strong and mild solutions to the system. The existence of variational solutions is then obtained by an approximation strategy. The cases $\eta > 0$ and $\eta = 0$ are treated separately in order to emphasize the influence of solid viscosity on the model. But let us start with some general notations and definitions valid for both cases.

### 2.1 Semigroup framework

The system (2) can be rewritten as a first order system as follows

$$ 
\left\{
\begin{array}{ll}
\partial_t u_s - v_s = 0, & \text{in } \Omega \times (0, T), \qquad (18a) \\
\rho_s(1-\phi) \partial_t v_s - \text{div} (\sigma_s(u_s)) - \text{div} (\sigma_s^{\text{vis}}(v_s)) & \\
\qquad - \phi^2 k_f^{-1}(v_f - v_s) + (\alpha - \phi) \nabla p = \rho_s(1-\phi) f, & \text{in } \Omega \times (0, T), \qquad (18b) \\
\rho_f \phi \partial_t v_f - \text{div} (\phi \sigma_f(v_f)) & \\
\qquad + \phi^2 k_f^{-1}(v_f - v_s) - \theta v_f + \phi \nabla p = \rho_f \phi f, & \text{in } \Omega \times (0, T), \qquad (18c) \\
\frac{\alpha - \phi}{\kappa} \partial_t p + \text{div} ((\alpha - \phi) v_s + \phi v_f) = \frac{\theta}{\rho_f}, & \text{in } \Omega \times (0, T). \qquad (18d)
\end{array}
\right.
$$

Let $z = (u_s, v_s, v_f, p)$ and $z_0 = (u_{s0}, v_{s0}, v_{f0}, p_0)$ denote respectively the unknown variable and the initial
condition of (18). We formulate (18) as an abstract evolution problem

$$ \begin{cases} \dot{z}(t) + A_{\eta}^{\kappa}z(t) + G(t)z(t) = g(t), & t \in [0, T], \\ z(0) = z_0, & \end{cases} \qquad (19) $$

where $A_{\eta}^{\kappa}$ is an unbounded operator specified with respect to the solid viscosity $\eta$ and the bulk modulus $\kappa$, and $G(t)$ is a bounded perturbation defined below.

Let us first define the energy space

$$ Z = [H_0^1(\Omega)]^d \times [L^2(\Omega)]^d \times [L^2(\Omega)]^d \times L^2(\Omega) $$

associated with (18). Since the functions $\rho_s(1-\phi)$, $\rho_f\phi$ and $\frac{\alpha-\phi}{\kappa}$ are bounded and bounded from below by strictly positive constants, the space $Z$ can be endowed with the scalar product defined by

$$ (z, y)_Z = \int_{\Omega} \sigma_s(u_s) : \varepsilon(d_s) \, dx + \int_{\Omega} \rho_s(1-\phi) v_s \cdot w_s \, dx + \int_{\Omega} \rho_f \phi v_f \cdot w_f \, dx + \int_{\Omega} \frac{\alpha-\phi}{\kappa} pq \, dx, $$

for any $z = (u_s, v_s, v_f, p)$ and $y = (d_s, w_s, w_f, q)$ belonging to $Z$. The associated norm reads

$$ \|z\|_Z^2 = \|u_s\|_s^2 + \int_{\Omega} \rho_s(1-\phi) |v_s|^2 \, dx + \int_{\Omega} \rho_f \phi |v_f|^2 \, dx + \int_{\Omega} \frac{\alpha-\phi}{\kappa} p^2 \, dx, \qquad (20) $$

with

$$ \|u_s\|_s^2 = \int_{\Omega} \sigma_s(u_s) : \varepsilon(u_s) \, dx. \qquad (21) $$

This norm is equivalent to the canonical norm on $Z$ according to Korn inequality (13).

Setting

$$ Y = [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times L^2(\Omega) $$

as an intermediate space, we introduce the bilinear form $a_{\eta}^{\kappa}$ defined for all $z = (u_s, v_s, v_f, p) \in Y$ and $y = (d_s, w_s, w_f, q) \in Y$ by

$$ \begin{aligned} a_{\eta}^{\kappa}(z, y) = & - \int_{\Omega} \sigma_s(v_s) : \varepsilon(d_s) \, dx + \int_{\Omega} \sigma_s(u_s) : \varepsilon(w_s) \, dx + 2\eta \int_{\Omega} \varepsilon(v_s) : \varepsilon(w_s) \, dx \\ & + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f - v_s) \cdot (w_f - w_s) \, dx \\ & + \int_{\Omega} \operatorname{div} ((\alpha - \phi) v_s + \phi v_f) \, q \, dx - \int_{\Omega} p \operatorname{div} ((\alpha - \phi) w_s + \phi w_f) \, dx. \end{aligned} \qquad (22) $$

The bilinear form $a_{\eta}^{\kappa}$ is continuous over $Y \times Y$.

Associated with this bilinear form, we introduce the unbounded operator $(A_{\eta}^{\kappa}, D(A_{\eta}^{\kappa}))$ defined by

$$ (A_{\eta}^{\kappa}z, y)_Z = a_{\eta}^{\kappa}(z, y), \quad \forall z \in D(A_{\eta}^{\kappa}), \forall y \in Y, \qquad (23) $$

in the domain

$$ D(A_{\eta}^{\kappa}) = \{z \in Y : \exists g \in Z, a_{\eta}^{\kappa}(z, y) = (g, y)_Z, \quad \forall y \in Y\}. \qquad (24) $$

Finally, for all $t \in [0, T]$, we define the time-dependent operator

$$ G(t) : z = (u_s, v_s, v_f, p) \in Z \mapsto (0, 0, -\frac{\theta(t)}{\rho_f \phi} v_f, 0). \qquad (25) $$

Taking $g = (0, f, f, \frac{\kappa}{\alpha-\phi} \cdot \frac{\theta}{\rho_f})$, the state-space formulation (19) is equivalent to (18) in a sense that will be specified in Corollary 2.4.
**Remark 2.1.** Note that in the domain of operator the equation stating that the time derivative of the structure displacement is equal to the structure velocity (that comes from the first order rewriting of a second order in time problem) will hold true in $[H_0^1(\Omega)]^d$ in the space variable. This is the reason of the presence of the term $-\int_{\Omega} \sigma_s(v_s) : \varepsilon(d_s) dx$ in (22). Yet, even if the solid velocity is considered in $[H_0^1(\Omega)]^d$ in the latter integral, we will see that when $\eta = 0$ the resulting weak solution does not satisfy (18a) in $[H_0^1(\Omega)]^d$ but only in $[L^2(\Omega)]^d$. The same issue appears when studying the wave equation.

For $z = (u_s, v_s, v_f, p) \in D(A_\eta^\kappa)$, we can write

$$A_{\eta}^{\kappa} z = \begin{pmatrix} -v_s \\ (\rho_s(1-\phi))^{-1}(-\operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) + \phi^2 k_f^{-1}(v_s - v_f) + (\alpha - \phi)\nabla p) \\ (\rho_f\phi)^{-1}(-\operatorname{div}(\phi\sigma_f(v_f)) + \phi^2 k_f^{-1}(v_f - v_s) + \phi\nabla p) \\ \frac{\kappa}{\alpha-\phi}\operatorname{div}((\alpha-\phi)v_s + \phi v_f) \end{pmatrix}, \quad (26)$$

so that the operator $A_\eta^\kappa$ can be expressed in matrix form as

$$A_{\eta}^{\kappa} = N_{0}^{-1} \begin{pmatrix} 0 & -\mathbb{I} & 0 & 0 \\ -\operatorname{div}(\sigma_{s}(\cdot)) & -\operatorname{div}(\sigma_{s}^{\text{vis}}(\cdot)) + \phi^{2} k_{f}^{-1} & -\phi^{2} k_{f}^{-1} & (\alpha - \phi) \nabla \\ 0 & -\phi^{2} k_{f}^{-1} & -\operatorname{div}(\phi \sigma_{f}(\cdot)) + \phi^{2} k_{f}^{-1} & \phi \nabla \\ 0 & \operatorname{div}((\alpha - \phi) \cdot) & \operatorname{div}(\phi \cdot) & 0 \end{pmatrix},$$

where $\mathbb{I}$ denotes the identity operator of the space $[H_0^1(\Omega)]^d$ endowed with the norm (21), and

$$N_0 = \begin{pmatrix} \mathbb{I} & 0 & 0 & 0 \\ 0 & \rho_s(1-\phi) & 0 & 0 \\ 0 & 0 & \rho_f\phi & 0 \\ 0 & 0 & 0 & \frac{\alpha-\phi}{\kappa} \end{pmatrix}.$$

Moreover, from (24) and (26), it follows that

$$D(A_{\eta}^{\kappa}) = \left\{ \begin{array}{ll} u_s, v_s, v_f \in [H_0^1(\Omega)]^d & -\operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) + (\alpha - \phi)\nabla p \in [L^2(\Omega)]^d \\ p \in L^2(\Omega) & -\operatorname{div}(\phi\sigma_f(v_f)) + \phi\nabla p \in [L^2(\Omega)]^d \\ & \operatorname{div}((\alpha - \phi)v_s + \phi v_f) \in L^2(\Omega) \end{array} \right\}. \quad (27)$$

Note that belonging to $D(A_\eta^\kappa)$ does not mean that all the above terms individually belong to $L^2(\Omega)$, but only that their sum does. For instance, $(v_f, p)$ does not necessarily belong to $[H^2(\Omega)]^d \times H^1(\Omega)$, but we know that $-\operatorname{div}(\phi\sigma_f(v_f)) + \phi\nabla p \in [L^2(\Omega)]^d$. Specifying $D(A_\eta^\kappa)$ in terms of classical Sobolev spaces requires to study the regularity of the solution to the static problem $A_\eta^\kappa z = g$ with $g \in Z$. This issue, delicate from a theoretical point of view, will be explored in more details in numerical experiments, see Section 5.

In what follows, we exploit the previous framework to prove that Problem (18) has a unique strong and mild solution for $\kappa < +\infty$. We also recover the existence of variational solutions as the limit of a sequence of strong solutions. If $\eta > 0$, the solid equation (18b) is parabolic, while it becomes hyperbolic when $\eta = 0$. For this reason, we distinguish the cases $\eta > 0$ and $\eta = 0$.

## 2.2 The case $\eta > 0$

Let us start with the parabolic-parabolic coupling configuration. This case was treated in [38] for $\theta = 0$ and $\alpha = 1$. Here, we propose a proof of existence and uniqueness which is valid for a time-dependent $\theta$. Note that considering $\alpha \neq 1$ does not induce additional difficulties.

**Theorem 2.2.** Assume that (h1), (h2), (h3) and (h6) hold true and that $\eta > 0$.

(i) If $\theta \in C^1([0, T]; L^\infty(\Omega))$, $z_0 \in D(A_\eta^\kappa)$ and $f \in H^1(0, T; [L^2(\Omega)]^d)$, then there exists a unique strong solution $z \in C^1([0, T]; Z) \cap C^0([0, T]; D(A_\eta^\kappa))$ satisfying (19).
(ii) If $\theta \in C^0([0, T] \times \Omega)$, $z_0 \in Z$ and $f \in L^2(0, T; [L^2(\Omega)]^d)$, then Problem (19) has a unique mild solution $z \in C^0([0, T]; Z)$ such that $z(0) = z_0$ and

$$ \int_{0}^{T} z(t) \psi(t) \, dt \in D(A_{\eta}^{\kappa}), \qquad (28) $$

$$ -\int_{0}^{T} z(t) \dot{\psi}(t) \, dt + A_{\eta}^{\kappa} \left( \int_{0}^{T} z(t) \psi(t) \, dt \right) + \int_{0}^{T} G(t) z(t) \psi(t) \, dt = \int_{0}^{T} g(t) \psi(t) \, dt, \qquad (29) $$

for all $\psi \in C_c^1([0, T]; \mathbb{R})$. Moreover, $z$ verifies the Duhamel formula

$$ z(t) = \Phi_{\eta}^{\kappa}(t)z_0 + \int_{0}^{t} \Phi_{\eta}^{\kappa}(t-s) (-G(s)z(s) + g(s)) \, ds, \qquad (30) $$

where $\Phi_{\eta}^{\kappa}$ denotes the continuous semigroup generated by $A_{\eta}^{\kappa}$ in the sense that

$$ A_{\eta}^{\kappa} x = -\frac{d}{dt} (\Phi_{\eta}^{\kappa}(t)x) \Big|_{t=0^{+}}, \quad x \in Z. \qquad (31) $$

*Proof.* Let us prove (ii). We shall first show that the operator $A_{\eta}^{\kappa}$ defined by (23) is maximal-accretive, namely:

* $(A_{\eta}^{\kappa}z, z)_Z \ge 0, \quad \forall z \in D(A_{\eta}^{\kappa});$

* $A_{\eta}^{\kappa} + \lambda_0 I$ is surjective from $D(A_{\eta}^{\kappa})$ to $Z$, for all $\lambda_0 > 0$.

For any $z = (u_s, v_s, v_f, p) \in D(A_{\eta}^{\kappa})$, we have by definition of the bilinear form $a_{\eta}^{\kappa}$ and the operator $A_{\eta}^{\kappa}$

$$ (A_{\eta}^{\kappa} z, z)_Z = a(z, z) = 2\eta \int_{\Omega} |\varepsilon(v_s)|^2 \, dx + \int_{\Omega} \phi^2 k_f^{-1} (v_f - v_s) \cdot (v_f - v_s) \, dx + \int_{\Omega} \phi \sigma_f (v_f) : \varepsilon(v_f) \, dx. $$

Since $k_f^{-1}(v_f - v_s) \cdot (v_f - v_s) \ge 0$, we find that $(A_{\eta}^{\kappa}z, z)_Z \ge 0$.

Let $\lambda_0 > 0$ be a positive real number and let $g$ be an element of $Z$. To prove that $A_{\eta}^{\kappa} + \lambda_0 I$ is surjective from $D(A_{\eta}^{\kappa})$ to $Z$, we consider the variational problem

$$ \left\{ \begin{array}{l} \text{Find } z \in Y \text{ such that} \\ \forall y \in Y, \quad a_{\eta}^{\kappa}(z, y) + \lambda_0(z, y)_Z = (g, y)_Z. \end{array} \right. \qquad (32) $$

Using Poincaré inequality, we see that the linear form $y \mapsto (g, y)_Z$ is continuous over $Y$ and that the bilinear form $a_{\eta}^{\kappa}(\cdot, \cdot) + \lambda_0(\cdot, \cdot)_Z$ is continuous over $Y \times Y$. Moreover,

$$ \begin{aligned} a_{\eta}^{\kappa}(z, z) + \lambda_0(z, z)_Z = & 2\eta \int_{\Omega} |\varepsilon(v_s)|^2 \, dx + \int_{\Omega} \phi^2 k_f^{-1} (v_f - v_s) \cdot (v_f - v_s) \, dx + \int_{\Omega} \phi \sigma_f (v_f) : \varepsilon(v_f) \, dx \\ & + \lambda_0 \left( \|u_s\|_s^2 + \int_{\Omega} \rho_s (1-\phi) |v_s|^2 \, dx + \int_{\Omega} \rho_f \phi |v_f|^2 \, dx + \int_{\Omega} \frac{\alpha-\phi}{\kappa} p^2 \, dx \right) \\ & \ge \lambda_0 \|u_s\|_s^2 + 2\eta \|\varepsilon(v_s)\|^2 + 2\mu_f \phi_{\min} \|\varepsilon(v_f)\|^2 + \lambda_0 \frac{(\alpha-\phi)_{\min}}{\kappa} \|p\|^2, \end{aligned} \qquad (33) $$

where $\|\cdot\|$ denotes the $L^2$ norm indifferently in $[L^2(\Omega)]^d$ or $L^2(\Omega)$. Consequently, the bilinear form $a_{\eta}^{\kappa}(\cdot, \cdot) + \lambda_0(\cdot, \cdot)_Z$ is coercive on $Y$ thanks to Korn inequality (12).

From Lax-Milgram theorem, we deduce that there exists a unique $z \in Y$ solution of (32). Since by construction $a_{\eta}^{\kappa}(z, y) = (g - \lambda_0 z, y)_Z$ for all $y \in Y$ and $g - \lambda_0 z \in Z$, we finally get that $z \in D(A_{\eta}^{\kappa})$ in view of (24).
Hence, $A_{\eta}^{\kappa}$ is maximal-accretive and Lumer-Phillips theorem (see for instance [82, Chapter 1, Theorem 4.3]) implies that $A_{\eta}^{\kappa}$ is the infinitesimal generator – in the sense of (31) – of a $C^0$-semigroup of contraction $(\Phi_{\eta}^{\kappa}(t))_{t \ge 0}$. In particular, we have

$$ \|\Phi_{\eta}^{\kappa}(t)\|_{\mathcal{L}(Z)} \le 1, \quad t \in [0, T]. \qquad (34) $$

Then, we observe that $G(t)$ is a bounded perturbation of $A_{\eta}^{\kappa}$. Indeed, for any $z \in Z$,

$$ \|G(t)z\|_{Z}^{2} = \int_{\Omega} \rho_{f} \phi \left( \frac{\theta(t)}{\rho_{f} \phi} \right)^{2} |v_{f}|^{2} \, dx \le \omega^{2} \|z\|_{Z}^{2}, $$

with $(\rho_f \phi_{\min})^{-1} \omega = \|\theta\|_{L^{\infty}((0,T) \times \Omega)}$. Thus $G \in C^0([0, T]; \mathcal{L}(Z))$ and

$$ \|G(t)\|_{\mathcal{L}(Z)} \le \omega, \quad t \in [0, T]. \qquad (35) $$

Therefore, the assertion (ii) follows from [14, Part II, Chapter 1, Proposition 3.4] and [36, Corollary 2.19].

If $\theta \in C^1([0, T]; L^\infty(\Omega))$ then $G \in C^1([0, T]; \mathcal{L}(Z))$, which proves (i) by an application of [14, Part II, Chapter 1, Proposition 3.5]. $\square$

**Remark 2.3.** The bilinear form $a_{\eta}^{\kappa}(\cdot, \cdot) + \lambda_0(\cdot, \cdot)_Z$ is coercive on $Y$ precisely because $\eta > 0$. It will not be the case when $\eta = 0$. In particular, this implies that, here, $(\Phi_{\eta}^{\kappa}(t))_{t \ge 0}$ is an analytic semigroup [14, Part II, Chapter 1, Theorem 2.12].

The solution $z \in C^1([0, T]; Z) \cap C^0([0, T]; D(A_{\eta}^{\kappa}))$, called strong solution in the foregoing, is sometimes referred to as *strict* solution to account for the $C^1$ regularity in time – see for instance [14, Part II, Chapter 1, Definition 3.1]. The next result clarifies in which sense this solution satisfies the original equation under study.

**Corollary 2.4.** If $\theta \in C^1([0, T]; L^\infty(\Omega))$, $z_0 \in D(A_{\eta}^{\kappa})$ and $f \in H^1(0, T; [L^2(\Omega)]^d)$, then the strong solution defined above satisfies (18) almost everywhere in $(0, T) \times \Omega$.

*Proof.* The strong solution satisfies

$$ \begin{cases} \dot{z}(t) + A_{\eta}^{\kappa}z(t) + G(t)z(t) = g(t), & t \in [0, T], \\ z(0) = z_0. \end{cases} $$

Since $z \in C^1([0, T]; Z)$, we have $\partial_t u_s \in C^0([0, T; [H_0^1(\Omega)]^d])$, $\partial_t v_s \in C^0([0, T; [L^2(\Omega)]^d)$, $\partial_t v_f \in C^0([0, T; [L^2(\Omega)]^d)$ and $\partial_t p \in C^0([0, T; L^2(\Omega))$. In view of (27), the regularity $z \in C^0([0, T]; D(A_{\eta}^{\kappa}))$ implies that

$$ -\mathrm{div}(\sigma_s(u_s)) - \mathrm{div}(\sigma_s^{\mathrm{vis}}(v_s)) + (\alpha - \phi)\nabla p \in C^0([0, T]; [L^2(\Omega)]^d), \quad -\mathrm{div}(\phi\sigma_f(v_f)) + \phi\nabla p \in C^0([0, T]; [L^2(\Omega)]^d) $$

$$ \text{and } \mathrm{div}((\alpha - \phi)v_s + \phi v_f) \in C^0([0, T]; L^2(\Omega)). \text{ Thus for every } t \in [0, T], \text{(18b), (18c) and (18d) are verified in } [L^2(\Omega)]^d, \text{ and in particular almost everywhere. } \square $$

In other words, Corollary 2.4 does not mean that each individual term appearing in (18) is defined almost everywhere. However, each line of (18) is satisfied almost everywhere since $\dot{z}$ and $A_{\eta}^{\kappa}z$ are both defined almost everywhere.

**Theorem 2.2** provides the existence and uniqueness of two types of solutions: the *strong* solution and the *mild* solution. The strong solution is regular since it belongs to $C^1([0, T]; Z) \cap C^0([0, T]; D(A_{\eta}^{\kappa}))$ but it requires high regularity assumptions on the source terms and on the initial conditions, in particular $z_0 \in D(A_{\eta}^{\kappa})$. The mild solution requires weaker assumptions, but the Duhamel formula (30) is quite abstract. The next theorem establishes the existence and uniqueness of a third notion of solution: the *variational* solution, that satisfies a weak formulation in the following sense.
**Theorem 2.5.** Assume that (h1) – (h6) hold true and that $\eta > 0$. If $z_0 = (u_{s0}, v_{s0}, v_{f0}, p_0) \in Z$, then there exists a variational solution $u_s \in C^0([0, T]; [H_0^1(\Omega)]^d)$, $\partial_t u_s \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$, $v_f \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ and $p \in C^0([0, T]; L^2(\Omega))$ such that

$$ (u_s(0), \partial_t u_s(0), v_f(0), p(0)) = (u_{s0}, v_{s0}, v_{f0}, p_0) \quad (36) $$

and such that the following equations hold true in $D'(0, T)$:

$$ \left\{ \begin{array}{l} \forall (w_s, w_f, q) \in [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times L^2(\Omega), \\ \frac{d^2}{dt^2} \int_{\Omega} \rho_s(1-\phi) u_s(t) \cdot w_s \,dx + \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(w_s) \,dx + 2\eta \int_{\Omega} \varepsilon(\partial_t u_s(t)) : \varepsilon(w_s) \,dx \\ \qquad - \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - \partial_t u_s(t)) \cdot w_s \,dx - \int_{\Omega} p(t) \operatorname{div}((\alpha - \phi) w_s) \,dx = \int_{\Omega} \rho_s(1-\phi)f(t) \cdot w_s \,dx, \quad (37a) \end{array} \right. $$

$$ \left\{ \begin{array}{l} \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \,dx + \int_{\Omega} \phi \sigma_f(v_f(t)) : \varepsilon(w_f) \,dx + \int_{\Omega} \phi^2 k_f^{-1} (v_f(t) - \partial_t u_s(t)) \cdot w_f \,dx \\ \qquad - \int_{\Omega} \theta(t) v_f(t) \cdot w_f \,dx - \int_{\Omega} p(t) \operatorname{div}(\phi w_f) \,dx = \int_{\Omega} \rho_f \phi f(t) \cdot w_f \,dx, \quad (37b) \end{array} \right. $$

$$ \frac{d}{dt} \int_{\Omega} \frac{\alpha - \phi}{\kappa} p(t)q \,dx + \int_{\Omega} \operatorname{div}((\alpha - \phi) \partial_t u_s(t) + \phi v_f(t))q \,dx = \int_{\Omega} \frac{\theta(t)}{\rho_f} q \,dx. \quad (37c) $$

Furthermore, the energy estimate (11) holds true and, if we assume that $\partial_t \theta \in L^\infty((0, T) \times \Omega)$, this solution is unique.

**Proof.** To show the existence of variational solutions verifying (37), we proceed as follows. First, we approximate the data by sequences of regular functions and we consider the sequence of strong solutions associated with these regular data. Then, we show that these strong solutions satisfy a variational formulation and we pass to the limit on this formulation after having established some a priori estimates and strong convergences of the sequences.

As $A_\eta^\kappa$ is maximal, $D(A_\eta^\kappa)$ is dense in $Z$. Let $z_0^n$ be a sequence of elements of $D(A_\eta^\kappa)$ converging towards $z_0$ strongly in $Z$. Let $f^n$ denote a sequence of $H^1(0, T; L^2(\Omega))$ converging towards $f$ in $L^2(0, T; L^2(\Omega))$ and $\theta^n$ denote a sequence of $C^1([0, T]; L^\infty(\Omega))$ converging towards $\theta$ in $C^0([0, T] \times \Omega)$. From Theorem 2.2, we know that there exists a unique strong solution $z^n = (u_s^n, v_s^n, v_f^n, p^n) \in C^1([0, T]; Z) \cap C^0([0, T]; D(A_\eta^\kappa))$ to the problem

$$ \begin{cases} \dot{z}^n(t) + A_\eta^\kappa z^n(t) + G^n(t)z^n(t) = g^n(t), & t \in [0, T], \\ z^n(0) = z_0^n. \end{cases} \quad (38) $$

Multiplying (38) by $y = (d_s, w_s, w_f, p) \in Y$, we see from (23) that $(A_\eta^\kappa z^n(t), y)_Z = a_\eta^\kappa(z^n(t), y)$. Hence $z^n$ satisfies the following variational formulation: for all $s \in [0, T]$,

$$ (VF)^n \left\{ \begin{array}{l} \forall (d_s, w_s, w_f, p) \in Y = [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times L^2(\Omega), \\ \int_{\Omega} \sigma_s(\partial_t u_s^n(s)) : \varepsilon(d_s) \,dx = \int_{\Omega} \sigma_s(v_s^n(s)) : \varepsilon(d_s) \,dx, \\ \int_{\Omega} \rho_s(1-\phi) \partial_t v_s^n(s) \cdot w_s \,dx + \int_{\Omega} \sigma_s(u_s^n(s)) : \varepsilon(w_s) \,dx + 2\eta \int_{\Omega} \varepsilon(v_s^n(s)) : \varepsilon(w_s) \,dx \\ \qquad - \int_{\Omega} \phi^2 k_f^{-1}(v_f^n(s) - v_s^n(s)) \cdot w_s \,dx - \int_{\Omega} p^n(s) \operatorname{div}((\alpha - \phi) w_s) \,dx = \int_{\Omega} \rho_s(1-\phi) f^n(s) \cdot w_s \,dx, \\ \int_{\Omega} \rho_f \phi \partial_t v_f^n(s) \cdot w_f \,dx + \int_{\Omega} \phi \sigma_f(v_f^n(s)) : \varepsilon(w_f) \,dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f^n(s) - v_s^n(s)) \cdot w_f \,dx \\ \qquad - \int_{\Omega} \theta^n(s) v_f^n(s) \cdot w_f \,dx - \int_{\Omega} p^n(s) \operatorname{div}(\phi w_f) \,dx = \int_{\Omega} \rho_f \phi f^n(s) \cdot w_f \,dx, \\ \int_{\Omega} \frac{\alpha - \phi}{\kappa} \partial_t p^n(s) q \,dx + \int_{\Omega} \operatorname{div}((\alpha - \phi) v_s^n(s) + \phi v_f^n(s)) q \,dx = \int_{\Omega} \frac{\theta^n(s)}{\rho_f} q \,dx. \end{array} \right. $$
Recalling that $z^n \in C^0([0, T]; D(A_\eta^k)) \subset C^0([0, T]; Y)$, we can choose $d_s = u_s^n(s)$, $w_s = v_s^n(s)$, $w_f = v_f^n(s)$ and $q = p^n(s)$ as test functions. Integrating in time from 0 to $t$ and applying Grönwall Lemma like in Section 1, we get the energy inequality

$$ 
\begin{aligned}
\frac{1}{2} \int_{\Omega} \sigma_s(u_s^n(t)) : \varepsilon(u_s^n(t)) \, dx &+ \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |v_s^n(t)|^2 \, dx + 2\eta \int_0^t \int_{\Omega} \varepsilon(v_s^n) : \varepsilon(v_s^n) \, dx \, ds \\
&+ \frac{\rho_f}{2} \int_{\Omega} \phi |v_f^n(t)|^2 \, dx + \int_0^t \int_{\Omega} \phi \sigma_f(v_f^n) : \varepsilon(v_f^n) \, dx \, ds + \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |p^n(t)|^2 \, dx \\
&\le \exp\left(\max\left(1, \frac{2\|\theta^n\|_{C^0([0,T]\times\Omega)}}{\rho_f \phi_{\min}}\right) t\right) \left( C \int_{\Omega} |f^n(t)|^2 \, dx \, ds + C \int_0^t \int_{\Omega} |\theta^n(t)|^2 \, dx \, ds + \frac{1}{2} \int_{\Omega} \sigma_s(u_s^n) : \varepsilon(u_s^n) \, dx \right. \\
&\qquad \left. + \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |v_s^n(0)|^2 \, dx + \frac{\rho_f}{2} \int_{\Omega} \phi |v_f^n(0)|^2 \, dx + \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |p_0^n(t)|^2 \, dx \right). \quad (39)
\end{aligned}
$$

Thanks to the assumptions done on the data, the right-hand side of the latter inequality is uniformly bounded with respect to $n$. Consequently, taking into account the assumptions (h2) and (h6) on $\phi$ and $\alpha$, Korn inequality (12), the coercivity of the elastic and fluid forms (13) and (14), we deduce that $u_s^n$ is uniformly bounded in $C^0([0, T]; [H_0^1(\Omega)]^d)$, $v_s^n$ in $C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$, $v_f^n$ in $C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ and $p^n$ in $C^0([0, T]; L^2(\Omega))$.

In the same way, one can show that $z^n$ is a Cauchy sequence in $C^0([0, T]; [H_0^1(\Omega)]^d) \times (C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)) \times (C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)) \times C^0([0, T]; L^2(\Omega))$. Indeed, denoting $z^{n,m} = z^n - z^m$, using the linearity of the coupled problem, the uniform bound we just obtained and taking into account the fact that $\|\theta^n\|_{L^\infty((0,T)\times\Omega)}$ is bounded uniformly in $n$, we obtain that there exists $C > 0$ independent of $n$ such that

$$ \frac{d}{dt} \chi^{n,m}(t) \le C \chi^{n,m}(t) + \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |f^{n,m}(t)|^2 \, dx + \frac{\rho_f}{2} \int_{\Omega} \phi |f^{n,m}(t)|^2 \, dx + C \|\theta^{n,m}\|_{L^\infty((0,T)\times\Omega)}, $$

with

$$ 
\begin{aligned}
\chi^{n,m}(t) = \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |v_s^{n,m}(t)|^2 \, dx &+ \frac{1}{2} \int_{\Omega} \sigma_s(u_s^{n,m}(t)) : \varepsilon(u_s^{n,m}(t)) \, dx + \frac{\rho_f}{2} \int_{\Omega} \phi |v_f^{n,m}(t)|^2 \, dx \\
&+ \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |p^{n,m}(t)|^2 \, dx + 2\eta \int_0^t \int_{\Omega} \varepsilon(v_s^{n,m}) : \varepsilon(v_s^{n,m}) \, dx \, ds + \int_0^t \int_{\Omega} \phi \sigma_f(v_f^{n,m}) : \varepsilon(v_f^{n,m}) \, dx \, ds.
\end{aligned}
$$

Grönwall Lemma and the fact that the sequences associated with the data are Cauchy sequences imply that

$$ \forall \varepsilon > 0, \exists N \in \mathbb{N}, \forall n \ge N, \forall m \ge N, \chi^{n,m}(t) \le \varepsilon \exp(Ct). $$

As a consequence, there exists $u_s \in C^0([0, T]; [H_0^1(\Omega)]^d)$, $v_s \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$, $v_f \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ and $p \in C^0([0, T]; L^2(\Omega))$ such that

$$ 
\begin{aligned}
u_s^n &\to u_s \quad \text{in } C^0([0, T]; [H_0^1(\Omega)]^d), \quad v_s^n &\to v_s \quad \text{in } C^0([0, T]; [L^2(\Omega)]^d), \quad v_f^n &\to v_f \quad \text{in } C^0([0, T]; [L^2(\Omega)]^d), \\
p^n &\to p \quad \text{in } C^0([0, T]; L^2(\Omega)), \quad v_s^n &\to v_s \quad \text{in } L^2(0, T; [H_0^1(\Omega)]^d), \quad v_f^n &\to v_f \quad \text{in } L^2(0, T; [H_0^1(\Omega)]^d).
\end{aligned}
$$

Note that since $v_s^n = \partial_t u_s^n$, it also holds true in the limit and $\partial_t u_s = v_s \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$. These convergences enable to pass to the limit in $(VF)^n$, except for the inertial terms. Yet, these terms can be rewritten thanks to an integration by parts in time as follows: for $\psi \in \mathcal{D}(0, T)$, we have for example

$$ 
\begin{aligned}
\int_0^T \int_{\Omega} \rho_s (1-\phi) \partial_t v_s^n(t) \cdot w_s \psi(t) \, dx \, dt &= - \int_0^T \int_{\Omega} \rho_s (1-\phi) v_s^n(t) \cdot w_s \dot{\psi}(t) \, dx \, dt \\
&\longrightarrow_{n \to \infty} - \int_0^T \int_{\Omega} \rho_s (1-\phi) v_s(t) \cdot w_s \dot{\psi}(t) \, dx \, dt.
\end{aligned}
$$
By similar arguments and thanks to the strong convergences, we get, in $\mathcal{D}'(0, T)$,

$$ \forall (d_s, w_s, w_f, p) \in Y, \qquad \frac{d}{dt} \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(d_s) \, dx = \int_{\Omega} \sigma_s(v_s(t)) : \varepsilon(d_s) \, dx, \qquad (40a) $$

$$ \frac{d}{dt} \int_{\Omega} \rho_s(1-\phi)v_s(t) \cdot w_s \, dx + \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(w_s) \, dx + 2\eta \int_{\Omega} \varepsilon(v_s(t)) : \varepsilon(w_s) \, dx \\ - \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - v_s(t)) \cdot w_s \, dx - \int_{\Omega} p(t) \operatorname{div}((\alpha - \phi)w_s) \, dx = \int_{\Omega} \rho_s(1-\phi)f(t) \cdot w_s \, dx, \quad (40b) $$

$$ \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \, dx + \int_{\Omega} \phi \sigma_f(v_f(t)) : \varepsilon(w_f) \, dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - v_s(t)) \cdot w_f \, dx \\ - \int_{\Omega} \theta(t) v_f(t) \cdot w_f \, dx - \int_{\Omega} p(t) \operatorname{div}(\phi w_f) \, dx = \int_{\Omega} \rho_f \phi f(t) \cdot w_f \, dx, \quad (40c) $$

$$ \frac{d}{dt} \int_{\Omega} \frac{\alpha - \phi}{\kappa} p(t) q \, dx + \int_{\Omega} \operatorname{div}((\alpha - \phi) v_s(t) + \phi v_f(t)) q \, dx = \int_{\Omega} \frac{\theta(t)}{\rho_f} q \, dx. \quad (40d) $$

To obtain (37), it only remains to rewrite (40a) and (40b) as a second order equation in time, which holds true since $v_s = \partial_t u_s$ in $C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$. Lastly, we recover the initial conditions (36) by simply passing to the limit in the second line of (38).

To ensure uniqueness, we observe that every variational solution satisfying the energy estimate is unique. Indeed, for $f = 0$, $\theta = 0$ and $z_0 = 0$, we obtain $z = 0$ in virtue of (11). Therefore, it is sufficient to prove that every variational solution satisfying (37) verifies the energy identity (10) and thus the energy estimate. To do so, let us first derive a bound on $(\partial_t u_s, \partial_t v_s, \partial_t v_f, \partial_t p)$. From (40), we deduce

$$ \forall y \in Y, \quad - \int_{0}^{T} (z(t), y)_{Z} \dot{\psi}(t) \, dt + \int_{0}^{T} a_{\eta}^{\kappa}(z(t), y) \psi(t) \, dt \\ \qquad + \int_{0}^{T} (G(t)z(t), y)_{Z} \psi(t) \, dt = \int_{0}^{T} (g(t), y)_{Z} \psi(t) \, dt. \quad (41) $$

Since $f \in L^2(0, T; [L^2(\Omega)]^d)$, $\theta \in C^0((0, T) \times \Omega)$ and by continuity of the bilinear form $a_{\eta}^{\kappa}$ over $Y \times Y$, we have

$$ \forall y \in Y, \quad - \int_{0}^{T} (z(t), y)_{Z} \dot{\psi}(t) \, dt = \int_{0}^{T} (h(t), y)_{Z} \psi(t) \, dt \quad \text{with } h \in L^2(0, T; Y'). $$

Thus $\dot{z} \in L^2(0, T; Y')$, namely

$$ (\partial_t u_s, \partial_t v_s, \partial_t v_f, \partial_t p) \in L^2(0, T; [H^{-1}(\Omega)]^d) \times L^2(0, T; [H^{-1}(\Omega)]^d) \times L^2(0, T; [H^{-1}(\Omega)]^d) \times L^2(0, T; L^2(\Omega)). $$

Since $\partial_t u_s = v_s$, finally, it holds

$$ \begin{align} \partial_t u_s \in L^2(0, T; [H_0^1(\Omega)]^d) \quad \text{and} \quad \partial_{tt} u_s \in L^2(0, T; [H^{-1}(\Omega)]^d), \notag \\ \partial_t u_s \in L^2(0, T; [H_0^1(\Omega)]^d) \quad \text{and} \quad \partial_t v_f \in L^2(0, T; [H^{-1}(\Omega)]^d), \tag{42} \\ p \in L^2(0, T; L^2(\Omega)) \quad \text{and} \quad \partial_t p \in L^2(0, T; L^2(\Omega)). \notag \end{align} $$

Using a standard result of functional analysis (see for instance [54, Chapter XVIII, Proposition 7]), the previous regularities imply that the following relations hold in $\mathcal{D}'(0, T)$:

$$ \forall w_s \in [H_0^1(\Omega)]^d, \quad \frac{d^2}{dt^2} \int_{\Omega} \rho_s(1-\phi) u_s(t) \cdot w_s \, dx = \langle \rho_s(1-\phi) \partial_{tt} u_s(t), w_s \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d}, $$

$$ \forall w_f \in [H_0^1(\Omega)]^d, \quad \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \, dx = \langle \rho_f \phi \partial_t v_f(t), w_f \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d}, $$

$$ \forall q \in L^2(\Omega), \quad \frac{d}{dt} \int_{\Omega} \frac{\alpha - \phi}{\kappa} p(t) q \, dx = \int_{\Omega} \frac{\alpha - \phi}{\kappa} \partial_t p(t) q \, dx. $$
Moreover, since functions in $[H_0^1(\Omega)]^d \otimes \mathcal{D}(0, T)$ and $L^2(\Omega) \otimes \mathcal{D}(0, T)$ generate respectively $L^2(0, T; [H_0^1(\Omega)]^d)$ and $L^2(0, T; L^2(\Omega))$, we obtain the space-time variational formulation

$$ 
\left\{
\begin{array}{l}
\forall (w_s, w_f, q) \in L^2(0, T; [H_0^1(\Omega)]^d) \times L^2(0, T; [H_0^1(\Omega)]^d) \times L^2(0, T; L^2(\Omega)), \\
\int_0^T \left\langle \rho_s(1-\phi) \partial_{tt} u_s, w_s \right\rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} \mathrm{d}t + \int_0^T \int_{\Omega} \sigma_s(u_s) : \varepsilon(w_s) \mathrm{d}x \mathrm{d}t + 2\eta \int_0^T \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(w_s) \mathrm{d}x \mathrm{d}t \\
\quad - \int_0^T \int_{\Omega} \phi^2 k_f^{-1}(v_f - \partial_t u_s) \cdot w_s \mathrm{d}x \mathrm{d}t - \int_0^T \int_{\Omega} p \operatorname{div}((\alpha - \phi) w_s) \mathrm{d}x \mathrm{d}t = \int_0^T \int_{\Omega} \rho_s(1-\phi)f \cdot w_s \mathrm{d}x \mathrm{d}t, \\
\int_0^T \left\langle \rho_f \phi \partial_t v_f, w_f \right\rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} \mathrm{d}t + \int_0^T \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(w_f) \mathrm{d}x \mathrm{d}t + \int_0^T \int_{\Omega} \phi^2 k_f^{-1}(v_f - \partial_t u_s) \cdot w_f \mathrm{d}x \mathrm{d}t \\
\quad - \int_0^T \int_{\Omega} \theta v_f \cdot w_f \mathrm{d}x \mathrm{d}t - \int_0^T \int_{\Omega} p \operatorname{div}(\phi w_f) \mathrm{d}x \mathrm{d}t = \int_0^T \int_{\Omega} \rho_f \phi f \cdot w_f \mathrm{d}x \mathrm{d}t, \\
\int_0^T \int_{\Omega} \frac{\alpha - \phi}{\kappa} \partial_t p q \mathrm{d}x \mathrm{d}t + \int_0^T \int_{\Omega} \operatorname{div}((\alpha - \phi) \partial_t u_s + \phi v_f) q \mathrm{d}x \mathrm{d}t = \int_0^T \int_{\Omega} \frac{\theta}{\rho_f} q \mathrm{d}x \mathrm{d}t.
\end{array}
\right.
$$

Now, since we know that $\partial_t u_s \in L^2(0, T; [H_0^1(\Omega)]^d)$, $v_f \in L^2(0, T; [H_0^1(\Omega)]^d)$ and $p \in L^2(0, T; L^2(\Omega))$, we can choose $(w_s, w_f, q) = (\partial_t u_s, v_f, p)$ as test functions in the above formulation, which provides the energy identity (10) and thus the energy estimate (11). $\square$

**Remark 2.6.** The method used to prove Theorem 2.5 is standard and close to the Faedo-Galerkin method. The difference with Faedo-Galerkin method is that the approximated sequence is directly recovered from the existence of strong solutions instead of being constructed on a suitable finite dimensional space. This allows us to obtain strong convergence for the whole sequence, whereas Faedo-Galerkin method provides only weak convergence of subsequences. In addition, it directly provides the continuity with respect to time of the solution and the strong convergence of the initial condition $z(0) = z_0$ in $Z$.

**Remark 2.7.** The variational solution could also be defined without assuming that it is continuous with respect to time, but only assuming that the regularities (42) are satisfied. The time continuity of the solution can then be recovered using the existence of a continuous and linear mapping of the space

$$ W(0, T) = \{u \in L^2(0, T; [H_0^1(\Omega)]^d) \text{ such that } \partial_t u \in L^2(0, T; [H^{-1}(\Omega)]^d)\} $$
into $C^0([0, T]; [L^2(\Omega)]^d)$, see [71, Chapter 1, Theorem 3.1].

The mild solution and the variational solution are two notions of solution whose existence and uniqueness here require the same hypotheses on the data. In fact, the following result states that these two types of solution are the same whenever $f \in L^2((0, T) \times \Omega)$ and $\theta \in C^0([0, T] \times \Omega)$. Thus, they can be used indifferently depending on the context. For instance, the mild solution is widely used in control theory because of the practical aspects of Duhamel formula, whereas the variational solution formulation is usually the one implemented at the discrete level when considering finite element discretization.

**Proposition 2.8.** If $f \in L^2((0, T) \times \Omega)$ and $\theta \in C^0([0, T] \times \Omega)$, then the mild solution given by (30) and the variational solution satisfying (37) coincide.

*Proof.* The mild solution and the variational solution are both unique. Hence, it is sufficient to show that the variational solution defined in Theorem 2.5 is also a mild solution, namely that it satisfies (28) and (29). Let $\psi$ be given in $\mathcal{D}(0, T)$. Since $v_s \in L^2(0, T; [H_0^1(\Omega)]^d)$ and $v_f \in L^2(0, T; [H_0^1(\Omega)]^d)$, it holds that

$$ a_{\eta}^{\kappa} \left( \int_{0}^{T} z(t) \psi(t) \mathrm{d}t, y \right) = \int_{0}^{T} a_{\eta}^{\kappa}(z(t), y) \psi(t) \mathrm{d}t $$

for all $y \in Y$, so that we can rewrite (41) as

$$ \forall y \in Y, \quad - \int_{0}^{T} (z(t), y)_{Z} \dot{\psi}(t) \mathrm{d}t + a_{\eta}^{\kappa} \left( \int_{0}^{T} z(t) \psi(t) \mathrm{d}t, y \right) + \int_{0}^{T} (G(t) z(t), y)_{Z} \psi(t) \mathrm{d}t = \int_{0}^{T} (g(t), y)_{Z} \psi(t) \mathrm{d}t. $$
From the definition of $D(A_{\eta}^{\kappa})$, it follows that

$$ \int_{0}^{T} z(t)\psi(t) \, dt \in D(A_{\eta}^{\kappa}) $$

and thus

$$ -\int_{0}^{T} (z(t), y)_{Z} \dot{\psi}(t) \, dt + \left( A_{\eta}^{\kappa} \left( \int_{0}^{T} z(t)\psi(t) \, dt \right), y \right)_{Z} + \int_{0}^{T} (G(t)z(t), y)_{Z} \psi(t) \, dt = \int_{0}^{T} (g(t), y)_{Z} \psi(t) \, dt \quad (44) $$

for all $y \in Y$. As $Y$ is dense in $Z$, (44) is also true for all $y \in Z$, which proves (29). □

**Remark 2.9.** Note that the mild and variational solutions coincide only under the assumption $f \in L^2((0, T) \times \Omega)$ and $\theta \in C^0([0, T] \times \Omega)$. If $f \in L^2(0, T; [H^{-1}(\Omega)]^d)$, we can readily extend the existence and uniqueness of variational solutions proved in Theorem 2.5, but the existence of a mild solution is not guaranteed.

## 2.3 The case $\eta = 0$

Without solid viscosity, the solid formulation becomes hyperbolic. This hyperbolic-parabolic coupling was studied in [10], where existence of variational solutions is derived. In [10], the fluid mass input $\theta$ is supposed to be small enough, namely there exists $C_f > 0$ such that

$$ \forall v_f \in [H_0^1(\Omega)]^d, \quad \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) \, dx - \int_{\Omega} \theta |v_f|^2 \, dx \geq C_f \|v_f\|^2_{[H_0^1(\Omega)]^d}. $$

Here, we are going to prove existence results of strong and mild solutions thanks to semigroup theory and deduce existence of variational solutions directly, without any smallness assumption on $\theta$.

The main issue in this case is that the underlying bilinear form is not coercive. Indeed if $\eta = 0$, then the bilinear form introduced in the proof of Theorem 2.2 is no more coercive on $Y$ in view of (33). Despite this lack of coercivity, we are going to show that Problem (32) is still well-posed when $\eta = 0$. To do so, we use the $T$-coercivity approach [50, 47], which is an alternative to Banach-Nečas-Babuška theory and that has been designed especially for non-coercive problems. For the sake of completeness, the definition and properties of $T$-coercivity are recalled below.

**Definition 2.1.** Let $V$ be an Hilbert space and let $a(\cdot, \cdot)$ be a continuous bilinear form over $V \times V$. We say that a is $T$-coercive if there exists a bijective application $T \in \mathcal{L}(V)$ and $\alpha > 0$ such that

$$ |a(z, Tz)| \geq \alpha \|z\|_{V}^{2}, \quad z \in V. $$

**Proposition 2.10.** Let $V$ be an Hilbert space. Let $\ell(\cdot)$ be a continuous linear form over $V$ and $a(\cdot, \cdot)$ be a continuous bilinear form over $V \times V$. The problem

$$ \begin{cases} \text{Find } z \in V \text{ such that} \\ \forall y \in V, \quad a(z, y) = \ell(y) \end{cases} $$

is well-posed if and only if a is $T$-coercive.

The following theorem states existence and uniqueness of solutions to Problem (19) in the case $\eta = 0$.

**Theorem 2.11.** If $\eta = 0$, then the conclusions of Theorem 2.2 remain true.

*Proof.* Let us show that $A_0^\kappa$ is maximal by proving that the variational problem

$$ \begin{cases} \text{Find } z \in Y \text{ such that} \\ \forall y \in Y, \quad a_0^\kappa(z, y) + \lambda_0(z, y)z = (g, y)z \end{cases} $$
is well-posed, where

$$
\begin{align*}
a_0^\kappa(z, y) + \lambda_0(z, y)_Z &= \int_\Omega \lambda_0 \sigma_s(u_s) : \varepsilon(d_s) \, dx - \int_\Omega \sigma_s(v_s) : \varepsilon(d_s) \, dx \\
&\quad + \int_\Omega \lambda_0 \rho_s (1-\phi) v_s \cdot w_s \, dx + \int_\Omega \sigma_s(u_s) : \varepsilon(w_s) \, dx + \int_\Omega \phi^2 k_f^{-1} (v_f - v_s) \cdot (w_f - w_s) \, dx \\
&\quad + \int_\Omega \lambda_0 \rho_f \phi v_f \cdot w_f \, dx + \int_\Omega \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx + \int_\Omega \lambda_0 \frac{\alpha-\phi}{\kappa} pq \, dx \\
&\quad - \int_\Omega p \operatorname{div} ((\alpha-\phi) w_s + \phi w_f) \, dx + \int_\Omega \operatorname{div} ((\alpha-\phi) v_s + \phi v_f) \, q \, dx
\end{align*}
$$

for any $z = (u_s, v_s, v_f, p)$ and $y = (d_s, w_s, w_f, q)$ in $Y$. From Proposition 2.10, it is sufficient to show that $a_0^\kappa(\cdot, \cdot) + \lambda_0(\cdot, \cdot)_Z$ is T-coercive.

For a given $z$, we look for a $y^*$ depending continuously on $z$ such that $a_0^\kappa(z, y^*) + \lambda_0(z, y^*)_Z \ge \underline{\alpha}\|z\|_Y^2$ for some constant $\underline{\alpha} > 0$. Choosing $w_s^* = v_s$, $w_f^* = v_f$, $q^* = p$ and $d_s^* = d_f^* = 1$ in the form $\beta u_s + \gamma v_s$ yields

$$
\begin{align*}
a_0^\kappa(z, y^*) + \lambda_0(z, y^*)_Z &= \int_\Omega \lambda_0 \beta \sigma_s(u_s) : \varepsilon(u_s) \, dx + \int_\Omega \lambda_0 \gamma \sigma_s(u_s) : \varepsilon(v_s) \, dx \\
&\quad - \int_\Omega \beta \sigma_s(v_s) : \varepsilon(u_s) \, dx - \int_\Omega \gamma \sigma_s(v_s) : \varepsilon(v_s) \, dx \\
&\quad + \int_\Omega \lambda_0 \rho_s (1-\phi) |v_s|^2 \, dx + \int_\Omega \sigma_s(u_s) : \varepsilon(v_s) \, dx + \int_\Omega \phi^2 k_f^{-1} (v_f - v_s) \cdot (v_f - v_s) \, dx \\
&\quad + \int_\Omega \lambda_0 \rho_f \phi |v_f|^2 \, dx + \int_\Omega \phi \sigma_f(v_f) : \varepsilon(v_f) \, dx + \int_\Omega \lambda_0 \frac{\alpha-\phi}{\kappa} p^2 \, dx.
\end{align*}
$$

By setting $\beta = \frac{1}{2}$ and $\gamma = -\frac{1}{2\lambda_0}$, the terms of the form $\int_{\Omega} \sigma_s(u_s) : \varepsilon(v_s) \, dx$ vanish so that

$$
a_0^\kappa(z, y^*) + \lambda_0(z, y^*)_Z \geq \frac{\lambda_0}{2} \|u_s\|_s^2 + \frac{1}{2\lambda_0} \|v_s\|_s^2 + 2\mu_f \phi_{\min} \|\varepsilon(v_f)\|^2 + \lambda_0 \frac{(\alpha - \phi)_{\min}}{\kappa} \|p\|^2.
$$

Therefore, $a_0^\kappa(\cdot, \cdot) + \lambda_0(\cdot, \cdot)_Z$ is T-coercive for the mapping $T$ defined by

$$
T : (u_s, v_s, v_f, p) \mapsto \left( \frac{1}{2} u_s - \frac{1}{2\lambda_0} v_s, v_s, v_f, p \right), \quad (45)
$$

which is continuous and bijective on $Y$.

The remainder of the proof follows the very same lines as for the viscous case.

<figure><img src="image_3.png" /></figure>

Next we recover the existence of variational solutions from the existence of strong solutions. However, we obtain a variational formulation that slightly differs from (37) because of the hyperbolic-parabolic coupling between the solid and fluid equations.

**Theorem 2.12.** Assume that (h1)–(h6) hold true and that $\eta = 0$. If $z_0 = (u_{s0}, v_{s0}, v_{f0}, p_0) \in Z$, then there exists a variational solution $u_s \in C^0([0, T]; [H_0^1(\Omega)]^d)$, $\partial_t u_s \in C^0([0, T]; [L^2(\Omega)]^d)$, $v_f \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ and $p \in C^0([0, T]; L^2(\Omega))$ such that

$$
(u_s(0), \partial_t u_s(0), v_f(0), p(0)) = (u_{s0}, v_{s0}, v_{f0}, p_0) \quad (46)
$$
and the following equations hold true, in $\mathcal{D}'(0, T)$,

$$ \left\{ \begin{array}{l} \forall (w_s, w_f, q) \in [\mathbf{H}_0^1(\Omega)]^d \times [\mathbf{H}_0^1(\Omega)]^d \times \mathbf{L}^2(\Omega), \\ \frac{d^2}{dt^2} \int_{\Omega} \rho_s(1-\phi) u_s(t) \cdot w_s \,dx + \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(w_s) \,dx \\ \qquad - \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - \partial_t u_s(t)) \cdot w_s \,dx - \int_{\Omega} p(t) \operatorname{div}((\alpha-\phi)w_s) \,dx = \int_{\Omega} \rho_s(1-\phi)f(t) \cdot w_s \,dx, \quad (47a) \end{array} \right. $$

$$ \left\{ \begin{array}{l} \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \,dx + \int_{\Omega} \phi \sigma_f(v_f(t)) : \varepsilon(w_f) \,dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - \partial_t u_s(t)) \cdot w_f \,dx \\ \qquad - \int_{\Omega} \theta(t) v_f(t) \cdot w_f \,dx - \int_{\Omega} p(t) \operatorname{div}(\phi w_f) \,dx = \int_{\Omega} \rho_f \phi f(t) \cdot w_f \,dx, \quad (47b) \end{array} \right. $$

$$ \frac{d}{dt} \int_{\Omega} \frac{\alpha - \phi}{\kappa} p(t) q \,dx + \frac{d}{dt} \int_{\Omega} \operatorname{div}((\alpha - \phi) u_s(t)) q \,dx + \int_{\Omega} \operatorname{div}(\phi v_f(t)) q \,dx = \int_{\Omega} \frac{\theta(t)}{\rho_f} q \,dx. \quad (47c) $$

This variational solution is unique, and coincides with the mild solution. Furthermore, the energy estimate (11) with $\eta = 0$ holds true.

**Remark 2.13.** Theorem 2.12 sheds light on the influence of solid viscosity on the model. Since $\eta = 0$, $\partial_t u_s$ does not belong to $\mathbf{L}^2(0, T; [\mathbf{H}_0^1(\Omega)]^d)$ but only to $\mathbf{C}^0([0, T]; [\mathbf{L}^2(\Omega)]^d)$. For this reason, equations (37c) and (47c) are not similar because, when $\eta = 0$, the term $\operatorname{div}((1-\phi)\partial_t u_s(t))$ is not in $\mathbf{L}^2(\Omega)$ in the space variable. One has only $\operatorname{div}((1-\phi)\partial_t u_s(t)) \in \mathbf{C}^0([0, T]; \mathbf{H}^{-1}(\Omega))$. This confirms that viscoelastic effects have an impact on the regularity of the solution, as it was already observed for other linear or non-linear poroelastic models [94, 22, 105].

*Proof.* We follow the same steps as for the proof of Theorem 2.5. The input data are approximated by regular functions, a priori estimates are established for the approximated solutions and we pass to the limit on the variational formulation $(VF)^n$ with $\eta = 0$.

The estimate (39) still holds true even if $\eta = 0$ because $z^n \in \mathbf{C}^0([0, T]; D(A_\kappa^\kappa)) \subset \mathbf{C}^0([0, T]; Y)$, in particular $v_s^n \in \mathbf{C}^0([0, T]; [\mathbf{H}_0^1(\Omega)]^d)$, which justifies that $z^n$ is regular enough to reproduce the formal calculations made in Section 1. As previously, this estimate implies that $z^n$ is a Cauchy sequence in $\mathbf{C}^0([0, T]; Z)$. However, since $\eta = 0$, estimate (39) only implies that $v_f^n$ is a Cauchy sequence in $\mathbf{L}^2(0, T; [\mathbf{H}_0^1(\Omega)]^d)$. Hence, the convergence

$$ v_f^n \to v_f \quad \text{in } \mathbf{L}^2(0, T; [\mathbf{H}_0^1(\Omega)]^d) $$

is still valid but now $v_s^n = \partial_t u_s^n$ does not converge in $\mathbf{L}^2(0, T; [\mathbf{H}_0^1(\Omega)]^d)$ but only in $\mathbf{C}^0([0, T]; [\mathbf{L}^2(\Omega)]^d)$. This changes the way to pass to the limit on $(VF)^n$ and in particular in the first equation. Let $\psi$ be an element of $\mathcal{D}(0, T)$. For any $d_s \in [\mathbf{L}^2(\Omega)]^d$, we consider the unique solution $\eta_s \in [\mathbf{H}_0^1(\Omega)]^d$ of $-\operatorname{div}(\sigma_s(\eta_s)) = d_s$ as a test function, so that

$$ \begin{aligned} \int_0^T \int_{\Omega} \sigma_s(v_s^n(t)) : \varepsilon(\eta_s) \psi(t) \,dx \,dt &= \int_0^T \int_{\Omega} \varepsilon(v_s^n(t)) : \sigma_s(\eta_s) \psi(t) \,dx \,dt \\ &= \int_0^T \int_{\Omega} v_s^n(t) \cdot d_s \psi(t) \,dx \,dt \xrightarrow{n \to \infty} \int_0^T \int_{\Omega} v_s(t) \cdot d_s \psi(t) \,dx \,dt. \end{aligned} $$

This proves that for all $d_s \in [\mathbf{L}^2(\Omega)]^d$, we have

$$ \frac{d}{dt} \int_{\Omega} u_s(t) \cdot d_s \,dx - \int_{\Omega} v_s(t) \cdot d_s \,dx = 0 \quad (48) $$

and we recover that $v_s = \partial_t u_s$ in $\mathcal{D}'((0, T) \times \Omega)$. In particular, it holds that $\partial_t u_s \in \mathbf{C}^0([0, T]; [\mathbf{L}^2(\Omega)]^d)$.
We can obtain (47a) and (47b) in a similar way as for the viscous case. Finally, to get (47c), we observe that

$$ \int_0^T \int_\Omega \operatorname{div}((\alpha - \phi) v_s^n(t)) q \psi(t) \, dx \, dt = - \int_0^T \int_\Omega \operatorname{div}((\alpha - \phi) u_s^n(t)) q \dot{\psi}(t) \, dx \, dt \\ \xrightarrow{n \to \infty} - \int_0^T \int_\Omega \operatorname{div}((\alpha - \phi) u_s(t)) q \dot{\psi}(t) \, dx \, dt, $$

where we have integrated by parts in time and used that $u_s^n$ converges in $C^0([0, T]; [H_0^1(\Omega)]^d)$.

In the same way as for the viscous case, we notice that the weak formulation (47) provides some regularity on the time derivative of the solution. For instance, the fluid equation (47b) implies that for any $\psi \in D(0, T)$ and $w_f \in [H_0^1(\Omega)]^d$, we have

$$ - \int_0^T \int_\Omega \rho_f \phi v_f(t) \cdot w_f \dot{\psi}(t) \, dx \, dt = \int_0^T \int_\Omega F(t) \cdot w_f \psi(t) \, dx \, dt, $$

with $F = \rho_f \phi f + \operatorname{div}(\phi \sigma_f(v_f)) - \phi^2 k_f^{-1}(v_f - v_s) + \theta v_f - \phi \nabla p \in L^2(0, T; [H^{-1}(\Omega)]^d)$. Since functions in $[H_0^1(\Omega)]^d \otimes D(0, T)$ generate $L^2(0, T; [H_0^1(\Omega)]^d)$, it follows that $\partial_t v_f \in L^2(0, T; [H^{-1}(\Omega)]^d)$ and that, for any test function $w_f \in L^2(0, T; [H_0^1(\Omega)]^d)$,

$$ \int_0^T \left\langle \rho_f \phi \partial_t v_f, w_f \right\rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} \, dt + \int_0^T \int_\Omega \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx \, dt + \int_0^T \int_\Omega \phi^2 k_f^{-1}(v_f - \partial_t u_s) \cdot w_f \, dx \, dt \\ - \int_0^T \int_\Omega \theta v_f \cdot w_f \, dx \, dt - \int_0^T \int_\Omega p \operatorname{div}(\phi w_f) \, dx \, dt = \int_0^T \int_\Omega \rho_f \phi f \cdot w_f \, dx \, dt. \quad (49) $$

Similarly, we infer from (47a) and (47c) that $\partial_{tt}u_s \in L^2(0, T; [H^{-1}(\Omega)]^d)$ and $\frac{\alpha - \phi}{\kappa} \partial_t p + \operatorname{div}((\alpha - \phi)\partial_t u_s) \in L^2(0, T; L^2(\Omega))$, and that for any $w_s \in L^2(0, T; [H_0^1(\Omega)]^d)$

$$ \int_0^T \left\langle \rho_s(1-\phi)\partial_{tt}u_s, w_s \right\rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} \, dt + \int_0^T \int_\Omega \sigma_s(u_s) : \varepsilon(w_s) \, dx \, dt \\ - \int_0^T \int_\Omega \phi^2 k_f^{-1}(v_f - \partial_t u_s) \cdot w_s \, dx \, dt - \int_0^T \int_\Omega p \operatorname{div}((\alpha - \phi)w_s) \, dx \, dt = \int_0^T \int_\Omega \rho_s(1-\phi)f \cdot w_s \, dx \, dt, \quad (50) $$

and for any $q \in L^2(0, T; L^2(\Omega))$

$$ \int_0^T \int_\Omega \left( \frac{\alpha - \phi}{\kappa} \partial_t p + \operatorname{div}((\alpha - \phi)\partial_t u_s) \right) q \, dx \, dt + \int_0^T \int_\Omega \operatorname{div}(\phi v_f) q \, dx \, dt = 0. \quad (51) $$

Note that the main difference compared to the viscous case is that $\partial_t p$ is not in $L^2(0, T; L^2(\Omega))$ any more since the structure velocity $\partial_t u_s$ does not belong to $L^2(0, T; [H_0^1(\Omega)]^d)$. Yet $\partial_t \left( \frac{\alpha - \phi}{\kappa} p + \operatorname{div}((\alpha - \phi)u_s) \right) \in L^2(0, T; L^2(\Omega))$ and $\partial_t p \in L^2(0, T; H^{-1}(\Omega))$.

Let us now prove that the weak solution is unique. Let $(u_s, v_f, p)$ be a solution to (47) with zero initial conditions and source terms, and let $\tau$ be given in $(0, T)$. Contrary to the viscous case, we can not take $w_s = \partial_t u_s$ as test functions in (50) because $\partial_t u_s \notin L^2(0, T; [H_0^1(\Omega)]^d)$. To overcome this lack of regularity, we consider the so-called Ladyzhenskaya tests functions [68]. For the solid and pressure equations, we use the same tests functions that were considered in [12, Theorem 3] and [91, Section 4.2.2.] for Biot's consolidation model, namely

$$ \psi_s(t) = \begin{cases} - \int_t^\tau u_s(\sigma) \, d\sigma & \text{if } \tau \ge t \\ 0 & \text{if } \tau \le t \end{cases} \quad \text{and} \quad \psi_p(t) = \begin{cases} - \int_t^\tau \int_0^v p(\sigma) \, d\sigma \, dv & \text{if } \tau \ge t \\ 0 & \text{if } \tau \le t. \end{cases} $$
To these tests functions, we have to add a fluid test function which is built in the very same manner and corresponds to the fluid counterpart of the previous structure and pressure test functions. Therefore we consider

$$ \psi_f(t) = \begin{cases} - \int_t^\tau \int_0^v v_f(\sigma) \, d\sigma \, dv & \text{if } \tau \ge t \\ 0 & \text{if } \tau \le t. \end{cases} $$

The functions $\psi_s$ and $\psi_f$ belong to $C^1([0, T]; [H_0^1(\Omega)]^d)$, $\psi_p$ belongs to $C^1([0, T]; L^2(\Omega))$ and they are admissible tests functions. For $t \le \tau$, remembering that the considered solution is associated with zero initial conditions, they satisfy

$$ \psi_s(\tau) = 0, \quad \partial_t \psi_s(t) = u_s(t), \quad \partial_t \psi_s(0) = 0. \qquad (52) $$

$$ \psi_p(\tau) = 0, \quad \partial_t \psi_p(t) = \int_0^t p(\sigma) \, d\sigma, \quad \partial_{tt} \psi_p(t) = p(t), \quad \partial_{tt} \psi_p(0) = 0, \qquad (53) $$

and

$$ \psi_f(\tau) = 0, \quad \partial_t \psi_f(t) = \int_0^t v_f(\sigma) \, d\sigma, \quad \partial_{tt} \psi_f(t) = v_f(t), \quad \partial_{tt} \psi_f(0) = 0. \qquad (54) $$

Taking $\psi_s$ as a test function in (50) we compute the different terms. Due to (52), we have in a standard way (see [71] in the case of an abstract second order equation or [91] for the Biot's consolidation model)

$$ \int_0^\tau \langle \rho_s(1-\phi) \partial_{tt} u_s, \psi_s \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} dt = -\frac{1}{2} \int_{\Omega} \rho_s(1-\phi) |u_s(\tau)|^2 dx, $$

$$ \int_0^\tau \int_{\Omega} \sigma_s(u_s) : \varepsilon(\psi_s) dx dt = -\frac{1}{2} \int_{\Omega} \sigma_s(\psi_s(0)) : \varepsilon(\psi_s(0)) dx, $$

Moreover, since $v_f(t) = \partial_{tt} \psi_f(t)$, $\partial_t u_s(t) = \partial_{tt} \psi_s(t)$, $\partial_t \psi_f(0) = \partial_t \psi_s(0) = 0$ and $\psi_s(\tau) = 0$, the friction term writes, after integration by parts in time,

$$ - \int_0^\tau \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot \psi_s \, dx \, dt = \int_0^\tau \int_{\Omega} \phi^2 k_f^{-1} (\partial_t \psi_f - \partial_t \psi_s) \cdot \partial_t \psi_s \, dx \, dt. $$

Finally we obtain the following identity

$$ -\frac{1}{2} \int_{\Omega} \rho_s(1-\phi) |u_s(\tau)|^2 dx - \frac{1}{2} \int_{\Omega} \sigma_s(\psi_s(0)) : \varepsilon(\psi_s(0)) dx \\ + \int_{0}^{\tau} \int_{\Omega} \phi^2 k_f^{-1} (\partial_t \psi_f - \partial_t \psi_s) \cdot \partial_t \psi_s dx dt - \int_{0}^{\tau} \int_{\Omega} p \operatorname{div} ((\alpha - \phi)\psi_s) dx dt = 0. \qquad (55) $$

Let us now focus on the fluid equation. We take $\psi_f$ as a test function in (49). Due to the properties (54), the fluid inertial and viscous terms become respectively

$$ \int_0^\tau \langle \rho_f \phi \partial_t v_f, \psi_f \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} dt = - \int_0^\tau \int_{\Omega} \rho_f \phi v_f \cdot \partial_t \psi_f dx dt = - \frac{1}{2} \int_{\Omega} \rho_f \phi |\partial_t \psi_f(\tau)|^2 dx, $$

and

$$ \int_0^\tau \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(\psi_f) dx dt = - \int_0^\tau \int_{\Omega} \phi \sigma_f(\partial_t \psi_f) : \varepsilon(\partial_t \psi_f) dx dt. $$

Once again the friction term can be transformed as follows

$$ \int_0^\tau \int_{\Omega} \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot \psi_f \, dx \, dt = - \int_0^\tau \int_{\Omega} \phi^2 k_f^{-1} (\partial_t \psi_f - \partial_t \psi_s) \cdot \partial_t \psi_f \, dx \, dt, $$
and, thanks to (54), (53), an integration by parts in time in the pressure term yields

$$ - \int_0^\tau \int_\Omega p \operatorname{div}(\phi \psi_f) \, dx \, dt = \int_0^\tau \int_\Omega \partial_t \psi_p \operatorname{div}(\phi \partial_t \psi_f) \, dx \, dt. $$

The last term, involving $\theta$, writes

$$ - \int_0^\tau \int_\Omega \theta v_f \cdot \psi_f \, dx \, dt = \int_0^\tau \int_\Omega \theta |\partial_t \psi_f|^2 \, dx \, dt + \int_0^\tau \int_\Omega \partial_t \theta \psi_f \cdot \partial_t \psi_f \, dx \, dt. $$

Summing up all these contributions implies

$$ -\frac{1}{2} \int_{\Omega} \rho_f \phi |\partial_t \psi_f(\tau)|^2 \, dx - \int_{0}^{\tau} \int_{\Omega} \phi \sigma_f(\partial_t \psi_f) : \varepsilon(\partial_t \psi_f) \, dx \, dt - \int_{0}^{\tau} \int_{\Omega} \phi^2 k_f^{-1} (\partial_t \psi_f - \partial_t \psi_s) \cdot \partial_t \psi_f \, dx \, dt \\ + \int_{0}^{\tau} \int_{\Omega} \partial_t \psi_p \operatorname{div}(\phi \partial_t \psi_f) \, dx \, dt = - \int_{0}^{\tau} \int_{\Omega} \theta |\partial_t \psi_f|^2 \, dx \, dt - \int_{0}^{\tau} \int_{\Omega} \partial_t \theta \psi_f \cdot \partial_t \psi_f \, dx \, dt. \quad (56) $$

Next we take $\psi_p$ as a test function in (51). As in [91] we have

$$ \int_0^\tau \int_\Omega \left( \frac{\alpha - \phi}{\kappa} \partial_t p + \operatorname{div}((\alpha - \phi)\partial_t u_s) \right) \psi_p \, dx \, dt = -\frac{1}{2} \int_\Omega \frac{\alpha - \phi}{\kappa} |\partial_t \psi_p(\tau)|^2 \, dx + \int_0^\tau \int_\Omega p \operatorname{div}((\alpha - \phi)\psi_s) \, dx \, dt. $$

Moreover

$$ \int_0^\tau \int_\Omega \operatorname{div}(\phi v_f) \psi_p \, dx \, dt = - \int_0^\tau \int_\Omega \operatorname{div}(\phi \partial_t \psi_f) \partial_t \psi_p \, dx \, dt, $$

and thus we obtain the following identity

$$ -\frac{1}{2} \int_{\Omega} \frac{\alpha - \phi}{\kappa} |\partial_t \psi_p(\tau)|^2 \, dx + \int_{0}^{\tau} \int_{\Omega} p \operatorname{div}((\alpha - \phi)\psi_s) \, dx \, dt - \int_{0}^{\tau} \int_{\Omega} \operatorname{div}(\phi \partial_t \psi_f) \partial_t \psi_p \, dx \, dt = 0. \quad (57) $$

Summing up (55), (56) and (57), we obtain

$$ \begin{aligned} & \frac{1}{2} \int_{\Omega} \rho_s(1-\phi) |u_s(\tau)|^2 \, dx + \frac{1}{2} \int_{\Omega} \rho_f \phi |\partial_t \psi_f(\tau)|^2 \, dx + \frac{1}{2} \int_{\Omega} \frac{\alpha-\phi}{\kappa} |\partial_t \psi_p(\tau)|^2 \, dx \\ & + \frac{1}{2} \int_{\Omega} \sigma_s(\psi_s(0)) : \varepsilon(\psi_s(0)) \, dx + \int_{0}^{\tau} \int_{\Omega} \phi \sigma_f(\partial_t \psi_f) : \varepsilon(\partial_t \psi_f) \, dx \, dt + \int_{0}^{\tau} \int_{\Omega} \phi^2 k_f^{-1} (\partial_t \psi_f - \partial_t \psi_s)^2 \, dx \, dt \\ & = \int_{0}^{\tau} \int_{\Omega} \theta |\partial_t \psi_f|^2 \, dx \, dt + \int_{0}^{\tau} \int_{\Omega} \partial_t \theta \psi_f \cdot \partial_t \psi_f \, dx \, dt. \end{aligned} \quad (58) $$

To conclude, we observe that

$$ \int_0^\tau \int_\Omega \theta |\partial_t \psi_f|^2 \, dx \, dt \leq (\rho_f \phi_{\min})^{-1} \|\theta\|_{L^\infty((0,T)\times\Omega)} \int_0^\tau \int_\Omega \rho_f \phi |\partial_t \psi_f|^2 \, dx \, dt $$

and that, since $\psi_f(t) = - \int_t^\tau \partial_t \psi_f(\sigma) d\sigma$, we can estimate the last term of (58) as follows

$$ \begin{aligned} \int_0^\tau \int_\Omega \partial_t \theta \psi_f \cdot \partial_t \psi_f \, dx \, dt &\le \|\partial_t \theta\|_{L^\infty((0,T)\times\Omega)} \int_0^\tau \int_\Omega \left| \int_t^\tau \partial_t \psi_f(\sigma) \, d\sigma \right| |\partial_t \psi_f(t)| \, dx \, dt \\ &\le \|\partial_t \theta\|_{L^\infty((0,T)\times\Omega)} \int_\Omega \int_0^\tau \int_0^\tau |\partial_t \psi_f(\sigma)| |\partial_t \psi_f(t)| \, d\sigma \, dt \, dx \\ &\le T (\rho_f \phi_{\min})^{-1} \|\partial_t \theta\|_{L^\infty((0,T)\times\Omega)} \int_0^\tau \int_\Omega \rho_f \phi |\partial_t \psi_f|^2 \, dx \, dt. \end{aligned} $$

Consequently, using Grönwall Lemma, we deduce that $u_s = \partial_t \psi_f = \partial_t \psi_p = 0$. Hence $u_s = v_f = p = 0$, which proves the uniqueness of the variational solution.
Now that we know that the variational solution is unique, it follows that it is necessarily the one obtained by the approximation process built from $(VF)^n$. Since this approximation process is based on the energy estimate (39) with $\eta = 0$, we can pass to the limit in this estimation to get (11).

In particular, to show that the mild solution is equal to the variational solution, it is sufficient to prove that it also derives from this approximation process. Let us denote by $z$ the mild solution given by Theorem 2.11, and remind the notation $G^n(t)(u_s, v_s, v_f, p) = (0,0, -(\rho_f\phi)^{-1}\theta^n(t)v_f)$, with $\theta^n \in C^1([0, T]; L^\infty(\Omega))$ converging towards $\theta$ in $C^0([0, T] \times \Omega)$. From Duhamel formula, it holds that

$$z(t) - z^n(t) = \Phi_0^\kappa(t)(z_0 - z_0^n) + \int_0^t \Phi_0^\kappa(t-s)(g(s) - g^n(s)) \, ds - \int_0^t \Phi_0^\kappa(t-s)(G(s)z(s) - G^n(s)z(s)) \, ds$$

Writing $G(s)z(s) - G^n(s)z^n(s) = (G(s) - G^n(s))z^n(s) + G(s)(z(s) - z^n(s))$ and recalling that $\Phi_0^\kappa$ is a $C^0$-semigroup of contraction, we infer

$$\|z(t) - z^n(t)\|_Z \le \|z_0 - z_0^n\|_Z + \int_0^t \|g(s) - g^n(s)\|_Z \, ds \\ + (\rho_f \phi_{\min})^{-1} \|\theta - \theta^n\|_{C^0([0,T] \times \Omega)} \int_0^t \|z^n(s)\|_Z \, ds + \omega \int_0^t \|z(s) - z^n(s)\|_Z \, ds,$$

where $\omega$ is defined in (35). Thus, for any $\delta > 0$, we can find $n$ large enough such that

$$\|z(t) - z^n(t)\|_Z \le \delta + \omega \int_0^t \|z(s) - z^n(s)\|_Z \, ds.$$

Using Grönwall Lemma, we conclude that $\|z(t) - z^n(t)\|_Z \le \delta e^{\omega T}$, and hence $z^n \to z$ in $C^0([0, T]; Z)$. $\square$

**Remark 2.14.** In the previous proof, we took advantage of the semigroup framework to show the existence of the variational solution. Note that it could also be shown by regularization of the viscous case, see [12, Theorem 2] where such a regularization is performed on Biot's consolidation model.

**Remark 2.15.** As in the viscous case, we could also define the variational solution without assuming that it is continuous with respect to time, but rather by seeking for

$$\begin{align*}
u_s &\in L^\infty(0, T; [H_0^1(\Omega)]^d), & \partial_t u_s &\in L^\infty(0, T; [L^2(\Omega)]^d) & \text{and} & \partial_{tt} u_s &\in L^2(0, T; [H^{-1}(\Omega)]^d), \\
v_f &\in L^2(0, T; [H_0^1(\Omega)]^d) \cap L^\infty(0, T; [L^2(\Omega)]^d) & \partial_t v_f &\in L^2(0, T; [H^{-1}(\Omega)]^d), \\
p &\in L^\infty(0, T; L^2(\Omega)) & \partial_t \left( \frac{\alpha - \phi}{\kappa} p + \operatorname{div}((\alpha - \phi)u_s) \right) &\in L^2(0, T; L^2(\Omega)),
\end{align*}$$

such that (46) and (47) are verified. With this definition, the continuity in time of the solution can then be recovered using, for instance, a parabolic regularization, while it is obtained directly in the above proof.

In the next section, we analyze the poromechanics problem for an incompressible elastic skeleton, modeled by the assumption $\kappa = +\infty$. This assumption is crucial for targeting biomedical applications since the tissues in our body are mostly composed of water, and thus are close to being incompressible.

# 3 Existence of solutions for an incompressible skeleton $\kappa = +\infty$

When $\kappa = +\infty$ – and thus $\alpha = 1$, see Remark A.1 – the system of equations (18) reads

$$\partial_t u_s - v_s = 0, \tag{59a}$$

$$\rho_s(1-\phi)\partial_t v_s - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) - \operatorname{div}(\sigma_s(u_s)) - \phi^2 k_f^{-1}(v_f - v_s) + (1-\phi)\nabla p = \rho_s(1-\phi)f, \tag{59b}$$

$$\rho_f \phi \partial_t v_f - \operatorname{div}(\phi \sigma_f(v_f)) + \phi^2 k_f^{-1}(v_f - v_s) - \theta v_f + \phi \nabla p = \rho_f \phi f, \tag{59c}$$

$$\operatorname{div}((1-\phi)v_s + \phi v_f) = \frac{\theta}{\rho_f}. \tag{59d}$$
It has to be completed with boundary conditions (3) and initial conditions (4). Note that, in the present case, there is no initial condition for the pressure anymore.

Equation (59d) traduces the mixture's incompressibility, which comes from the assumption that the solid and the fluid phases are both incompressible. It takes the form of a constraint on the divergence of the mixture's velocity.

But, as already noticed in Section 1, it is sufficient to consider the case

$$ \mathrm{div}((1-\phi)v_s + \phi v_f) = 0, \quad \text{in } \Omega \times (0, T). \qquad (60) $$

Indeed, assuming for instance that $\theta \in H^1(0, T; L^2(\Omega)) \cap L^2(0, T; H^1(\Omega))$ and that the compatibility condition $\int_{\Omega} \theta(t) dx = 0$ is satisfied for all $t \in [0, T]$, we consider the Bogovskii's operator [26] and we build $v_{\theta} \in H^1(0, T; [H_0^1(\Omega)]^d) \cap L^2(0, T; [H^2(\Omega)]^d)$ such that

$$ \mathrm{div} v_{\theta} = \frac{\theta}{\rho_f}. \qquad (61) $$

The change of variables $\hat{v}_s = v_s - v_\theta$ and $\hat{v}_f = v_f - v_\theta$ gives $\mathrm{div}((1-\phi)\hat{v}_s + \phi\hat{v}_f) = 0$ by construction. Furthermore, $(u_s - \int_0^t v_\theta(s) ds, \hat{v}_s, \hat{v}_f, p)$ verifies (59a), (59b) and (59c) with right-hand sides that are different but still regular since $v_\theta \in H^1(0, T; [H_0^1(\Omega)]^d) \cap L^2(0, T; [H^2(\Omega)]^d)$.

The first part of this section is devoted to the functional analysis of the coupling constraint (60).

## 3.1 Functional framework

We consider the space

$$ V_{\phi} = \{ (v_s, v_f) \in [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d : \mathrm{div}((1-\phi)v_s + \phi v_f) = 0 \quad \text{in } \Omega \} $$

of functions in $[H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d$ satisfying the constraint (60). Let us also define the space $H_\phi$ as the closure of $V_\phi$ in $[L^2(\Omega)]^d \times [L^2(\Omega)]^d$.

Then, we introduce the mixture's divergence operator defined by

$$ B: \begin{array}{rcl} [H_0^1(\Omega)]^d & \times & [H_0^1(\Omega)]^d \\ (v_s, v_f) & \longrightarrow & \mathrm{div}((1-\phi)v_s + \phi v_f). \end{array} $$

The bounded operator $B$ satisfies the following inf-sup condition.

**Proposition 3.1.** Assume that $\phi \in H^{d/2+r}(\Omega)$ with $r > 0$. There exists $\underline{\beta} > 0$ such that, for all $p \in L_0^2(\Omega)$,

$$ \sup_{(v_s, v_f) \in [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d} \frac{\int_{\Omega} \mathrm{div}((1-\phi)v_s + \phi v_f) p \, dx}{\|(v_s, v_f)\|_{[H_0^1(\Omega)]^d} \times \| [H_0^1(\Omega)]^d \|} \geq \underline{\beta} \|p\|. \qquad (62) $$

*Proof.* There exists $C_{\mathrm{div}} > 0$ such that for any $p \in L_0^2(\Omega)$, there exists $v_p \in [H_0^1(\Omega)]^d$ satisfying

$$ \mathrm{div} v_p = p \quad \text{and} \quad \|\nabla v_p\| \leq C_{\mathrm{div}} \|p\|. \qquad (63) $$

Setting $v = (v_p, v_p)$, we have $Bv = p$ by construction and $\|v\|_{[H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d} \leq C\|p\|$ from the above inequality. $\square$

**Remark 3.2.** Note that the constant $\underline{\beta}$ of the above inf-sup condition does not depend on the porosity $\phi$ and is although valid for $\phi = 0$ (resp. $\bar{\phi} = 1$) which are the limit cases for which there is no fluid (resp. no structure).
This inf-sup condition allows us to state the following result, which is a generalization of De Rham Theorem [100, 63, 32]. It is a key ingredient to show the existence of pressure in the incompressible case.

**Theorem 3.3.** Assume that $\phi \in H^{d/2+r}(\Omega)$ with $r > 0$. If $f = (f_s, f_f) \in [H^{-1}(\Omega)]^d \times [H^{-1}(\Omega)]^d$ satisfies

$$ \langle f, w \rangle = \langle f_s, w_s \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} + \langle f_f, w_f \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} = 0, \quad \forall w = (w_s, w_f) \in V_\phi, $$

then there exists a unique $p \in L_0^2(\Omega)$ such that $f_s = -(1-\phi)\nabla p$ and $f_f = -\phi \nabla p$.

*Proof.* The proof follows standard arguments and is based on the Closed Range Theorem. Note that $V_\phi = \text{Ker } B$. Let us characterize the adjoint of $B$. Since $\phi \in H^{d/2+r}(\Omega)$ with $r > 0$, $\phi$ is a multiplier of $[H_0^1(\Omega)]^d$, namely

$$ \forall v \in [H_0^1(\Omega)]^d, \quad \phi v \in [H_0^1(\Omega)]^d \quad \text{and} \quad \|\nabla(\phi v)\| \le C_\phi \|\nabla v\| $$

for some positive constant $C_\phi$. Therefore, for all $p \in L_0^2(\Omega)$, $(1-\phi)\nabla p$ and $\phi \nabla p$ belong to $[H^{-1}(\Omega)]^d$ so that we can define the adjoint operator as

$$ B^* : \begin{array}{rcl} L_0^2(\Omega) & \to & [H^{-1}(\Omega)]^d \times [H^{-1}(\Omega)]^d \\ p & \mapsto & -(1-\phi) \nabla p, -\phi \nabla p. \end{array} $$

Thanks to Proposition 3.1, the Closed Range Theorem implies that $(\text{Ker } B)^\circ = \text{Rg } B^*$. Consequently, for any $f \in (\text{Ker } B)^\circ = (V_\phi)^\circ$, namely for any $f = (f_s, f_f) \in [H^{-1}(\Omega)]^d \times [H^{-1}(\Omega)]^d$ satisfying

$$ \langle f, w \rangle = 0, \quad \forall w \in V_\phi, $$

there exists a unique $p \in L_0^2(\Omega)$ such that $f_s = -(1-\phi)\nabla p$ and $f_f = -\phi \nabla p$.

<figure><img src="image_4.png" /></figure>

Theorem 3.3 allows us to characterize the space $H_\phi$ in the following way.

**Proposition 3.4.** The space $H_\phi$ can be expressed as

$$ H_{\phi} = \{ (v_s, v_f) \in [L^2(\Omega)]^d \times [L^2(\Omega)]^d : \begin{array}{l} \text{div}((1-\phi)v_s + \phi v_f) = 0 \quad \text{in } \mathcal{D}'(\Omega) \\ \text{and} \quad ((1-\phi)v_s + \phi v_f) \cdot n = 0 \quad \text{on } \partial\Omega \end{array} \}. $$

*Proof.* We denote by $\mathcal{H}$ the space

$$ \{ (v_s, v_f) \in [L^2(\Omega)]^d \times [L^2(\Omega)]^d : \begin{array}{l} \text{div}((1-\phi)v_s + \phi v_f) = 0 \quad \text{in } \mathcal{D}'(\Omega) \\ \text{and} \quad ((1-\phi)v_s + \phi v_f) \cdot n = 0 \quad \text{on } \partial\Omega \end{array} \}. $$

Let $v = (v_s, v_f)$ be an element of $H_\phi$. By definition, $H_\phi$ is the closure of $V_\phi$ in $[L^2(\Omega)]^d \times [L^2(\Omega)]^d$, so there exists a sequence $(v_s^n, v_f^n)$ belonging to $V_\phi$ that converges towards $v$ in $[L^2(\Omega)]^d \times [L^2(\Omega)]^d$. Since $\text{div}((1-\phi)v_s^n + \phi v_f^n) = 0$, the equality $\text{div}((1-\phi)v_s + \phi v_f) = 0$ holds true in the limit. Further, $(1-\phi)v_s + \phi v_f \in \text{H}_{\text{div}}(\Omega)$. The continuity of the normal trace operator then implies that

$$ \| ((1-\phi)v_s + \phi v_f) \cdot n - ((1-\phi)v_s^n + \phi v_f^n) \cdot n \|_{H^{-1/2}(\partial\Omega)} \xrightarrow{n \to \infty} 0, $$

which implies that $((1-\phi)v_s + \phi v_f) \cdot n = 0$, since $(1-\phi)v_s^n + \phi v_f^n) \cdot n = 0$. Hence, $H_\phi \subset \mathcal{H}$.

Now, let us prove the other embedding. Let denote by $\mathcal{H}^*$ the orthogonal complement of $H_\phi$ into $\mathcal{H}$ and let $f = (f_s, f_f)$ be an element of $\mathcal{H}^*$. Noting that $\mathcal{H}^* \subset H_\phi^\perp$ and $V_\phi \subset H_\phi$, it follows from Theorem 3.3 that there exists a pressure $p \in L_0^2(\Omega)$ such that $f_s = -(1-\phi)\nabla p$ and $f_f = -\phi \nabla p$. Moreover since $\nabla p = -(f_s + f_f)$ and $(f_s, f_f) \in [L^2(\Omega)]^d \times [L^2(\Omega)]^d$, we get $p \in H^1(\Omega)$. Since $f$ belongs to $\mathcal{H}$, we have $\text{div}((1-\phi)^2 + \phi^2)\nabla p) = \text{div}((1-\phi)^2 \nabla p + \phi^2 \nabla p) = 0$ in $\mathcal{D}'(\Omega)$ and $((1-\phi)^2 \nabla p + \phi^2 \nabla p) \cdot n = 0$. Thus $p$ is equal to zero (up to a constant) as the unique solution of an elliptic Neumann problem, so $\nabla p = 0$ and $f = 0$. In conclusion, $\mathcal{H}^* = \{0\}$, which proves that $H_\phi = \mathcal{H}$.

<figure><img src="image_5.png" /></figure>
**Remark 3.5.** If $\phi \in C^\infty(\Omega)$, one can show that $V_\phi$ and $H_\phi$ correspond to the closures of the space $V_\phi$ in $[H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d$ and $[L^2(\Omega)]^d \times [L^2(\Omega)]^d$, where

$$V_\phi = \{(v_s, v_f) \in [D(\Omega)]^d \times [D(\Omega)]^d : \operatorname{div}((1-\phi)v_s + \phi v_f) = 0\}.$$

We are now going to combine this functional framework adapted to the constraint (60) with the semigroup approach in order to study Problem (59). Here again, we investigate the cases $\eta > 0$ and $\eta = 0$ separately. To simplify the proof we consider that $\theta$, which appears now only in the term $-\theta v_f$ in the fluid equation, does not depend on time: $\theta \in L^\infty(\Omega)$. It simplifies the proof, but it could be easily modified to include the time-dependent case as in the proofs of Section 2.

## 3.2 The case $\eta > 0$

We formulate the problem in the functional framework established previously. We seek for a solution $z = (u_s, v_s, v_f)$ in the energy space $H = [H_0^1(\Omega)]^d \times H_\phi$ endowed with the scalar product

$$(z, y)_H = \int_{\Omega} \sigma_s(u_s) : \varepsilon(d_s) + \int_{\Omega} \rho_s(1-\phi)v_s \cdot w_s \, dx + \int_{\Omega} \rho_f \phi v_f \cdot w_f \, dx,$$

for any $z = (u_s, v_s, v_f)$, $y = (d_s, w_s, w_f)$ belonging to $H$, and with the corresponding norm

$$\|z\|_H^2 = \|u_s\|_s^2 + \int_{\Omega} \rho_s(1-\phi) |v_s|^2 \, dx + \int_{\Omega} \rho_f \phi |v_f|^2 \, dx.$$

Setting $V = [H_0^1(\Omega)]^d \times V_\phi$, we consider the bilinear form

$$a_{\eta}^{\infty}(z, y) = - \int_{\Omega} \sigma_s(v_s) : \varepsilon(d_s) \, dx + \int_{\Omega} \sigma_s(u_s) : \varepsilon(w_s) \, dx + 2\eta \int_{\Omega} \varepsilon(v_s) : \varepsilon(w_s) \, dx \\ + \int_{\Omega} \phi^2 k_f^{-1}(v_f - v_s) \cdot (w_f - w_s) \, dx + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx - \int_{\Omega} \theta v_f \cdot w_f \, dx \quad (64)$$

defined for all $z = (u_s, v_s, v_f)$ and $y = (d_s, w_s, w_f)$ in $V$. This bilinear form is the same as (22) but without the terms involving the pressure because of the test functions in $V$. Note that here, since we have assumed that $\theta$ does not depend on time, we can include the term $\int_{\Omega} \theta v_f \cdot w_f$ in the definition of the bilinear form associated with our coupled problem. When $\theta$ depends on time, one can not and we have to introduce the operator $G(t)$ which is a bounded perturbation, see (25). As in (23) and (24), we define the unbounded operator $(A_\eta^\infty, D(A_\eta^\infty))$ associated with the bilinear form (64) by

$$(A_{\eta}^{\infty} z, y)_H = a_{\eta}^{\infty}(z, y), \quad \forall z \in D(A_{\eta}^{\infty}), \forall y \in V,$$

in the domain

$$D(A_{\eta}^{\infty}) = \{z \in V : \exists g \in H, a_{\eta}^{\infty}(z, y) = (g, y)_H, \ y \in V\}.$$

The above definitions are quite abstract, in particular because they rely on test functions in the constrained space $V_\phi$. In the next proposition, we recover a more explicit expression of $A_\eta^\infty$ and $D(A_\eta^\infty)$ thanks to the generalization of De Rham's Theorem established previously.

**Proposition 3.6.** The operator's domain can be characterized as

$$D(A_{\eta}^{\infty}) = \left\{ u_s, v_s, v_f \in [H_0^1(\Omega)]^d \left| \begin{array}{l} \exists! p \in L_0^2(\Omega) \text{ such that } \\ -\operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) + (1-\phi)\nabla p \in [L^2(\Omega)]^d \\ -\operatorname{div}(\phi\sigma_f(v_f)) + \phi\nabla p \in [L^2(\Omega)]^d \\ \operatorname{div}((1-\phi)v_s + \phi v_f) = 0, \end{array} \right. \right\}. \quad (65)$$
In addition, for all $z = (u_s, v_s, v_f) \in D(A_{\eta}^{\infty})$ and $g = (g_u, g_s, g_f) \in H$, we have

$$A_{\eta}^{\infty} z = g \Leftrightarrow \exists! p \in L_{0}^{2}(\Omega), \quad \begin{cases} g_u = -v_s, \\ \rho_s(1-\phi)g_s = -\operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) - \phi^2 k_f^{-1}(v_f - v_s) + (1-\phi)\nabla p, \\ \rho_f\phi g_f = -\operatorname{div}(\phi\sigma_f(v_f)) + \phi^2 k_f^{-1}(v_f - v_s) - \theta v_f + \phi \nabla p. \end{cases} \tag{66}$$

*Proof.* Let $z = (u_s, v_s, v_f)$ be an element of $D(A_{\eta}^{\infty})$. By definition, there exists $g = (g_u, g_s, g_f) \in H$ such that $a_{\eta}^{\infty}(z, y) = (g, y)_{H}$ for all $y = (d_s, w_s, w_f) \in \bar{V}$, namely

$$\int_{\Omega} \sigma_s(g_u) : \varepsilon(d_s) \, dx = - \int_{\Omega} \sigma_s(v_s) : \varepsilon(d_s) \, dx, \quad \forall d_s \in [H_0^1(\Omega)]^d, \tag{67}$$

and

$$\int_{\Omega} \rho_s(1-\phi) g_s \cdot w_s \, dx + \int_{\Omega} \rho_f \phi g_f \cdot w_f \, dx = \int_{\Omega} \sigma_s(u_s) : \varepsilon(w_s) \, dx + 2\eta \int_{\Omega} \varepsilon(v_s) : \varepsilon(w_s) \, dx \\ + \int_{\Omega} \phi^2 k_f^{-1}(v_f - v_s) \cdot (w_f - w_s) \, dx + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx - \int_{\Omega} \theta v_f \cdot w_f \, dx, \quad \forall (w_s, w_f) \in V_{\phi}. \tag{68}$$

The relation (67) implies that $g_u = -v_s$ in $[H_0^1(\Omega)]^d$. From (68), we deduce that

$$((g_s, g_f), (w_s, w_f))_{H_{\phi}} = \langle -\operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) - \phi^2 k_f^{-1}(v_f - v_s), w_s \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d} \\ + \langle -\operatorname{div}(\phi \sigma_f(v_f)) + \phi^2 k_f^{-1}(v_f - v_s) - \theta v_f, w_f \rangle_{[H^{-1}(\Omega)]^d, [H_0^1(\Omega)]^d}, \quad \forall (w_s, w_f) \in V_{\phi}.$$

Applying Theorem 3.3, we get the existence of a pressure $p \in L_0^2(\Omega)$ such that

$$\begin{cases} \rho_s(1-\phi)g_s = -\operatorname{div}(\sigma_s(u_s)) - \operatorname{div}(\sigma_s^{\text{vis}}(v_s)) - \phi^2 k_f^{-1}(v_f - v_s) + (1-\phi)\nabla p \quad \text{in } [H^{-1}(\Omega)]^d, \\ \rho_f\phi g_f = -\operatorname{div}(\phi\sigma_f(v_f)) + \phi^2 k_f^{-1}(v_f - v_s) - \theta v_f + \phi \nabla p \quad \text{in } [H^{-1}(\Omega)]^d, \end{cases}$$

which proves (66). Since $g_s$ and $g_f$ belong to $[L^2(\Omega)]^d$, the above relation essentially holds true in $[L^2(\Omega)]^d$, which yields (65). □

**Remark 3.7.** The characterization of the Lagrange multiplier $p$ associated with the constraint on the mixture velocity as the weak solution of an elliptic problem, as done in [8] in the context of fluid-structure interaction problems, is not straightforward precisely because the constraint involves the mixture velocity which is not a natural unknown of our coupled problem.

Lastly, we set $g = (0, f, f)$ and we denote by $\Pi$ the Leray projection operator from $[H_0^1(\Omega)]^d \times [L^2(\Omega)]^d \times [L^2(\Omega)]^d$ into $H$. We are now ready to state the following existence and uniqueness result.

**Theorem 3.8.** Assume that (h1) – (h3) hold true, that $\eta > 0$ and that $\int_{\Omega} \theta \, dx = 0$.

(i) If $z_0 \in D(A_{\eta}^{\infty})$ and $f \in H^1(0, T; [L^2(\Omega)]^d)$ so that $\Pi g \in H^1(0, T; H)$, then there exists a unique strong solution $z \in C^1([0, T]; H) \cap C^0([0, T]; D(A_{\eta}^{\infty}))$ satisfying

$$\left\{ \begin{array}{ll} \dot{z}(t) + A_{\eta}^{\infty} z(t) = \Pi g(t), & t \in [0, T], \\ z(0) = z_0. & \end{array} \right. \tag{69}$$
(ii) If $z_0 \in H$ and $f \in L^2(0, T; [L^2(\Omega)]^d)$ so that $\Pi g \in L^2(0, T; H)$, then Problem (69) has a unique mild solution $z \in C^0([0, T]; H)$ such that $z(0) = z_0$ and

$$ \int_{0}^{T} z(t)\psi(t) \, dt \in D(A_{\eta}^{\infty}), \\ -\int_{0}^{T} z(t)\dot{\psi}(t) \, dt + A_{\eta}^{\infty}\left(\int_{0}^{T} z(t)\psi(t) \, dt\right) = \int_{0}^{T} \Pi g(t)\psi(t) \, dt, $$

for all $\psi \in C_c^1([0, T]; \mathbb{R})$. Moreover, $z$ is given by the Duhamel formula

$$ z(t) = \Phi_{\eta}^{\infty}(t)z_0 + \int_{0}^{t} \Phi_{\eta}^{\infty}(t-s)\Pi g(s) \, ds, \quad (70) $$

where $\Phi_{\eta}^{\infty}$ denotes the continuous semigroup generated by $A_{\eta}^{\infty}$ in the sense that

$$ A_{\eta}^{\infty} x = -\frac{d}{dt} (\Phi_{\eta}^{\infty}(t)x) |_{t=0^{+}}, \quad x \in H. \quad (71) $$

*Proof.* The proof of this result is almost similar to the proof of Theorem 2.2, replacing $Z$ by $H$ and $Y$ by $V$. The only difference is that the term $-\theta v_f$ is treated within the operator $A_{\eta}^{\infty}$ instead of being considered as a perturbation.

For any $z = (u_s, v_s, v_f) \in D(A_{\eta}^{\infty})$, we observe that

$$ a_{\eta}^{\infty}(z, z) = 2\eta \int_{\Omega} |\varepsilon(v_s)|^2 \, dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f - v_s) \cdot (v_f - v_s) \, dx + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(v_f) \, dx - \int_{\Omega} \theta |v_f|^2 \, dx. $$

Thus $(A_{\eta}^{\infty}z, z)_H \geq -\omega\|z\|_H^2$ with $\omega = (\rho_f\phi_{\min})^{-1}\|\theta\|_{L^{\infty}(\Omega)}$.

Moreover, for all $\lambda_0 > \omega$, the operator $A_{\eta}^{\infty} + \lambda_0 I$ is surjective from $D(A_{\eta}^{\infty})$ to $H$ because

$$ a_{\eta}^{\infty}(z, z) + \lambda_0(z, z)_H \geq \lambda_0 \|u_s\|_s^2 + 2\eta \|\varepsilon(v_s)\|^2 + 2\mu_f \phi_{\min} \|\varepsilon(v_f)\|^2 + (\lambda_0 \rho_f \phi_{\min} - \|\theta\|_{L^{\infty}(\Omega)}) \|v_f\|^2. $$

From Lumer-Phillips theorem, we deduce that $A_{\eta}^{\infty}$ is the generator – in the sense of (71) – of a strongly continuous semigroup and the conclusion follows from [14, Part II, Chapter 1, Propositions 3.1–3.3] and [36, Corollary 2.25]. □

**Remark 3.9.** By reproducing the proof of Theorem 2.2, we can extend the result of Theorem 3.8 for a time-dependent $\theta$ satisfying the compatibility condition $\int_{\Omega} \theta(t) \, dx = 0$ for all $t \in [0, T]$. The existence of a mild solution then requires that $\theta \in C^0([0, T] \times \Omega)$ and the existence of a strong solution is guaranteed under the assumption $\theta \in C^1([0, T]; L^{\infty}(\Omega))$. However, because of the lifting (61), more regularity on $\theta$ is needed for the original problem to be well-posed. More precisely, when performing the change of variables $(u_s, v_f, p) \mapsto (u_s - \int_0^t v_{\theta}(s) \, ds, v_f - v_{\theta}, p)$, the right-hand sides of (59b) and (59c) become respectively

$$ \rho_s(1-\phi)f - \rho_s(1-\phi)\partial_t v_{\theta} + \text{div} \left( \sigma_s \left( \int_0^t v_{\theta}(s) \, ds \right) \right) + \text{div} \left( \sigma_s^{\text{vis}}(v_{\theta}) \right) \quad \text{and} \quad \rho_f \phi f - \rho_f \phi \partial_t v_{\theta} + \text{div} \left( \phi \sigma_f(v_{\theta}) \right) + \theta v_{\theta}. $$

The existence of a mild solution requires that all these terms belong to $L^2(0, T; [L^2(\Omega)]^d)$ and thus that the lifting $v_{\theta} \in H^1(0, T; [L^2(\Omega)]^d) \cap L^2(0, T; [H^2(\Omega)]^d)$, which is ensured if $\theta \in H^1(0, T; L^2(\Omega)) \cap L^2(0, T; H^1(\Omega))$. Similarly, the existence of a strong solution is guaranteed under the assumption $\theta \in H^2(0, T; L^2(\Omega)) \cap H^1(0, T; H^1(\Omega))$.

The previous theorem only involves the displacement and velocity fields. The existence of pressure and the relation between (69) and the original system (59) are precised below.
**Corollary 3.10.** Assume that $z_0 \in D(A_\eta^\infty)$, $f \in H^1(0, T; [L^2(\Omega)]^d)$ and let $z = (u_s, v_s, v_f) \in C^1([0, T]; H) \cap C^0([0, T]; D(A_\eta^\infty))$ be the strong solution of (69). There exists a unique pressure $p \in C^0([0, T]; L_0^2(\Omega))$ such that $(z, p)$ satisfies (59) pointwise almost everywhere.

*Proof.* Let $z = (u_s, v_s, v_f)$ be the solution of (69). Since $z \in C^1([0, T]; H) \cap C^0([0, T]; D(A_\eta^\infty))$, for almost every $t \in (0, T)$, the equation

$$A_\eta^\infty z(t) = \Pi g(t) - \dot{z}(t)$$

holds true in the energy space $H$, where we recall that $g = (0, f, f)$. Thus for almost every $t \in (0, T)$, Proposition 3.6 ensures the existence of a pressure $p(t) \in L_0^2(\Omega)$ such that

$$\begin{cases} -\partial_t u_s = -v_s, \\ \rho_s(1-\phi)f - \rho_s(1-\phi)\partial_t v_s = -\operatorname{div}(\sigma_s^{\text{vis}}(v_s)) - \operatorname{div}(\sigma_s(u_s)) - \phi^2 k_f^{-1}(v_f - v_s) + (1-\phi)\nabla p, \\ \rho_f \phi f - \rho_f \phi \partial_t v_f = -\operatorname{div}(\phi \sigma_f(v_f)) + \phi^2 k_f^{-1}(v_f - v_s) + \phi \nabla p. \end{cases}$$

In virtue of (65), the two last lines of the above system are verified at least in $[L^2(\Omega)]^d$ (the first one is satisfied in $[H_0^1(\Omega)]^d$). Hence (59) is satisfied almost everywhere. $\square$

If the input data are less regular, we get the existence of displacement, velocities and pressure in the following weak sense.

**Theorem 3.11.** Assume that (h1) – (h4) are satisfied, $\eta > 0$, $\int_\Omega \theta \, dx = 0$ and $z_0 = (u_{s0}, v_{s0}, v_{f0}) \in H$. Then there exists a unique variational solution $u_s \in C^0([0, T]; [H_0^1(\Omega)]^d)$, $\partial_t u_s \in C^0([0, T]; [L^2(\Omega)]^d)$ and $v_f \in C^0([0, T]; [L^2(\Omega)]^d)$ with $(\partial_t u_s, v_f) \in L^2(0, T; V_\phi)$ such that

$$(u_s(0), \partial_t u_s(0), v_f(0)) = (u_{s0}, v_{s0}, v_{f0})$$

and, for all $(w_s, w_f) \in V_\phi$, the following equation holds in $D'(0, T)$:

$$\begin{aligned} & \frac{d^2}{dt^2} \int_{\Omega} \rho_s(1-\phi) u_s(t) \cdot w_s \, dx + \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(w_s) \, dx + 2\eta \int_{\Omega} \varepsilon(\partial_t u_s(t)) : \varepsilon(w_s) \, dx \\ & + \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \, dx + \int_{\Omega} \phi \sigma_f(v_f(t)) : \varepsilon(w_f) \, dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - \partial_t u_s(t)) \cdot (w_f - w_s) \, dx \\ & - \int_{\Omega} \theta v_f(t) \cdot w_f \, dx = \int_{\Omega} \rho_s(1-\phi)f(t) \cdot w_s \, dx + \int_{\Omega} \rho_f \phi f(t) \cdot w_f \, dx. \end{aligned} \quad (72)$$

The energy estimate (17) holds true and the variational solution coincides with the mild solution given by (70). Furthermore, there exists a unique pressure $p$ such that $(u_s, v_s, v_f, p)$ satisfies (59) in the distribution sense, with $v_s = \partial_t u_s$.

**Remark 3.12.** The incompressibility constraint is satisfied since $(\partial_t u_s, v_f)$ belongs to $L^2(0, T; V_\phi)$. It can be written in variational form as

$$\forall q \in L^2(\Omega), \quad \int_{\Omega} \operatorname{div}((1-\phi)\partial_t u_s(t) + \phi v_f(t)) q \, dx = 0.$$

*Proof.* The proof of existence of $(u_s, v_f)$ follows exactly the same lines as in the compressible case. We build a sequence of strong solutions for smooth data. These solutions $(u_s^n, v_s^n, v_f^n)_n$ satisfy the energy estimate (17) and constitute a Cauchy sequence in $C^0([0, T]; H)$. Moreover $(v_s^n, v_f^n)$ is a Cauchy sequence in $L^2(0, T; V_\phi)$
and we can pass to the limit as we did in the proof of Theorem 2.5. We get the first order system:

$$ \left\{ \begin{array}{l} \forall t \in [0, T], \forall (d_s, w_s, w_f) \in V = [H_0^1(\Omega)]^d \times V_\phi, \\ \frac{d}{dt} \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(d_s) \, dx - \int_{\Omega} \sigma_s(v_s(t)) : \varepsilon(d_s) \, dx = 0, \end{array} \right. \qquad (73a) $$

$$ \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} \int_{\Omega} \rho_s(1-\phi)v_s(t) \cdot w_s \, dx + \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(w_s) \, dx + 2\eta \int_{\Omega} \varepsilon(v_s(t)) : \varepsilon(w_s) \, dx \\ \displaystyle \qquad + \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \, dx + \int_{\Omega} \phi \sigma_f(v_f(t)) : \varepsilon(w_f) \, dx + \int_{\Omega} \phi^2 k_f^{-1}(v_f(t) - v_s(t)) \cdot (w_f - w_s) \, dx \\ \displaystyle \qquad - \int_{\Omega} \theta v_f(t) \cdot w_f \, dx = \int_{\Omega} \rho_s(1-\phi)f(t) \cdot w_s \, dx + \int_{\Omega} \rho_f \phi f(t) \cdot w_f \, dx, \end{array} \right. \qquad (73b) $$

which can be rewritten in second order to obtain (72).

Apart for the regularity provided by the energy estimate (in particular $(v_s, v_f) \in L^2(0, T; V_\phi)$), like in the compressible regime, the previous system provides some regularity on the time derivatives of the solution. The first equation (73a) states that $\partial_t u_s = v_s$ in $L^2(0, T; [H_0^1(\Omega)]^d)$ and the second equation (73b) implies that

$$ (\partial_{tt}u_s, \partial_t v_f) \in L^2(0, T; V'_\phi). $$

These regularities, together with the density of $V_\phi$ in $H_\phi$, yield

$$ \frac{d^2}{dt^2} \int_{\Omega} \rho_s(1-\phi) u_s(t) \cdot w_s \, dx + \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \, dx = \left\langle (\rho_s(1-\phi) \partial_{tt} u_s(t), \rho_f \phi \partial_t v_f(t)), (w_s, w_f) \right\rangle_{V'_\phi, V_\phi} $$

in $D'(0, T)$, for all $(w_s, w_f) \in V_\phi$. Since functions in $V_\phi \otimes D(0, T)$ generate the space $L^2(0, T; V_\phi)$, we get

$$ \left\{ \begin{array}{l} \forall (w_s, w_f) \in L^2(0, T; V_\phi), \\ \displaystyle \int_0^T \left\langle (\rho_s(1-\phi) \partial_{tt} u_s, \rho_f \phi \partial_t v_f), (w_s, w_f) \right\rangle_{V'_\phi, V_\phi} \, dt + \int_0^T \int_\Omega \sigma_s(u_s) : \varepsilon(w_s) \, dx \, dt \\ \displaystyle \qquad + 2\eta \int_0^T \int_\Omega \varepsilon(v_s) : \varepsilon(w_s) \, dx \, dt + \int_0^T \int_\Omega \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx \, dt + \int_0^T \int_\Omega \phi^2 k_f^{-1}(v_f - v_s) \cdot (w_f - w_s) \, dx \, dt \\ \displaystyle \qquad - \int_0^T \int_\Omega \theta v_f \cdot w_f \, dx \, dt = \int_0^T \int_\Omega \rho_s(1-\phi) f \cdot w_s \, dx \, dt + \int_0^T \int_\Omega \rho_f \phi f \cdot w_f \, dx \, dt. \end{array} \right. $$

The energy estimate (17) then follows by choosing $(w_s, w_f) = (\partial_t u_s, v_f)$ above and guarantees uniqueness.

The equivalence between the variational and mild solutions can be proved in the same way as in the compressible case (see Proposition 2.8). We only need to notice that $V$ is dense in $H$.

To show the existence of pressure, we integrate (73b) in time over $(0, t)$ (see for instance [100] for a similar argument for the Stokes system). Let us define

$$ U_s(t) = \int_0^t u_s(s) \, ds, \quad V_s(t) = \int_0^t v_s(s) \, ds, \quad V_f(t) = \int_0^t v_f(s) \, ds \quad \text{and} \quad F(t) = \int_0^t f(s) \, ds, $$

it follows that

$$ \left\langle \rho_s(1-\phi)(v_s(t) - v_{s0}) - \operatorname{div}(\sigma_s^{\text{vis}}(V_s(t))) - \operatorname{div}(\sigma_s(U_s(t))) - \phi^2 k_f^{-1}(V_f(t) - V_s(t)) - \rho_s(1-\phi)F(t), w_s \right\rangle + \left\langle \rho_f\phi(v_f(t) - v_{f0}) - \operatorname{div}(\phi\sigma_f(V_f(t))) + \phi^2 k_f^{-1}(V_f(t) - V_s(t)) - \rho_f\phi F(t), w_f \right\rangle = 0 $$

for all $(w_s, w_f) \in V_\phi$. Combining Theorem 3.3 and Nečas Lemma provides the existence of $P \in C^0([0, T]; L_0^2(\Omega))$ such that

$$ \begin{aligned} & \rho_s(1-\phi)(v_s - v_{s0}) - \operatorname{div}(\sigma_s^{\text{vis}}(V_s)) - \operatorname{div}(\sigma_s(U_s)) - \phi^2 k_f^{-1}(V_f - V_s) - \rho_s(1-\phi)F = -(1-\phi)\nabla P, \\ & \rho_f\phi(v_f - v_{f0}) - \operatorname{div}(\phi\sigma_f(V_f)) + \phi^2 k_f^{-1}(V_f - V_s) - \rho_f\phi F = -\phi\nabla P. \end{aligned} $$

As a consequence, $p = \partial_t P$ satisfies (59) in the distribution sense. $\square$
**Remark 3.13.** We could also define the variational solution without assuming time continuity. Time continuity would then follows from the continuous injection of the space

$$W_{\phi}(0, T) = \{ (v_s, v_f) \in L^2(0, T; V_{\phi}) \text{ such that } (\partial_t v_s, \partial_t v_f) \in L^2(0, T; V'_{\phi}) \}$$

into $C^0([0, T]; H_{\phi})$.

## 3.3 The case $\eta = 0$

This case combines the two difficulties encountered earlier: the incompressibility constraint and the absence of solid viscosity. To handle this case, an option is to combine the functional framework adapted to the incompressibility constraint with the T-coercivity approach used in Section 2.3. This method provides the following result.

**Theorem 3.14.** If $\eta = 0$, then the conclusions of Theorem 3.8 and Corollary 3.10 remain true.

*Proof.* To prove that the operator $A_0^\infty + \lambda_0 I$ is surjective from $D(A_0^\infty)$ to $H$ for all $\lambda_0 > \omega$, we show that $a_0^\infty(\cdot, \cdot) + \lambda_0(\cdot, \cdot)_H$ is T-coercive for the mapping $T: (u_s, v_s, v_f) \mapsto (\frac{1}{2}u_s - \frac{1}{2\lambda_0}v_s, v_s, v_f)$ defined by (45). To do so, we reproduce exactly the same calculations as in the compressible case (see the proof of Theorem 2.11), but replacing $Y$ and $Z$ by $V$ and $H$ respectively. This mapping is bijective from $V$ into itself because it does not affect the velocity components and thus the mixture's divergence constraint. The rest of the proof is the same as in Theorem 3.8 and Corollary 3.10. □

Yet, we present here another approach to prove that $A_0^\infty + \lambda_0 I$ is surjective. This proof is based on a mixed formulation of the problem and is more suitable for numerical approximation. Indeed, the formulation (64) involves a constrained space $V_\phi$ that we would like to relax for numerical purpose, most numerical strategies relying on mixed formulations. Note moreover that the space $V_\phi$ depends on the porosity and thus on a specific data set. The mixed problem we would like to solve writes

$$\begin{cases} \text{Find } z \in Y_0 \text{ such that} \\ \forall y \in Y_0, \quad a_{\lambda_0}(z, y) = (g, y)_{Z_0}, \end{cases} \qquad (75)$$

with

$$\begin{aligned} a_{\lambda_0}(z, y) ={}& \lambda_0 \int_{\Omega} \sigma_s(u_s) : \varepsilon(d_s) \, dx - \int_{\Omega} \sigma_s(v_s) : \varepsilon(d_s) \, dx + \lambda_0 \int_{\Omega} \rho_s(1-\phi) v_s \cdot w_s \, dx + \int_{\Omega} \sigma_s(u_s) : \varepsilon(w_s) \, dx \\ &+ \int_{\Omega} \phi^2 k_f^{-1} (v_f - v_s) \cdot (w_f - w_s) \, dx + \lambda_0 \int_{\Omega} \rho_f \phi v_f \cdot w_f \, dx + \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx - \int_{\Omega} \theta v_f \cdot w_f \, dx \\ &- \int_{\Omega} p \operatorname{div} ((1-\phi) w_s + \phi w_f) \, dx + \int_{\Omega} \operatorname{div} ((1-\phi) v_s + \phi v_f) \, q \, dx \end{aligned}$$

for any $z = (u_s, v_s, v_f, p)$ and $y = (d_s, w_s, w_f, q)$ in $Y_0$, where $Z_0 = [H_0^1(\Omega)]^d \times [L^2(\Omega)]^d \times [L^2(\Omega)]^d \times L_0^2(\Omega)$ and $Y_0 = [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times L_0^2(\Omega)$. The spaces $Z_0$ and $Y_0$ are almost similar to $Z$ and $Y$ but include the additional condition $\int_{\Omega} p \, dx = 0$ that is required to ensure pressure uniqueness.

**Proposition 3.15.** Let $g \in Z_0$. If $\lambda_0 > (\rho_f \phi_{\min})^{-1} \|\theta\|_{L^\infty(\Omega)}$, then Problem (75) is well-posed in $Y_0$.

*Proof.* According to Proposition 2.10, it is sufficient to find $y^*$ depending continuously on $z$ such that the inequality $a_{\lambda_0}(z, y^*) \ge \alpha \|z\|_{Y_0}^2$ is satisfied for any $z \in Z_0$, with $\alpha > 0$.

From the properties of the divergence operator, we know, as already stated in (63), that for any $p \in L_0^2(\Omega)$, there exists $v_p \in [H_0^1(\Omega)]^d$ and $C_{\text{div}} > 0$ such that

$$\operatorname{div} v_p = p \quad \text{and} \quad \|\nabla v_p\|^2 \le C_{\text{div}} \|p\|^2.$$
For some constants $\alpha, \beta$ and $\gamma$ to be adjusted, we choose $d_s^* = \beta u_s + \gamma v_s$, $w_s^* = \alpha v_s - v_p$, $w_f^* = \alpha v_f - v_p$ and $q^* = \alpha p$. Thus

$$ 
\begin{aligned}
a_{\lambda_0}(z, y^*) ={}& \lambda_0 \int_{\Omega} \beta \sigma_s(u_s) : \varepsilon(u_s) \, dx + \lambda_0 \int_{\Omega} \gamma \sigma_s(u_s) : \varepsilon(v_s) \, dx - \int_{\Omega} \beta \sigma_s(v_s) : \varepsilon(u_s) \, dx - \int_{\Omega} \gamma \sigma_s(v_s) : \varepsilon(v_s) \, dx \\
& + \lambda_0 \int_{\Omega} \rho_s(1-\phi)(\alpha |v_s|^2 - v_s \cdot v_p) \, dx + \int_{\Omega} \alpha \sigma_s(u_s) : \varepsilon(v_s) \, dx - \int_{\Omega} \sigma_s(u_s) : \varepsilon(v_p) \, dx \\
& + \int_{\Omega} \alpha \phi^2 k_f^{-1}(v_f - v_s) \cdot (v_f - v_s) \, dx + \int_{\Omega} (\lambda_0 \rho_f \phi - \theta)(\alpha |v_f|^2 - v_f \cdot v_p) \, dx \\
& + \int_{\Omega} \phi (\alpha \sigma_f(v_f) : \varepsilon(v_f) - \sigma_f(v_f) : \varepsilon(v_p)) \, dx + \int_{\Omega} p \operatorname{div}((1-\phi)v_p + \phi v_p) \, dx.
\end{aligned}
$$

Note that the term $\int_{\Omega} p \operatorname{div}((1-\phi)v_p + \phi v_p) \, dx$ is equal to $\|p\|^2$ thanks to the choice of $v_p$. As in the proof of Theorem 2.11, we set $\beta = \frac{\alpha}{2}$ and $\gamma = -\frac{\alpha}{2\lambda_0}$ in order to remove the terms in the form $\int_{\Omega} \sigma_s(u_s) : \varepsilon(v_s) \, dx$. Consequently, we have

$$ 
\begin{aligned}
a_{\lambda_0}(z, y^*) \ge{}& \frac{\lambda_0 \alpha}{2} \int_{\Omega} \sigma_s(u_s) : \varepsilon(u_s) \, dx - \int_{\Omega} \sigma_s(u_s) : \varepsilon(v_p) \, dx + \frac{\alpha}{2\lambda_0} \int_{\Omega} \sigma_s(v_s) : \varepsilon(v_s) \, dx \\
& + \lambda_0 \rho_s(1-\phi_{\max}) \int_{\Omega} (\alpha |v_s|^2 - v_s \cdot v_p) \, dx + (\lambda_0 \rho_f \phi_{\min} - \|\theta\|_{L^{\infty}(\Omega)}) \int_{\Omega} (\alpha |v_f|^2 - v_f \cdot v_p) \, dx \\
& + \phi_{\min} \int_{\Omega} (\alpha \sigma_f(v_f) : \varepsilon(v_f) - \sigma_f(v_f) : \varepsilon(v_p)) \, dx + \int_{\Omega} p^2 \, dx.
\end{aligned}
\quad (76) $$

We choose $\lambda_0$ such that $\lambda_0 \rho_f \phi_{\min} - \|\theta\|_{L^{\infty}(\Omega)} > 0$. Next, for all $\delta > 0$, Young inequality yields

$$ 
\begin{aligned}
& - \int_{\Omega} \sigma_s(u_s) : \varepsilon(v_p) \, dx \geq -\frac{\delta}{2} \int_{\Omega} \sigma_s(u_s) : \varepsilon(u_s) \, dx - \frac{1}{2\delta} \int_{\Omega} \sigma_s(v_p) : \varepsilon(v_p) \, dx, \\
& - \int_{\Omega} \sigma_f(v_f) : \varepsilon(v_p) \, dx \geq -\frac{\delta}{2} \int_{\Omega} \sigma_f(v_f) : \varepsilon(v_f) \, dx - \frac{1}{2\delta} \int_{\Omega} \sigma_f(v_p) : \varepsilon(v_p) \, dx, \\
& \qquad - \int_{\Omega} v_s \cdot v_p \, dx \geq -\frac{\delta}{2} \int_{\Omega} |v_s|^2 \, dx - \frac{1}{2\delta} \int_{\Omega} |v_p|^2 \, dx, \\
& \qquad - \int_{\Omega} v_f \cdot v_p \, dx \geq -\frac{\delta}{2} \int_{\Omega} |v_f|^2 \, dx - \frac{1}{2\delta} \int_{\Omega} |v_p|^2 \, dx.
\end{aligned}
\quad (77) $$

Furthermore, it holds

$$ 
\begin{aligned}
\|v_p\|^2 &\le C_p \|\nabla v_p\|^2 \le C_p C_{\operatorname{div}} \|p\|^2, \\
\int_{\Omega} \sigma_f(v_p) : \varepsilon(v_p) \, dx &= \lambda_f \|\operatorname{div} v_p\|^2 + 2\mu_f \|\varepsilon(v_p)\|^2 \le (\lambda_f + 2\mu_f C_{\operatorname{div}}) \|p\|^2, \\
\int_{\Omega} \sigma_s(v_p) : \varepsilon(v_p) \, dx &= \lambda \|\operatorname{div} v_p\|^2 + 2\mu \|\varepsilon(v_p)\|^2 \le (\lambda + 2\mu C_{\operatorname{div}}) \|p\|^2,
\end{aligned}
\quad (78) $$

where $C_p$ denotes the constant of Poincaré inequality.

Using (77) and (78) to bound from below the right-hand side of (76) and rearranging terms, we obtain

$$ 
\begin{aligned}
a_{\lambda_0}(z, y^*) \ge{}& \left(\frac{\lambda_0 \alpha}{2} - \frac{\delta}{2}\right) \|u_s\|_s^2 + \frac{\alpha}{2\lambda_0} \|v_s\|_s^2 + \lambda_0 \rho_s(1-\phi_{\max}) \left(\alpha - \frac{\delta}{2}\right) \|v_s\|^2 \\
& + (\lambda_0 \rho_f \phi_{\min} - \|\theta\|_{L^{\infty}(\Omega)}) \left(\alpha - \frac{\delta}{2}\right) \|v_f\|^2 + 2\mu_f \phi_{\min} \left(\alpha - \frac{\delta}{2}\right) \|\varepsilon(v_f)\|^2 + \left(1 - \frac{\delta^*}{2\delta}\right) \|p\|^2,
\end{aligned}
$$

where $\delta^* = \lambda + 2\mu C_{\operatorname{div}} + \lambda_0 \rho_s(1-\phi_{\max})C_p C_{\operatorname{div}} + (\lambda_0 \rho_f \phi_{\min} - \|\theta\|_{L^{\infty}(\Omega)})C_p C_{\operatorname{div}} + \phi_{\min}(\lambda_f + 2\mu_f C_{\operatorname{div}})$, $\delta^* > 0$.
Hence, setting $\delta = \delta^*$ and $\alpha = \alpha^* = \max(\delta^*, \frac{2\delta^*}{\lambda_0})$, we get

$$a_{\lambda_0}(z, y^*) \geq \frac{\delta^*}{2} \|u_s\|_s^2 + \frac{\alpha^*}{2\lambda_0} \|v_s\|_s^2 + \mu_f \phi_{\min} \delta^* \|\varepsilon(v_f)\|^2 + \frac{1}{2} \|p\|^2.$$

Finally, we infer that $a_{\lambda_0}$ is T-coercive for the mapping

$$T : (u_s, v_s, v_f, p) \mapsto \left( \frac{\alpha^*}{2} u_s - \frac{\alpha^*}{2\lambda_0} v_s, \alpha^* v_s - v_p, \alpha^* v_f - v_p, \alpha^* p \right), \quad (79)$$

which is bijective since $p \mapsto v_p$ is a bijection. $\square$

**Remark 3.16.** This mixed formulation is also applicable to the case $\kappa = +\infty$ and $\eta > 0$. In that case, the proof can be simplified by considering the mapping $T : (u_s, v_s, v_f, p) \mapsto (u_s, \alpha^* v_s - v_p, \alpha^* v_f - v_p, \alpha^* p)$.

**Remark 3.17.** The mixed formulation is equivalent to the constrained formulation thanks to the inf-sup property proved in Proposition 3.1. Note moreover that, as for the proof of Proposition 3.1, the T-coercivity only relies on the standard inf-sup condition for the divergence operator and therefore is independent of the porosity $\phi$.

Finally, as for the compressible inviscid case, we can prove the existence of a variational solution.

**Theorem 3.18.** Assume that (h1) – (h4) are satisfied, $\eta = 0$, $\int_{\Omega} \theta \, dx = 0$ and $z_0 = (u_{s0}, v_{s0}, v_{f0}) \in H$. Then there exists a unique variational solution $u_s \in C^0([0, T]; [H_0^1(\Omega)]^d)$ and $(\partial_t u_s, v_f) \in C^0([0, T]; H_\phi)$ with $v_f \in L^2(0, T; [H_0^1(\Omega)]^d)$ such that

$$(u_s(0), \partial_t u_s(0), v_f(0)) = (u_{s0}, v_{s0}, v_{f0})$$

and the following equations holds, in $D'(0, T)$,

$$\left\{
\begin{array}{l}
\displaystyle \frac{d^2}{dt^2} \int_{\Omega} \rho_s (1-\phi) u_s(t) \cdot w_s \, dx + \int_{\Omega} \sigma_s(u_s(t)) : \varepsilon(w_s) \, dx \\
\displaystyle \quad + \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f(t) \cdot w_f \, dx + \int_{\Omega} \phi \sigma_f(v_f(t)) : \varepsilon(w_f) \, dx + \int_{\Omega} \phi^2 k_f^{-1} (v_f(t) - \partial_t u_s(t)) \cdot (w_f - w_s) \, dx \\
\displaystyle \quad - \int_{\Omega} \theta v_f(t) \cdot w_f \, dx = \int_{\Omega} \rho_s (1-\phi) f(t) \cdot w_s \, dx + \int_{\Omega} \rho_f \phi f(t) \cdot w_f \, dx, \quad \forall (w_s, w_f) \in V_{\phi}.
\end{array}
\right.
\quad (80)$$

The energy estimate (17) holds true (with $\eta = 0$) and the variational solution coincides with the mild solution. Furthermore, there exists a unique pressure $p$ such that $(u_s, v_s, v_f, p)$ satisfies (59) in the distribution sense, with $v_s = \partial_t u_s$.

*Proof.* We follow the same steps as in Theorem 2.12, but within the functional framework adapted to the incompressibility constraint. Existence of solutions is obtained by an approximated sequence of strong solutions $z^n = (u_s^n, v_s^n, v_f^n) \in C^1([0, T]; H) \cap C^0([0, T]; D(A_0^\infty))$ verifying

$$\left\{
\begin{array}{l}
\dot{z}^n(t) + A_0^\infty z^n(t) = \Pi g^n(t), \quad t \in [0, T], \\
z^n(0) = z_0^n,
\end{array}
\right.
$$

where $\Pi g^n = \Pi(0, f^n, f^n) \in H^1(0, T; H)$ and $z_0^n \in D(A_0^\infty)$ denote respectively an approximation of source
terms and initial conditions. This sequence of solution satisfies the variational formulation

$$ (VF)^n_\infty \left\{ \begin{array}{l} \forall t \in [0, T], \forall (d_s, w_s, w_f) \in V = [H_0^1(\Omega)]^d \times V_\phi, \\ \int_\Omega \sigma_s(\partial_t u_s^n(t)) : \varepsilon(d_s) \, dx = \int_\Omega \sigma_s(v_s^n(t)) : \varepsilon(d_s) \, dx, \qquad (81a) \\ \int_\Omega \rho_s(1-\phi) \partial_t v_s^n(t) \cdot w_s \, dx + \int_\Omega \sigma_s(u_s^n(t)) : \varepsilon(w_s) \, dx \\ \qquad + \int_\Omega \rho_f \phi \partial_t v_f^n(t) \cdot w_f \, dx + \int_\Omega \phi \sigma_f(v_f^n(t)) : \varepsilon(w_f) \, dx + \int_\Omega \phi^2 k_f^{-1}(v_f^n(t) - v_s^n(t)) \cdot (w_f - w_s) \, dx \\ \qquad - \int_\Omega \theta^n v_f^n(t) \cdot w_f \, dx = \int_\Omega \rho_s(1-\phi) f^n(t) \cdot w_s \, dx + \int_\Omega \rho_f \phi f^n(t) \cdot w_f \, dx. \qquad (81b) \end{array} \right. $$

Moreover, by taking $(u_s^n, v_s^n, v_f^n) \in C^0([0, T]; V)$ as test functions in (81a), we get that it satisfies the energy estimate (17) with $\eta = 0$. Hence $z^n$ is a Cauchy sequence in $C^0([0, T]; H)$ and $v_f^n$ is a Cauchy sequence in $L^2(0, T; [H_0^1(\Omega)]^d)$, which allows us to pass to the limit in (81b). For a given $d_s \in [L^2(\Omega)]^d$, we choose the unique solution $\eta_s \in [H_0^1(\Omega)]^d$ of $-\mathrm{div}(\sigma_s(\eta_s)) = d_s$ as a test function in (81a), which yields

$$ \forall t \in [0, T], \forall d_s \in [L^2(\Omega)]^d, \quad \frac{d}{dt} \int_{\Omega} u_s(t) \cdot d_s \, dx - \int_{\Omega} v_s(t) \cdot d_s \, dx = 0 $$

after passing to the limit. Putting these two limit formulations together gives (80). As for the viscous case the equation (80) implies that

$$ (\partial_{tt} u_s, \partial_t v_f) \in L^2(0, T; V'_\phi), $$

and for any test functions $(w_s, w_f)$ in $V_\phi$

$$ \frac{d^2}{dt^2} \int_{\Omega} \rho_s(1-\phi) u_s \cdot w_s \, dx + \frac{d}{dt} \int_{\Omega} \rho_f \phi v_f \cdot w_f \, dx = \langle (\rho_s(1-\phi) \partial_{tt} u_s, \rho_f \phi \partial_t v_f), (w_s, w_f) \rangle_{V'_\phi, V_\phi}. $$

To show uniqueness, we are going to use the same Ladyzhenskaya test functions as in the compressible case. The difficulty then lies in justifying that the calculations done in the compressible case remain valid in the constrained functional setting. Let $(u_s, v_f)$ be a solution to (80) with zero initial conditions and source terms, and let $\tau$ be given in $(0, T)$. We first write the weak space-time variational formulation satisfied by $(u_s, v_f)$. In a standard way, by multiplying the weak formulation (80) by a $\psi \in H^1(0, T)$ such that $\phi(T) = 0$ and integrating over $(0, T)$ and by parts in time we obtain

$$ \begin{aligned} & - \int_0^T \int_\Omega \rho_s(1-\phi) \partial_t u_s \cdot \partial_t \psi(t) w_s \, dx \, dt + \int_0^T \int_\Omega \sigma_s(u_s) : \varepsilon(w_s) \psi(t) \, dx \, dt \\ & \qquad - \int_0^T \int_\Omega \rho_f \phi v_f \cdot \partial_t \psi(t) w_f \, dx \, dt \\ & + \int_0^T \int_\Omega \phi \sigma_f(v_f) : \varepsilon(w_f) \psi(t) \, dx \, dt \\ & \qquad + \int_0^T \int_\Omega \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (w_f - w_s) \psi(t) \, dx \, dt \\ & \qquad - \int_0^T \int_\Omega \theta v_f \cdot w_f \psi(t) \, dx \, dt = 0. \end{aligned} $$

Since linear combinations of functions of the type $(\psi(t)w_s, \psi(t)w_f)$ with $\psi \in H^1(0, T)$ such that $\phi(T) = 0$ and $(w_s, w_f) \in V_\phi$ are dense in the space of functions $w$ of $H^1(0, T; V_\phi)$ such that $w(T) = 0$, we obtain

$$ \left\{ \begin{array}{l} - \int_0^T \int_\Omega \rho_s(1-\phi) \partial_t u_s \cdot \partial_t w_s \, dx \, dt + \int_0^T \int_\Omega \sigma_s(u_s) : \varepsilon(w_s) \, dx \, dt \\ \\ - \int_0^T \int_\Omega \rho_f \phi v_f \cdot \partial_t \psi(t) w_f \, dx \, dt + \int_0^T \int_\Omega \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx \, dt \\ \\ + \int_0^T \int_\Omega \phi^2 k_f^{-1} (v_f - \partial_t u_s) \cdot (w_f - w_s) \, dx \, dt - \int_0^T \int_\Omega \theta v_f \cdot w_f \, dx \, dt = 0, \\ \\ \forall (w_s, w_f) \in H^1(0, T; V_\phi) \text{ such that } w_s(T) = w_f(T) = 0. \end{array} \right. \qquad (82) $$
Then, we consider the same tests functions as in the compressible case, namely

$$ \psi_s(t) = \begin{cases} -\int_t^\tau u_s(\sigma) \, d\sigma & \text{if } \tau \ge t \\ 0 & \text{if } \tau \le t \end{cases} \quad \text{and} \quad \psi_f(t) = \begin{cases} -\int_t^\tau \int_0^v v_f(\sigma) \, d\sigma \, dv & \text{if } \tau \ge t \\ 0 & \text{if } \tau \le t. \end{cases} $$

These test functions are still admissible here. Indeed, we know that $(\partial_t u_s, v_f) \in C^0([0, T]; H_\phi)$. Recalling the characterization of the space $H_\phi$ established in Proposition 3.4, it follows that

$$ \mathrm{div}((1-\phi)\partial_t u_s + \phi v_f) = 0, \quad \text{in } C^0([0, T]; \mathcal{D}'(\Omega)). $$

Note that, as in the inviscid compressible case, $\mathrm{div}((1-\phi)\partial_t u_s + \phi v_f)$ belongs to $L^2((0, T); H^{-1}(\Omega))$. Next, by integrating two times in time, we obtain

$$ \mathrm{div} \left( (1-\phi) \left( -\int_t^\tau \int_0^v \partial_t u_s(\sigma) \, d\sigma \, dv \right) + \phi \left( -\int_t^\tau \int_0^v v_f(\sigma) \, d\sigma \, dv \right) \right) = 0. $$

Since $u_s(0) = 0$, we conclude that $(\psi_s, \psi_f) \in C^1([0, T]; V_\phi)$.

Choosing $(\psi_s, \psi_f)$ as test functions, the calculations are exactly the same as in the compressible case – see (58) – but without the pressure terms, and with $\theta$ independent of time. We get

$$ \frac{1}{2} \int_{\Omega} \rho_s (1-\phi) |u_s(\tau)|^2 \, dx + \frac{1}{2} \int_{\Omega} \rho_f \phi |\partial_t \psi_f(\tau)|^2 \, dx + \frac{1}{2} \int_{\Omega} \sigma_s(\psi_s(0)) : \varepsilon(\psi_s(0)) \, dx \\ + \int_{0}^{\tau} \int_{\Omega} \phi \sigma_f(\partial_t \psi_f) : \varepsilon(\partial_t \psi_f) \, dx \, dt + \int_{0}^{\tau} \int_{\Omega} \phi^2 k_f^{-1} (\partial_t \psi_f - \partial_t \psi_s)^2 \, dx \, dt = \int_{0}^{\tau} \int_{\Omega} \theta |\partial_t \psi_f|^2 \, dx \, dt. $$

Estimating the right-hand side as in the viscous case shows that $u_s = v_f = 0$ by an application of Grönwall Lemma.

Using Duhamel formula (70), one shows exactly as in the compressible case that the sequence $z^n$ also converges towards the mild solution in $C^0([0, T]; H)$, so that the mild solution coincides with the variational solution built from the same approximation process. Finally, the existence of a pressure $p$ such that $(u_s, v_s, v_f, p)$ satisfies (59) in the distribution sense is obtained by combining Theorem 3.3 and Nečas Lemma, like in the proof of Theorem 3.11. $\square$

# 4 Incompressible limit

In this section, we show how to pass to the limit in the weak formulation for $\kappa < +\infty$ as $\kappa$ goes to infinity and obtain the incompressible system. Similar incompressible limits were considered for Biot's consolidation model, both in linear [94] or non-linear [24] regimes, and the influence of compressibility was analyzed in the 1D linear case [21].

For this purpose, we need to get an energy estimate independent of $\kappa$ in the compressible case. This can be achieved by lifting the right-hand side of the pressure equation. We consider $v_{\theta,\alpha} = \frac{1}{\alpha}v_\theta$ where $v_\theta$ is defined by (61), so that

$$ \mathrm{div} (\alpha v_{\theta,\alpha}) = \rho_f^{-1} \theta. \qquad (83) $$

Note that to ensure that $v_{\theta,\alpha}$ is in $[H^2(\Omega)]^d$, we need more assumptions and more regularity on the Biot-Willis coefficient $\alpha$. Therefore, (h6) becomes

$$ (h6)_{\text{bis}} \begin{cases} \alpha \in H^{d/2+r}(\Omega) & \text{with } d/2+r \ge 2, \\ \forall x \in \Omega, \quad 0 < (\alpha - \phi)_{\min} \le \alpha(x) - \phi(x) \le (\alpha - \phi)_{\max} < 1. \end{cases} $$
As already noticed in Remark 3.9, in order to have the adequate regularity for the right-hand side of the equation verified by the new unknowns, such a lifting requires additional assumptions on the fluid mass input $\theta$, namely

$$ (h5)_{\text{bis}} \begin{cases} \theta \in C^0([0, T] \times \Omega) \cap H^1(0, T; L^2(\Omega)) \cap L^2(0, T; H^1(\Omega)), \\ \forall t \in [0, T], \quad \int_{\Omega} \theta(t) dx = 0. \end{cases} $$

Consequently, under $(h5)_{\text{bis}}$ and $(h6)_{\text{bis}}$, the lifting $v_{\theta,\alpha}$ satisfying (83) belongs to $H^1(0, T; [H_0^1(\Omega)]^d) \cap L^2(0, T; [H^2(\Omega)]^d)$, and we make the change of variables $(u_s, v_f, p) \mapsto (u_s - \int_0^t v_{\theta,\alpha}(s) ds, v_f - v_{\theta,\alpha}, p)$ so that the right-hand side of the pressure equation reduces to zero.

In order to recover the initial conditions in the incompressible limit, we are not going to pass to the limit in (37) and (47) which are written in $D'(0, T)$, but rather in the following weak formulation: for any $(u_{s0}, v_{s0}, v_{f0}, p_0) \in Z$ and $f \in L^2(0, T; [L^2(\Omega)]^d)$, find $u_s^\kappa \in C^0([0, T]; [H_0^1(\Omega)]^d)$, $\partial_t u_s^\kappa \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ if $\eta > 0$ or $\partial_t u_s^\kappa \in C^0([0, T]; [L^2(\Omega)]^d)$ if $\eta = 0$, $v_f^\kappa \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ and $p^\kappa \in C^0([0, T]; L^2(\Omega))$ such that

$$ \left\{ \begin{aligned} & \int_0^T \int_{\Omega} \rho_s(1-\phi) u_s^\kappa \cdot \partial_{tt} w_s \,dx \,dt + \int_0^T \int_{\Omega} \sigma_s(u_s^\kappa) : \varepsilon(w_s) \,dx \,dt + 2\eta \int_0^T \int_{\Omega} \varepsilon(\partial_t u_s^\kappa) : \varepsilon(w_s) \,dx \,dt \\ & \qquad - \int_0^T \int_{\Omega} \phi^2 k_f^{-1}(v_f^\kappa - \partial_t u_s^\kappa) \cdot w_s \,dx \,dt - \int_0^T \int_{\Omega} p^\kappa \operatorname{div}((\alpha-\phi)w_s) \,dx \,dt \\ & \qquad = \int_0^T \int_{\Omega} \rho_s(1-\phi)f \cdot w_s \,dx \,dt + \int_{\Omega} \rho_s(1-\phi) v_{s0} \cdot w_s(0) \,dx - \int_{\Omega} \rho_s(1-\phi) u_{s0} \cdot \partial_t w_s(0) \,dx, \quad (84a) \\ & \int_0^T \int_{\Omega} \rho_f \phi v_f^\kappa \cdot \partial_t w_f \,dx \,dt + \int_0^T \int_{\Omega} \phi \sigma_f(v_f^\kappa) : \varepsilon(w_f) \,dx \,dt \\ & \qquad + \int_0^T \int_{\Omega} \phi^2 k_f^{-1}(v_f^\kappa - \partial_t u_s^\kappa) \cdot w_f \,dx \,dt - \int_0^T \int_{\Omega} \theta v_f^\kappa \cdot w_f \,dx \,dt - \int_0^T \int_{\Omega} p^\kappa \operatorname{div}(\phi w_f) \,dx \,dt \\ & \qquad = \int_0^T \int_{\Omega} \rho_f \phi f \cdot w_f \,dx \,dt + \int_{\Omega} \rho_f \phi v_{f0} \cdot w_f(0) \,dx, \quad (84b) \end{aligned} \right. $$

and

$$ \left\{ \begin{aligned} & - \int_0^T \int_{\Omega} \frac{\alpha - \phi}{\kappa} p^\kappa \partial_t q \,dx \,dt + \int_0^T \int_{\Omega} \operatorname{div}((\alpha - \phi) \partial_t u_s^\kappa) q \,dx \,dt \\ & \quad + \int_0^T \int_{\Omega} \operatorname{div}(\phi v_f^\kappa) q \,dx \,dt = \int_{\Omega} \frac{\alpha - \phi}{\kappa} p_0 q(0) \,dx \qquad \text{if } \eta > 0, \quad (85a) \\ & - \int_0^T \int_{\Omega} \frac{\alpha - \phi}{\kappa} p^\kappa \partial_t q \,dx \,dt - \int_0^T \int_{\Omega} \operatorname{div}((\alpha - \phi) u_s^\kappa) \partial_t q \,dx \,dt \\ & \quad + \int_0^T \int_{\Omega} \operatorname{div}(\phi v_f^\kappa) q \,dx \,dt = \int_{\Omega} \frac{\alpha - \phi}{\kappa} p_0 q(0) \,dx + \int_{\Omega} \operatorname{div}((\alpha - \phi) u_{s0}) q(0) \,dx \qquad \text{if } \eta = 0, \quad (85b) \end{aligned} \right. $$

for all admissible test functions

$$ \left\{ \begin{aligned} & w_s \in H^2(0, T; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d), \\ & w_f \in H^1(0, T; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d), \\ & q \in H^1(0, T; L^2(\Omega)), \\ & w_s(T) = \partial_t w_s(T) = w_f(T) = q(T) = 0. \end{aligned} \right. \quad (86) $$

The main difference between the weak formulation (84) – (85) and (37) or (47) is that the test functions depend on space but also on time. Besides, the initial conditions are weakly imposed in (84) – (85), while they are strongly imposed in (36) or (46). This space-time weak formulation can be obtained from (37) or (47) with the same arguments used to derive (82).
**Remark 4.1.** By choosing $(w_s, w_f, q) = (\hat{w}_s(x), \hat{w}_f(x), \hat{q}(x)) \psi(t)$ with $(\hat{w}_s, \hat{w}_f, \hat{q}) \in [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times L^2(\Omega)$ and $\psi \in D(0, T)$, we see that the weak formulation (84) – (85) implies the variational formulation (37) or (47). Hence, from the uniqueness of the variational solution, the solutions of these two formulations coincide.

We are now ready to establish how the solution in the compressible case converges towards the solution in the incompressible regime as $\kappa$ goes to infinity.

**Theorem 4.2.** Assume that (h1) – (h4), (h5)$_{bis}$ and (h6)$_{bis}$ are satisfied. For $z_0 = (u_{s0}, v_{s0}, v_{f0}) \in Z$, let $(u_s^\kappa, v_f^\kappa, p^\kappa)$ be the solution of (84) – (85). As $\kappa$ goes to infinity, $(u_s^\kappa, \partial_t u_s^\kappa, v_f^\kappa)$ converge weakly towards the solution of the following formulation: find $u_s \in C^0([0, T]; [H_0^1(\Omega)]^d)$, $\partial_t u_s \in C^0([0, T]; [L^2(\Omega)]^d)$ and $v_f \in C^0([0, T]; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$ or $(\partial_t u_s, v_f) \in L^2(0, T; V_\phi)$ if $\eta > 0$, such that

$$ egin{cases} 
\displaystyle \int_0^T \int_{\Omega} \rho_s(1-\phi) u_s \cdot \partial_{tt} w_s \, dx \, dt + \int_0^T \int_{\Omega} \sigma_s(u_s) : \varepsilon(w_s) \, dx \, dt + 2\eta \int_0^T \int_{\Omega} \varepsilon(\partial_t u_s) : \varepsilon(w_s) \, dx \, dt \\ 
- \int_0^T \int_{\Omega} \rho_f \phi v_f \cdot \partial_t w_f \, dx \, dt + \int_0^T \int_{\Omega} \phi \sigma_f(v_f) : \varepsilon(w_f) \, dx \, dt + \int_0^T \int_{\Omega} \phi^2 k_f^{-1}(v_f - \partial_t u_s) \cdot (w_f - w_s) \, dx \, dt \\ 
- \int_0^T \int_{\Omega} \theta v_f \cdot w_f \, dx \, dt = \int_0^T \int_{\Omega} \rho_s(1-\phi)f \cdot w_s \, dx \, dt + \int_0^T \int_{\Omega} \rho_f \phi f \cdot w_f \, dx \, dt \\ 
+ \int_{\Omega} \rho_s(1-\phi) v_{s0} \cdot w_s(0) \, dx - \int_{\Omega} \rho_s(1-\phi) u_{s0} \cdot \partial_t w_s(0) \, dx + \int_{\Omega} \rho_f \phi v_{f0} \cdot w_f(0) \, dx, 
\end{cases} 
\quad (87) $$

and

$$ - \int_0^T \int_{\Omega} \mathrm{div} \left( (1-\phi) u_s \right) \partial_t q \, dx \, dt + \int_0^T \int_{\Omega} \mathrm{div} \left( \phi v_f \right) q \, dx \, dt = \int_{\Omega} \mathrm{div} \left( (1-\phi) u_{s0} \right) q(0) \, dx \quad \text{if } \eta = 0, \quad (88) $$

for all admissible test functions

$$ egin{cases} 
w_s \in H^2(0, T; [L^2(\Omega)]^d), \\ w_f \in H^1(0, T; [L^2(\Omega)]^d), \\ (w_s, w_f) \in L^2(0, T; V_\phi), \\ q \in H^1(0, T; L^2(\Omega)), \\ w_s(T) = \partial_t w_s(T) = w_f(T) = q(T) = 0. 
\end{cases} 
\quad (89) $$

**Proof.** Let us prove this result in the inviscid case $\eta = 0$, the viscous case being similar. Since we can lift the mixture’s divergence constraint as in (83), let us consider the case where the right-hand side of the pressure equation is equal to zero. The resulting energy estimate reads

$$ 
\begin{aligned} 
& \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |\partial_t u_s^\kappa(t)|^2 dx + \frac{1}{2} \int_{\Omega} \sigma_s(u_s^\kappa(t)) : \varepsilon(u_s^\kappa(t)) dx + 2\eta \int_0^t \int_{\Omega} \varepsilon(\partial_t u_s^\kappa) : \varepsilon(\partial_t u_s^\kappa) dx ds + \frac{\rho_f}{2} \int_{\Omega} \phi |v_f^\kappa(t)|^2 dx \\ 
& + \int_0^t \int_{\Omega} \phi \sigma_f(v_f^\kappa) : \varepsilon(v_f^\kappa) dx ds + \int_0^t \int_{\Omega} \phi^2 k_f^{-1} (v_f^\kappa - \partial_t u_s^\kappa) \cdot (v_f^\kappa - \partial_t u_s^\kappa) dx ds + \frac{1}{2} \int_{\Omega} \frac{\alpha - \phi}{\kappa} |p^\kappa(t)|^2 dx \\ 
\le & \exp \left( \max \left( 1, \frac{2\|\theta\|_{C^0([0,T]\times\Omega)}}{\rho_f \phi_{\min}} \right) t \right) \left( \left( \frac{\rho_s}{2}(1-\phi_{\min}) + \frac{\rho_f}{2}\phi_{\max} \right) \int_0^t \int_{\Omega} |f|^2 dx ds + \frac{\rho_s}{2} \int_{\Omega} (1-\phi) |v_{s0}|^2 dx \right. \\ 
& \left. + \frac{1}{2} \int_{\Omega} \sigma_s(u_{s0}) : \varepsilon(u_{s0}) dx + \frac{\rho_f}{2} \int_{\Omega} \phi |v_{f0}|^2 dx + \frac{1}{2} \int_{\Omega} \frac{\alpha - \phi}{\kappa} |p_0|^2 dx \right). 
\end{aligned} 
\quad (90) $$
We deduce that, up to subsequences, the following weak convergences hold true as $\kappa$ goes to infinity:

$$
\begin{align*}
u_s^\kappa &\to u_s^\infty \quad \text{weakly star in } L^\infty(0, T; [H_0^1(\Omega)]^d), \\
\partial_t u_s^\kappa &\to \partial_t u_s^\infty \quad \text{weakly star in } L^\infty(0, T; [L^2(\Omega)]^d), \\
v_f^\kappa &\to v_f^\infty \quad \text{weakly star in } L^\infty(0, T; [L^2(\Omega)]^d), \\
v_f^\kappa &\to v_f^\infty \quad \text{weakly in } L^2(0, T; [H_0^1(\Omega)]^d),
\end{align*}
$$

for some elements $u_s^\infty \in L^\infty(0, T; [H_0^1(\Omega)]^d)$ with $\partial_t u_s^\kappa \in L^\infty(0, T; [L^2(\Omega)]^d)$ and $v_f^\infty \in L^\infty(0, T; [L^2(\Omega)]^d) \cap L^2(0, T; [H_0^1(\Omega)]^d)$.

We have no bound on the pressure $p^\kappa$ but (90) implies that $\frac{p^\kappa}{\sqrt{\kappa}}$ is bounded in $L^\infty(0, T; L^2(\Omega))$. Hence, we can select a subsequence (still denoted by $p^\kappa$) such that $\frac{p^\kappa}{\sqrt{\kappa}}$ converges in the weak-$^*$ topology of $L^\infty(0, T; L^2(\Omega))$.

By adding (84a) to (84b) and by restricting the velocities test functions ($w_s$, $w_f$) to functions in $L^2(0, T; V_\phi)$, it follows that $(u_s^\kappa, v_f^\kappa, p^\kappa)$ satisfies

$$
\left\{
\begin{aligned}
& \int_0^T \int_\Omega \rho_s (1-\phi) u_s^\kappa \cdot \partial_{tt} w_s \,dx\,dt + \int_0^T \int_\Omega \sigma_s (u_s^\kappa) : \varepsilon(w_s) \,dx\,dt + 2\eta \int_0^T \int_\Omega \varepsilon(\partial_t u_s^\kappa) : \varepsilon(w_s) \,dx\,dt \\
& - \int_0^T \int_\Omega \rho_f \phi v_f^\kappa \cdot \partial_t w_f \,dx\,dt + \int_0^T \int_\Omega \phi \sigma_f (v_f^\kappa) : \varepsilon(w_f) \,dx\,dt + \int_0^T \int_\Omega \phi^2 k_f^{-1} (v_f^\kappa - \partial_t u_s^\kappa) \cdot (w_f - w_s) \,dx\,dt \\
& - \int_0^T \int_\Omega \theta v_f^\kappa \cdot w_f \,dx\,dt = \int_0^T \int_\Omega \rho_s (1-\phi) f \cdot w_s \,dx\,dt + \int_0^T \int_\Omega \rho_f \phi f \cdot w_f \,dx\,dt \\
& + \int_\Omega \rho_s (1-\phi) v_{s0} \cdot w_s (0) \,dx - \int_\Omega \rho_s (1-\phi) u_{s0} \cdot \partial_t w_s (0) \,dx + \int_\Omega \rho_f \phi v_{f0} \cdot w_f (0) \,dx, \\
& - \int_0^T \int_\Omega \frac{\alpha-\phi}{\kappa} p^\kappa \partial_t q \,dx\,dt - \int_0^T \int_\Omega \operatorname{div} ((\alpha-\phi) u_s^\kappa) \partial_t q \,dx\,dt \\
& + \int_0^T \int_\Omega \operatorname{div} (\phi v_f^\kappa) q \,dx\,dt = \int_\Omega \frac{\alpha-\phi}{\kappa} p_0 q (0) \,dx + \int_\Omega \operatorname{div} ((\alpha-\phi) u_{s0}) q (0) \,dx,
\end{aligned}
\right.
$$

for all admissible test functions verifying (89).

Thus we can pass to the weak limit in this formulation by noting that $\alpha - \phi \to 1 - \phi$,

$$
\int_0^T \int_\Omega \frac{\alpha - \phi}{\kappa} p^\kappa \partial_t q \, dx \, dt = \frac{1}{\sqrt{\kappa}} \int_0^T \int_\Omega (\alpha - \phi) \frac{p^\kappa}{\sqrt{\kappa}} \partial_t q \, dx \, dt \to 0,
$$

and that

$$
\int_{\Omega} \frac{\alpha - \phi}{\kappa} p_0 q(0) \, dx \to 0
$$

as $\kappa$ goes to infinity.

In conclusion, $(u_s^\infty, v_f^\infty)$ satisfies exactly (87) and (88) in the incompressible limit. Moreover $u_s^\infty$, $\partial_t u_s^\infty$ and $v_f^\infty$ are continuous functions in time because they also satisfy (80) and hence coincide with the mild solution. Indeed, (80) can be recovered from (87) – (88) by taking admissible test functions of the form $(w_s, w_f, q) = (\hat{w}_s(x), \hat{w}_f(x), \hat{q}(x)) \psi(t)$ with $(\hat{w}_s, \hat{w}_f, \hat{q}) \in V_\phi \times L^2(\Omega)$ and $\psi \in D(0, T)$. $\square$

**Remark 4.3.** In the case where the right-hand side of the pressure equation is not equal to zero, we need to perform a lifting. Note that without this lifting step the energy estimate does not provide a uniform bound in $\kappa$ as $\kappa$ goes to infinity because of the coefficient $\frac{p^\kappa}{2\rho_f^2(\alpha-\phi)_{\min}}$ appearing in the right-hand side of (11). Moreover once the lifting is performed under assumptions (h5)$_{bis}$ and (h6)$_{bis}$, the new right-hand sides of the structure and fluid equations depend on $\alpha$. Yet, it is easy to verify that they converge strongly in the proper spaces ensuring the convergence of the right-hand sides as $\alpha$ goes to one in $H^{d/2+r}(\Omega)$ with $d/2 + r \ge 2$.
**Remark 4.4.** Theorem 4.2 provides the weak convergence of the displacement and velocities in the incompressible limit. If the incompressible regime solution is more regular, we can also obtain the pressure convergence and recover strong convergence for the displacement and velocities. More precisely, following the same guidelines as in [58, Lemma 75.1], we can show that

$$ \|u_s^\kappa - u_s^\infty\|_{L^\infty(0,T;[H_0^1(\Omega)]^d)}^2 + \|\partial_t u_s^\kappa - \partial_t u_s^\infty\|_{L^\infty(0,T;[L^2(\Omega)]^d)}^2 + \|v_f^\kappa - v_f^\infty\|_{L^\infty(0,T;[L^2(\Omega)]^d)}^2 + \frac{1}{\kappa} \|p^\kappa - p^\infty\|_{L^\infty(0,T;L^2(\Omega))}^2 \\ + \eta \|\partial_t u_s^\kappa - \partial_t u_s^\infty\|_{L^2(0,T;[H_0^1(\Omega)]^d)}^2 + \|v_f^\kappa - v_f^\infty\|_{L^2(0,T;[H_0^1(\Omega)]^d)}^2 \lesssim \frac{1}{\kappa^2} \|\partial_t p^\kappa - \partial_t p^\infty\|_{H^1(0,T;L^2(\Omega))}^2. $$

Thus, if $\partial_t p^\infty \in H^1(0, T; L^2(\Omega))$, the above error estimate specifies the convergence speed of $(u_s^\kappa, \partial_t u_s^\kappa, v_f^\kappa, p^\kappa)$ towards $(u_s^\infty, \partial_t u_s^\infty, v_f^\infty, p^\infty)$ as $\kappa$ goes to infinity.

# 5 Numerical experiments

In this section, we present some numerical examples to illustrate the theoretical results presented earlier. In particular, we numerically investigate the regularity of the solutions to the static and time-dependent problems. Note that there is an extensive literature on the numerical approximation of Biot-type systems, see [88, 112, 107, 83, 74, 77, 106, 81, 110, 28, 69, 98] and references therein. The numerical analysis of the specific model (2) presented in this work was performed in [38, 10], where the time discretization is performed with a monolithic backward Euler scheme. Moreover, an alternating minimization splitting scheme was proposed in [27], which leads to a solver closely related to the undrained and fixed-stress splits of Biot's equations. We follow here the monolithic scheme of [38, 10]. In addition, all simulations in this section were performed using the FEniCS finite element software [72, 2].

## 5.1 Spatial discretization

For small values of bulk modulus, our equations can be discretized with standard finite elements. However, when the coefficient $\kappa$ becomes large, we have to take into account the saddle-point structure of the problem involving the mixture's divergence constraint and to choose finite element spaces that satisfy the inf-sup condition (62) at the discrete level.

In the incompressible or nearly incompressible case, the expression of the mapping

$$ T : (u_s, v_s, v_f, p) \mapsto \left( \frac{\alpha^*}{2} u_s - \frac{\alpha^*}{2\lambda_0} v_s, \alpha^* v_s - v_p, \alpha^* v_f - v_p, \alpha^* p \right) $$

defined in (79) suggests us how to select convenient finite element spaces in order to discretize the problem. Indeed, to get a stable discretization, it is sufficient to reproduce the construction of $v_p$ at the discrete level. This is possible by choosing finite elements that are stable (in the Brezzi [25] sense) for Stokes equations.

More precisely, let us suppose that we use a conforming approximation of the space $Y = [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times [H_0^1(\Omega)]^d \times L^2(\Omega)$ by a finite dimensional space

$$ Y_h = V_{s,h} \times V_{s,h} \times V_{f,h} \times Q_h \subset Y, $$

where $V_{s,h}$, $V_{f,h}$ and $Q_h$ denote respectively the finite element spaces chosen to discretize the solid part, the fluid part and the pressure of the mixture. Assume further that $(V_{s,h}, Q_h)$ and $(V_{f,h}, Q_h)$ are two inf-sup stable pairs associated with the standard Stokes problem and verify Fortin Lemma [25, Proposition 5.4.2], i.e. there exists two operators $\Pi_{s,h} : [H_0^1(\Omega)]^d \mapsto V_{s,h}$ and $\Pi_{f,h} : [H_0^1(\Omega)]^d \mapsto V_{f,h}$ satisfying, for each $v \in [H_0^1(\Omega)]^d$,

• For all $q_h \in Q_h$,

$$ \int_{\Omega} \mathrm{div} v q_h \, dx = \int_{\Omega} \mathrm{div} (\Pi_{i,h}(v)) q_h \, dx, \quad i \in \{s, f\}; \qquad (91) $$
• There exists a constant $C_{i,\pi} > 0$ independent of $h$ such that

$$ \|\nabla(\Pi_{i,h}(v))\| \le C_{i,\pi} \|\nabla v\|, \quad i \in \{s, f\}. \qquad (92) $$

Under these hypotheses, we claim that the bilinear form $a_{\lambda_0}$ defined in the mixed formulation (75) is uniformly $T_h$-coercive. Namely, it holds

$$ \forall h > 0, \exists T_h \in \mathcal{L}(Y_h), \forall z_h \in Y_h, \quad |a_{\lambda_0}(z_h, T_h z_h)| \ge \underline{\alpha} \|z_h\|_{Y_h}^2 \quad \text{and} \quad |||T_h||| \le C, \qquad (93) $$

for some constants $\underline{\alpha} > 0$ and $C > 0$ independent of $h$.

Indeed, setting

$$ T_h : (u_{s,h}, v_{s,h}, v_{f,h}, p_h) \mapsto \left( \frac{\alpha^*}{2} u_{s,h} - \frac{\alpha^*}{2\lambda_0} v_{s,h}, \alpha^* v_{s,h} - \Pi_{s,h}(v_{p_h}), \alpha^* v_{f,h} - \Pi_{f,h}(v_{p_h}), \alpha^* p_h \right), \qquad (94) $$

where $v_{p_h}$ is defined by (63), the property (91) enables us to reproduce the calculations from the proof of Proposition 3.15 at the discrete level, so that the first condition of (93) holds true. The second condition then follows from (92) since

$$ \|\nabla(\Pi_{i,h}(v_{p_h}))\| \le C_{i,\pi} \|\nabla v_{p_h}\| \le C_{i,\pi} C_{\text{div}} \|p_h\| $$

for each $i \in \{s, f\}$.

Therefore, a stable discretization of the incompressible system is offered by standard inf-sup stable conforming finite elements associated with the Stokes system. For instance, as it was observed in [10], one can use Taylor-Hood elements $[\mathbb{P}_{k+1}]^d - \mathbb{P}_k$ ($k \ge 1$) for the pairs $(V_{s,h}, Q_h)$ and $(V_{f,h}, Q_h)$. More broadly, the previous $T_h$-coercivity argument implies the stability of the MINI element $\mathbb{P}_1^b - \mathbb{P}_1$, the $\mathbb{P}_2 - \mathbb{P}_0$ element, or also Scott-Vogelius elements $[\mathbb{P}_k]^d - \mathbb{P}_{k-1}^{-1}$ with $k \ge 4$ and $d = 2$. Besides, we can select different finite element spaces for the solid and fluid parts, as long as each of them form a Stokes-stable pair with the space chosen for pressure.

Finally, note that the mapping (94) is independent of the porosity $\phi$ and that the obtention of (93) does not require any assumption on the size of the permeability tensor $k_f$, as it was assumed in [10]. Hence, our approach provides a robust discretization regardless of porosity and permeability.

## 5.2 Regularity of the operator's domain

In both compressible and incompressible cases, we proved the existence and uniqueness of a strong solution in $C^0([0, T]; D(A_{\eta}^{\kappa}))$, with $\eta \ge 0$ and $0 < \kappa \le +\infty$. The operator's domain $D(A_{\eta}^{\kappa})$ was defined by extension from a continuous bilinear form (see e.g. (23) and (24)) but we did not express it as a standard Sobolev space. In what follows, we give some numerical evidences that the operator's domain is not regular, namely

$$ D(A_{\eta \ge 0}^{\kappa}) \neq [\mathbb{H}^2(\Omega)]^d \cap [\mathbb{H}_0^1(\Omega)]^d \times [\mathbb{H}^2(\Omega)]^d \cap [\mathbb{H}_0^1(\Omega)]^d \times [\mathbb{H}^2(\Omega)]^d \cap [\mathbb{H}_0^1(\Omega)]^d \times \mathbb{H}^1(\Omega) $$

and

$$ D(A_{\eta \ge 0}^{\infty}) \neq [\mathbb{H}^2(\Omega)]^d \cap [\mathbb{H}_0^1(\Omega)]^d \times ([\mathbb{H}^2(\Omega)]^d \times [\mathbb{H}^2(\Omega)]^d \cap V_{\phi}) \times \mathbb{H}^1(\Omega). $$

To do so, we compute numerically the solution of the static problem $z + A_{\eta}^{\kappa}z = g$ with $g \in Z$ or $g \in H$, viz.

$$ \left\{ \begin{array}{ll} u_s - v_s = g_u, & \\ \rho_s(1-\phi)v_s - \text{div}(\sigma_s(u_s)) - \text{div}(\sigma_s^{\text{vis}}(v_s)) + (\alpha-\phi)\nabla p = \rho_s(1-\phi)g_s, & \\ \rho_f\phi v_f - \text{div}(\phi\sigma_f(v_f)) + \phi\nabla p = \rho_f\phi g_f, & \\ \frac{\alpha-\phi}{\kappa}p + \text{div}((\alpha-\phi)v_s + \phi v_f) = \begin{cases} \frac{\alpha-\phi}{\kappa}g_p & \text{if } \kappa < +\infty, \\ 0 & \text{if } \kappa = +\infty, \end{cases} & \end{array} \right. \qquad (95) $$

with $(g_u, g_s, g_f, g_p) \in [\mathbb{L}^2(\Omega)]^d \times [\mathbb{L}^2(\Omega)]^d \times [\mathbb{L}^2(\Omega)]^d \times \mathbb{L}^2(\Omega)$ and where we have assumed that $k_f = 0$ and $\theta = 0$ without loss of generality.
We consider $\Omega = \{x \in \mathbb{R}^2, |x| \le 1\}$ a very smooth domain and $(\mathcal{T}_h)_h$ a regular family of meshes of $\bar{\Omega}$, made of triangles. The coarsest mesh size $H$ corresponds to a uniform mesh constructed with 8 subdivisions along each axis direction. Setting $\rho_s = \rho_f = \lambda_f = \mu_f = \eta = \lambda = \mu = 1$ and taking a constant porosity $\phi = 0.5$, we compute the error in $L^2$-norm between the approximated solution $(u_{s,h}, v_{s,h}, v_{f,h}, p_h)$ of (95) and a reference solution computed on a very refined mesh.

The resulting convergence graphs are presented in Figure 1 for $\kappa = 1$ and smooth data $g_s, g_f, g_p$. The convergence rates depend on the regularity of the solid displacement data $g_u$. If $g_u$ is smooth, we obtain optimal orders of convergence, as expected by the theory. However, if $g_u \in [H_0^1(\Omega)]^d \setminus [H^2(\Omega)]^d$, Figure 1 (right) exhibits suboptimal convergence rates, thus indicating that the solution of (95) is not in $H^2$. The same occurs in the incompressible case, as shown in Figure 2.

<figure><img src="image_6.png" /><figcaption>Figure 1: Approximation errors and computed convergence rates for the discretization of the compressible ($\kappa = 1$) steady-state problem (95) with $[P_1]^2 \times [P_1]^2 \times [P_1]^2 \times P_1$ elements.</figcaption></figure>

## 5.3 Regularity of solutions

Even if the operator's domain is not regular, we are now going to shed light on the regularizing effect for the solution of the unsteady problem. For this purpose, we use the time discretization introduced and fully analyzed in [38], which consists of a midpoint scheme for the solid fields and an implicit backward Euler scheme for the fluid and the pressure. The major interest of this scheme is that it preserves energy balance at the discrete level.

As before, we perform the simulation on the smooth domain $\Omega = \{x \in \mathbb{R}^2, |x| \le 1\}$ meshed by a regular family of triangulations $(\mathcal{T}_h)_h$. We set $\rho_s = \rho_f = \lambda_f = \mu_f = \eta = \lambda = \mu = 1$, $k_f = 0$, $\theta = 0$ and $\phi = 0.5$. We take $\kappa = 10^{10}$, but smallest values of $\kappa$ would lead to comparable results. All the simulations are run during a hundred of time iterations, with a time step $\Delta t = 10^{-2}$ and up to the final time $T = 1$.

Figure 3 illustrates the possible regularizing effect of time. Indeed, although the initial velocities and the applied exterior body force belong to $[H_0^1(\Omega)]^d \setminus [H^2(\Omega)]^d$, we recover optimal convergence rates in $L^\infty(0, T; Z)$-norm between the approximated and reference solutions. Indeed, when using $P_1^b \times P_1^b \times P_1^b \times P_1$ elements for spatial discretization, the optimal convergence rate in this norm is equal to 1 for the displacement and to 2 for the other quantities since $Z = [H_0^1(\Omega)]^d \times [L^2(\Omega)]^d \times [L^2(\Omega)]^d \times L^2(\Omega)$. Note that putting $\eta = 0$ slightly degrades the convergence order of the solid velocity, but does not affect the regularity of solid displacement, fluid velocity and pressure.

If the initial conditions and the right-hand side belong strictly to the energy space, Figure 4 highlights
<figure><img src="image_7.png" /><figcaption>Figure 2: Approximation errors and computed convergence rates for the discretization of the incompressible ($\kappa = +\infty$) steady-state problem (95) with $\mathbb{P}_1^b \times \mathbb{P}_1^b \times \mathbb{P}_1^b \times \mathbb{P}_1$ elements.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3: Approximation errors and convergence rates for the discretization of the time-dependent problem (59) with $\mathbb{P}_1^b \times \mathbb{P}_1^b \times \mathbb{P}_1^b \times \mathbb{P}_1$ elements and $v_{s0}, v_{f0}, f$ in $[H_0^1(\Omega)]^d \setminus [H^2(\Omega)]^d$.</figcaption></figure>

that the convergence rates are considerably diminished, which confirms the regularity found in Theorems 3.11 and 3.18. Moreover, the solution is less regular when $\eta = 0$, as predicted by our theoretical results.

## Conclusion

In this work, we study the well-posedness for a fully unsteady and strongly coupled poromechanics model. Using an original combination of semigroup theory and T-coercivity, we demonstrated the existence and uniqueness of strong solutions in the compressible and incompressible cases, with or without solid viscosity. By unifying semigroup and variational techniques, we also recovered the existence and uniqueness of weak
<figure><img src="image_9.png" /><figcaption>Figure 4: Approximation errors and convergence rates for the discretization of the time-dependent problem (59) with $P_1^b \times P_1^b \times P_1^b \times P_1$ elements and $v_{s0}, v_{f0}, f$ in $[L^2(\Omega)]^d \setminus [H_0^1(\Omega)]^d$.</figcaption></figure>

solutions. From the methodological point of view, this unified approach enabled us to take benefit of the best of both theories depending on the hyperbolic-parabolic or parabolic-parabolic nature of the coupling, in particular to prove uniqueness. To handle the incompressible case, we developed a functional framework and an extension of De Rham Theorem adapted to the mixture's divergence constraint. When the incompressible regime is reached, our analysis offers a spatial discretization of the problem with conforming finite elements, including Taylor-Hood elements but also any Stokes-stable elements such as the MINI element. Moreover, this choice of finite element discretization spaces is robust with respect to porosity and permeability. Finally, our theoretical results are corroborated by numerical experiments. We have shown numerically that the operator's domain is not regular, and illustrate the possible regularizing effect on the unsteady problem.

# A Model origin

The model studied in the core of this article results from the linearization of a general poromechanical model presented in [46], where the Biot theory was revisited for finite strains to obtain a non-linear formulation compatible with the principles of thermodynamics, leading to a generic energy balance. In [38] it was linearized for a porous material satisfying Terzaghi's effective stress principle. In what follows, we perform the same linearization without assuming that this principle is satisfied. This leads to a model where the Biot-Willis coefficient may be different from 1, which allows us to relate the resulting linearized model to the standard Biot equations in Section 1.1.

## A.1 The general formulation

We consider a deformable porous material that occupies the space domain $\Omega(t)$ at time $t$. The deformed domain is obtained from a reference configuration domain $\hat{\Omega}$, namely $\Omega(t) = \hat{\mathcal{A}}(\hat{\Omega})$, where

$$ \hat{\mathcal{A}}(\cdot, t) : \hat{\Omega} \to \Omega(t) \\ \hat{x} \mapsto x = \hat{x} + \hat{u}_s(\hat{x}, t) $$

denotes the deformation mapping and $\hat{u}_s$ is the displacement field defined in the reference configuration. We then introduce usual mechanical quantities in the reference configuration such as the deformation gradient
tensor $\hat{F} = \nabla \hat{\mathcal{A}} = \hat{\mathbb{I}} + \nabla \hat{u}_s$, the Cauchy-Green deformation tensor $\hat{\mathcal{C}} = \hat{F}^T \hat{F}$, the Green-Lagrange strain tensor $\hat{E} = \frac{1}{2}(\hat{\mathcal{C}} - \hat{\mathbb{I}})$, and the apparent change of volume of the material $\hat{J} = \det \hat{F}$. By convention, we use a hat superscript for lagrangian quantities (defined in the reference configuration), and no superscript for the corresponding eulerian quantities (defined in the deformed configuration). For instance, $J$ is the function satisfying

$$J(x,t) = J(\hat{\mathcal{A}}(\hat{x},t),t) = \hat{J}(\hat{x},t).$$

The porous material is modeled as a mixture of a fluid phase and a solid phase called the skeleton. At each point of the deformed domain, we assume that the material contains a volume fraction $\phi$ of fluid and $1 - \phi$ of solid, with $\phi(x,t)$ the porosity, and that the fluid and solid parts interact with each other. Moreover, we suppose that the fluid is incompressible, as it is the case in most of biomedical applications.

Let us denote by $\hat{\phi}_0(\hat{x}) = \hat{\phi}(\hat{x},0)$ the initial porosity in the reference configuration, by $\hat{\rho}_s$ the solid density in the reference configuration, by $\rho_f$ the fluid density in the deformed configuration, by $k_f$ the hydraulic conductivity tensor in the deformed configuration – and $\hat{k}_f$ in the reference configuration – which represents the friction between the fluid and solid phases, by $\hat{m} = \rho_f(\hat{J}\hat{\phi} - \hat{\phi}_0)$ the added fluid mass per unit volume in the reference configuration, by $f$ an exterior body force applied to the material and by $\theta$ a distributed fluid mass source term. From [46, Theorem 6], a general form of the fluid Cauchy stress tensor $\sigma_f^{\text{tot}}$ is given by $\sigma_f^{\text{tot}} = \sigma_f(v_f) - p\mathbb{I}$, with

$$\sigma_f(v) = \lambda_f \operatorname{Tr}(\varepsilon(v)) \mathbb{I} + 2\mu_f \varepsilon(v), \quad \varepsilon(v) = \frac{1}{2}(\nabla v + \nabla v^T),$$

and the skeleton contribution to the second Piola-Kirchhoff stress tensor is given by

$$\hat{\Sigma}_s = \frac{\partial \hat{\Psi}_s}{\partial \hat{E}} + \frac{\partial \hat{\Psi}_{\text{damp}}}{\partial \hat{E}} + \hat{\phi} \hat{p} \hat{J} \hat{C}^{-1}, \qquad (96)$$

where $\hat{\Psi}_s$ denotes the skeleton free energy density potential and $\hat{\Psi}_{\text{damp}}$ is a viscous pseudo-potential, here chosen as a function of $\hat{E}$ for simplicity.

The macroscopic governing equations derived in [46] for a general poromechanics formulation valid in large strains and adapted to soft tissue perfusion then read:

$$\left\{
\begin{aligned}
& \hat{\rho}_s(1 - \hat{\phi}) \frac{\partial \hat{v}_s}{\partial t} - \operatorname{div}(\hat{F} \hat{\Sigma}_s(\hat{u}_s, \hat{v}_s, \hat{p})) \\
& \qquad - \hat{J} \hat{\phi}^2 \hat{k}_f^{-1}(\hat{v}_f - \hat{v}_s) + \hat{p} \hat{J} \hat{F}^{-T} \nabla \hat{\phi} = \hat{\rho}_s(1 - \hat{\phi}) \hat{f}, && \forall (\hat{x}, t) \in \hat{\Omega} \times (0, T), \quad (97a) \\
& \frac{1}{J} \frac{\partial}{\partial t} (J \rho_f \phi v_f)|_{\hat{x}} + \operatorname{div}(\rho_f \phi v_f \otimes (v_f - v_s)) - \operatorname{div}(\phi \sigma_f^{\text{tot}}(v_f, p)) \\
& \qquad + \phi^2 \hat{k}_f^{-1}(v_f - v_s) - \theta v_f = \rho_f \phi \hat{f}, && \forall (x, t) \in \Omega(t) \times (0, T), \quad (97b) \\
& \frac{1}{J} \frac{\partial m}{\partial t}|_{\hat{x}} + \operatorname{div}(\rho_f \phi (v_f - v_s)) = \theta, && \forall (x, t) \in \Omega(t) \times (0, T). \quad (97c)
\end{aligned}
\right.$$

In this system, the first two equations correspond respectively to the solid and fluid phases momentum balance equations, and the third one to the mass conservation equation. Finally, noting $\hat{J}_s = (1 - \hat{\phi})\hat{J}$, the system is closed by the relation

$$\hat{p} = - \frac{\partial \hat{\Psi}_s}{\partial \hat{J}_s}. \qquad (98)$$

One fundamental property of such model is its thermodynamical compatibility. Mathematically speaking, this imply that the solution of (97) satisfies an energy balance as proved in [46, Theorem 7]. Assuming in
our case homogeneous Dirichlet conditions for the fluid and the solid, we have

$$ \frac{dK}{dt} + \frac{d\hat{W}_s}{dt} = - \int_{\hat{\Omega}} \frac{\partial \hat{\Psi}_{\text{damp}}}{\partial \hat{E}} : \dot{\hat{E}} \, d\hat{x} - \int_{\Omega(t)} \phi \sigma_f(v_f) : \varepsilon(v_f) \, dx - \int_{\Omega(t)} \phi^2 k_f^{-1} (v_f - v_f) \cdot (v_f - v_s) \, dx \\ + \mathcal{P}_{\text{ext}}^{\text{total}} + \mathcal{J}_{\mathcal{K}\theta} + \mathcal{J}_{\mathcal{G}\theta}, \quad (99) $$

with

$$ K = \frac{1}{2} \int_{\Omega(t)} \rho_s(1-\phi) |v_s|^2 \, dx + \frac{1}{2} \int_{\Omega(t)} \rho_f \phi |v_f|^2 \, dx, \quad \hat{W}_s = \int_{\hat{\Omega}} \hat{\Psi}_s \, d\hat{x} $$

the mixture's kinetic energy and the skeleton free energy,

$$ \mathcal{J}_{\mathcal{K}\theta} = \frac{1}{2} \int_{\Omega(t)} \theta |v_f|^2 \, dx, \quad \mathcal{J}_{\mathcal{G}\theta} = \int_{\Omega(t)} \frac{p}{\rho_f} \theta \, dx $$

the incoming rates of fluid kinetic energy and Gibbs free energy, and

$$ \mathcal{P}_{\text{ext}}^{\text{total}} = \int_{\Omega(t)} \rho_f \phi f \cdot v_r \, dx + \int_{\Omega(t)} \rho_f f \cdot v_s \, dx = \int_{\Omega(t)} \rho_f \phi f \cdot v_f \, dx + \int_{\Omega(t)} \rho_s (1-\phi) f \cdot v_s \, dx $$

the power of external forces, with $v_r = v_f - v_s$ the relative velocity between the fluid and the solid and $\rho = \rho_s(1-\phi) + \rho_f\phi$ the porous material total density. One benefit of considering the linearized version of (97) is to keep such energy balance after linearization.

But before linearizing the coupled system, let us specify the skeleton constitutive law. Following the guidelines of [46, Section 5.4], we consider a free energy density potential of the form

$$ \hat{\Psi}_s = \hat{W}^{\text{skel}}(\hat{E}) + \hat{W}^{\text{bulk}}\left(\hat{J}_s \frac{1-\hat{\phi}_0}{\hat{\chi}_s(\hat{J})}\right), \quad (100) $$

where $\hat{\chi}_s(\hat{J})$ is a function representing the variations of solid volume directly due to macroscopic volume changes in the absence of pore pressure. In other words, (100) means that the energy bulk term depends on the ratio between the change of volume for the solid part $\hat{J}_s$ and the change of volume $\hat{\chi}(\hat{J}) = \frac{\hat{\chi}_s(\hat{J})}{1-\hat{\phi}_0}$ occurring in each cell of the microstructure when assuming that the pore pressure is constant during the deformation. We suppose that $\hat{\chi}(\hat{J})$ is affine with respect to $\hat{J}$ and that $\hat{\chi}(\hat{J}) = 1$ when $\hat{J} = 1$, namely

$$ \hat{\chi}(\hat{J}) - 1 = \beta(\hat{J} - 1) \quad \Leftrightarrow \quad \hat{\chi}_s(\hat{J}) = (1 - \hat{\phi}_0)(1 + \beta(\hat{J} - 1)), $$

with $0 \le \beta \le 1$ a coefficient vanishing for incompressible materials.

Combining (98) and (100), we see that $\hat{p} = -\frac{\partial \hat{W}^{\text{bulk}}}{\partial \hat{J}_s}$. From (96), it follows that

$$ \hat{\Sigma}_s = \frac{\partial \hat{W}^{\text{skel}}}{\partial \hat{E}} - \hat{p} \frac{\partial \hat{J}_s}{\partial \hat{E}} + \frac{\partial \hat{W}^{\text{bulk}}}{\partial \hat{J}} \cdot \frac{\partial \hat{J}}{\partial \hat{E}} + \frac{\partial \hat{\Psi}_{\text{damp}}}{\partial \hat{E}} + \hat{\phi} \hat{p} \hat{J} \hat{C}^{-1}. $$

Observing that $\frac{\partial}{\partial \hat{E}}(\hat{J}\hat{\phi}) = \frac{\partial}{\partial \hat{E}}(\frac{\hat{m}}{\rho_f} + \hat{\phi}_0) = 0$ because $\hat{m}$ corresponds to the fluid mass entering into the pores in the reference configuration and thus does not depend on the deformation of the material, we get $\frac{\partial \hat{J}_s}{\partial \hat{E}} = \frac{\partial \hat{J}}{\partial \hat{E}} = \hat{J} \hat{C}^{-1}$. Since $\frac{\partial \hat{W}^{\text{bulk}}}{\partial \hat{J}} = \frac{\hat{J}_s \hat{\chi}'_s(\hat{J})}{\hat{\chi}_s(\hat{J})}$, we obtain

$$ \hat{\Sigma}_s = \frac{\partial \hat{W}^{\text{skel}}}{\partial \hat{E}} + \frac{\partial \hat{\Psi}_{\text{damp}}}{\partial \hat{E}} - \left(1 - \frac{\hat{J}_s \hat{\chi}'_s(\hat{J})}{\hat{\chi}_s(\hat{J})} - \hat{\phi}\right) \hat{p} \hat{J} \hat{C}^{-1}. \quad (101) $$
Note that this expression is valid for any potentials $\hat{W}^{skel}$ and $\hat{\Psi}_{damp}$. To further develop the computations, we can for example use a Ciarlet-Geymonat-like potential [49] for the bulk potential $\hat{W}^{bulk}$, which yields

$$ \hat{W}^{bulk} \left( \hat{J}_s \frac{1 - \hat{\phi}_0}{\hat{\chi}_s(\hat{J})} \right) = \hat{\gamma} \kappa \left( \frac{\hat{J}_s}{\hat{\chi}_s(\hat{J})} - 1 - \log \left( \frac{\hat{J}_s}{\hat{\chi}_s(\hat{J})} \right) \right), $$

where $\kappa$ denotes the solid grains bulk modulus and $\hat{\gamma}$ is a scaling factor. In order to recognize the storage coefficient in the pressure equation, we choose from now on $\hat{\gamma} = \frac{1-\hat{\phi}}{1-\beta}$. Using (98), we get

$$ \hat{p} = - \frac{(1 - \hat{\phi})\kappa}{1 - \beta} \left( \frac{1}{\hat{\chi}_s(\hat{J})} - \frac{1}{\hat{J}_s} \right) = \frac{(1 - \hat{\phi})\kappa}{1 - \beta} \cdot \frac{\hat{\chi}_s(\hat{J}) - \hat{J}_s}{\hat{J}_s \hat{\chi}_s(\hat{J})}. $$

Remarking that $\hat{\chi}_s(\hat{J}) - \hat{J}_s = \hat{J}\hat{\phi} - \hat{\phi}_0 + \beta(1 - \hat{\phi}_0)(\hat{J} - 1) + 1 - \hat{J} = \rho_f^{-1}\hat{m} - (1 - \beta(1 - \hat{\phi}_0))(\hat{J} - 1)$, we obtain

$$ \hat{p} = \frac{\kappa}{(1 - \beta)(1 - \hat{\phi}_0)} \cdot \frac{\rho_f^{-1} \hat{m} - (1 - \beta(1 - \hat{\phi}_0))(\hat{J} - 1)}{\hat{J}(1 + \beta(\hat{J} - 1))}. \qquad (102) $$

This closure relation will be the cornerstone of the linearization process. As a matter of fact, it involves the interstitial pressure $\hat{p}$, the fluid added mass $\hat{m}$ that is related to the porosity $\hat{\phi}$, and the change of volume $\hat{J}$ that is close to 1 when linearizing the coupled system.

## A.2 Linearization

We linearize (97) for infinitesimal transformations around the configuration $(\hat{u}_s, \hat{v}_s, v_f, \hat{\phi}) = (0, 0, 0, \hat{\phi}_0)$. This first amounts to assume that

$$ \hat{A} = \hat{I} + \mathcal{O}(|\nabla \hat{u}_s|), \quad \hat{E} = \hat{\varepsilon} + \mathcal{O}(|\nabla \hat{u}_s|^2), $$

where $|\nabla \hat{u}_s|^2 = \nabla \hat{u}_s : \nabla \hat{u}_s$. Thus, the reference and deformed configurations reduce to a single domain, which will be denoted by $\Omega$, allowing us to drop from now on the hat superscripts used previously to distinguish variables defined on $\hat{\Omega}$ or $\Omega(t)$. Besides, it holds that

$$ J = 1 + \mathrm{Tr} \, \varepsilon + \mathcal{O}(|\varepsilon|^2) = 1 + \mathrm{div} \, u_s + \mathcal{O}(|\nabla u_s|^2) $$

and, in virtue of (102), that

$$ \phi = \phi_0 + \mathcal{O}(|(\mathrm{div} \, u_s, p)|). $$

Furthermore, any choice of potentials $W^{skel}(E)$ and $\Psi_{damp}(\hat{E})$ satisfies

$$ \frac{\partial W^{skel}}{\partial E} = \lambda \mathrm{Tr}(\varepsilon) \mathbb{I} + 2\mu \varepsilon + \mathcal{O}(|\nabla u_s|^2) \quad \text{and} \quad \frac{\partial \Psi_{damp}}{\partial \hat{E}} = \nu \mathrm{Tr}(\dot{\varepsilon}) \mathbb{I} + 2\eta \dot{\varepsilon} + \mathcal{O}(|\nabla v_s|^2) $$

for some Lamé constants $\lambda, \mu, \nu$ and $\eta$, where $\eta$ represents the solid grains viscosity. To simplify, we suppose that $\nu = 0$, but note that choosing $\nu = \lambda^* > 0$ would mean taking into account secondary consolidation effects as in [78].

Since $\frac{J_s \chi'_s(J)}{\chi_s(J)} = \frac{J_s \beta(1-\phi_0)}{\chi_s(J)}$ and $\frac{J_s}{\chi_s(J)} = 1 + \mathcal{O}(|(\mathrm{div} \, u_s, p)|)$, (101) implies that

$$ \sigma_s = \lambda \mathrm{Tr}(\varepsilon) \mathbb{I} + 2\mu \varepsilon + 2\eta \dot{\varepsilon} - (1 - \beta(1 - \phi_0) - \phi_0) p \mathbb{I} + \mathcal{O}(|(\nabla u_s, \nabla v_s, p)|^2). $$

The Biot-Willis coefficient, denoted by $\alpha$, is then defined as the coefficient multiplying the pressure term in the porous material linearized total stress tensor $\Sigma = \sigma_s + \phi_0 \sigma_f^{\mathrm{tot}}$. Therefore

$$ \alpha = 1 - \beta(1 - \phi_0), \qquad (103) $$
and (97a) becomes

$$ \rho_s(1 - \phi_0) \partial_t v_s - \operatorname{div} (\lambda \operatorname{Tr}(\varepsilon(u_s)) \mathbb{I} + 2\mu \varepsilon(u_s) + 2\eta \varepsilon(v_s)) \\ - \phi_0^2 k_f^{-1} (v_f - v_s) + (\alpha - \phi_0) \nabla p = \rho_s(1 - \phi_0)f, \quad \forall (x,t) \in \Omega \times (0,T), $$

which corresponds exactly to (2a).

The fluid equation (2b), namely

$$ \rho_f \phi \partial_t v_f - \operatorname{div}(\phi \sigma_f(v_f)) + \phi^2 k_f^{-1} (v_f - \partial_t u_s) - \theta v_f + \phi \nabla p = \rho_f \phi f, \quad \forall (x,t) \in \Omega \times (0,T), $$

readily results from the linearization of (97b). To recover (2c), we infer from (102) that

$$ \frac{(1 - \beta)(1 - \phi_0)}{\kappa} p = \frac{\rho_f^{-1} m - \alpha \operatorname{div} u_s + \mathcal{O}(|\nabla u_s|^2)}{1 + \mathcal{O}(|\operatorname{div} u_s|)}. $$

Since $(1 - \beta)(1 - \phi_0) = \alpha - \phi_0$, it follows that

$$ \rho_f^{-1} m = \frac{\alpha - \phi_0}{\kappa} p + \alpha \operatorname{div} u_s + \mathcal{O}(|(\nabla u_s, p)|^2). $$

Differentiating this relation with respect to time and using (97c), we finally get

$$ \frac{\alpha - \phi_0}{\kappa} \partial_t p + \operatorname{div}((\alpha - \phi_0)v_s + \phi_0 v_f) = \rho_f^{-1} \theta, \quad \forall (x,t) \in \Omega \times (0,T), $$

which coincides with the pressure equation (2c) studied in the core of the article.

**Remark A.1.** The coefficient $\frac{\alpha - \phi_0}{\kappa}$ is known as the storage coefficient and is often denoted by $c_0$. The Biot-Willis coefficient, given by (103), is usually computed by the formula $\alpha = 1 - \frac{\kappa_0}{\kappa}$, with $\kappa_0$ the drained bulk modulus. Hence $\kappa_0 = \beta(1 - \phi_0)\kappa$ and we have in particular

$$ \kappa \to +\infty \Leftrightarrow \alpha \to 1 \Leftrightarrow \beta \to 0 \Leftrightarrow \chi_s(J) \to 1 - \phi_0, $$

so that (100) reduces to

$$ \Psi_s = W^{\text{skel}}(E) + W^{\text{bulk}}(J_s) $$

for nearly-incompressible materials. Such a decomposition for $\Psi_s$ corresponds to a material satisfying Terzagi's effective stress principle, for which the microstructure Poisson effects are not taken into account.

As previously mentioned, let us observe that the non-linear energy balance (99) also holds true after linearization. Indeed, replacing each quantity of (99) by its linearized counterpart, we recover exactly (10).

## Acknowledgement

The authors of this paper would like to express their gratitude to Dominique Chapelle (Inria) for his advices on modeling aspects related to Appendix A. Moreover, we are thankful to Patrick Ciarlet (ENSTA Paris) for insightful discussions about the use of T-coercivity to solve Stokes equations at the continuous and discrete levels. Finally, we express our gratitude to the reviewers for their valuable comments that allowed us to improve in particular the positioning of this work within the literature on poroelasticity and to clarify some steps.
References

[1] C. Ager, B. Schott, M. Winter, and W. A. Wall. A Nitsche-based cut finite element method for the coupling of incompressible fluid flow with poroelasticity. *Computer Methods in Applied Mechanics and Engineering*, 351:253-280, 2019.

[2] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The FEniCS project version 1.5. *Archive of Numerical Software*, 3(100):9-23, 2015.

[3] I. Ambartsumyan, E. Khattatov, I. Yotov, and P. Zunino. A Lagrange multiplier method for a Stokes-Biot fluid-poroelastic structure interaction model. *Numerische Mathematik*, 140(2):513-553, 2018.

[4] D. Ambrosi and L. Preziosi. On the closure of mass balance models for tumor growth. *Mathematical Models and Methods in Applied Sciences*, 12(05):737-754, 2002.

[5] P. Angot. Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions. *ESAIM: Mathematical Modelling and Numerical Analysis*, 52(5):1875-1911, 2018.

[6] J. L. Auriault. Dynamic behaviour of a porous medium saturated by a newtonian fluid. *International Journal of Engineering Science*, 18(6):775-785, Jan. 1980.

[7] J.-L. Auriault. Poroelastic media. In *Homogenization and porous media*, pages 163-182. Springer, 1997.

[8] G. Avalos and R. Triggiani. Semigroup well-posedness in the energy space of a parabolic-hyperbolic coupled Stokes-Lamé PDE system of fluid-structure interaction. *Discrete & Continuous Dynamical Systems - S*, 2(3):417-447, 2009.

[9] S. Badia, A. Quaini, and A. Quarteroni. Coupling Biot and Navier-Stokes equations for modelling fluid-poroelastic media interaction. *Journal of Computational Physics*, 228(21):7986-8014, 2009.

[10] N. Barnafi, P. Zunino, L. Dedè, and A. Quarteroni. Mathematical analysis and numerical approximation of a general linearized poro-hyperelastic model. *Computers & Mathematics with Applications*, 91:202-228, 2021.

[11] M. Barré and P. Ciarlet Jr. The T-coercivity approach for mixed problems. Submitted, 2022.

[12] H. Barucq, M. Madaune-Tort, and P. Saint-Macary. Theoretical aspects of wave propagation for Biot's consolidation problem. *Monografías del Seminario Mathemático García de Galdeano*, 31:449-458, 2004.

[13] P. J. Basser. Interstitial pressure, volume, and flow during infusion into brain tissue. *Microvascular Research*, 44(2):143-165, Sept. 1992.

[14] A. Bensoussan, G. Da Prato, M. C. Delfour, and S. K. Mitter. *Representation and control of infinite dimensional systems*. Springer Science & Business Media, 2007.

[15] L. Berger, R. Bordas, K. Burrowes, V. Grau, S. Tavener, and D. Kay. A poroelastic model coupled to a fluid network with applications in lung modelling. *International Journal for Numerical Methods in Biomedical Engineering*, 32(1), 2016.

[16] M. A. Biot. General theory of three-dimensional consolidation. *Journal of applied physics*, 12(2):155-164, 1941.

[17] M. A. Biot. Theory of elasticity and consolidation for a porous anisotropic solid. *Journal of applied physics*, 26(2):182-185, 1955.
[18] M. A. Biot. Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low-Frequency Range. *The Journal of the Acoustical Society of America*, 28(2):168-178, Mar. 1956.

[19] M. A. Biot and G. Temple. Theory of finite deformations of porous solids. *Indiana University Mathematics Journal*, 21(7):597-620, 1972.

[20] L. Bociu, S. Canic, B. Muha, and J. T. Webster. Multilayered Poroelasticity Interacting with Stokes Flow. *SIAM Journal on Mathematical Analysis*, 53(6):6243-6279, Jan. 2021.

[21] L. Bociu, G. Guidoboni, R. Sacco, and M. Verri. On the role of compressibility in poroviscoelastic models. *Mathematical Biosciences and Engineering*, 16(5):6167-6028, 2019.

[22] L. Bociu, G. Guidoboni, R. Sacco, and J. T. Webster. Analysis of Nonlinear Poro-Elastic and Poro-Visco-Elastic Models. *Archive for Rational Mechanics and Analysis*, 222(3):1445-1519, Dec. 2016.

[23] L. Bociu, B. Muha, and J. T. Webster. Weak solutions in nonlinear poroelasticity with incompressible constituents. *Nonlinear Analysis: Real World Applications*, 67:103563, Oct. 2022.

[24] L. Bociu and J. T. Webster. Nonlinear quasi-static poroelasticity. *Journal of Differential Equations*, 296:242-278, Sept. 2021.

[25] D. Boffi, F. Brezzi, M. Fortin, et al. *Mixed finite element methods and applications*, volume 44. Springer, 2013.

[26] M. E. Bogovskii. Solution of the first boundary value problem for the equation of continuity of an incompressible medium. In *Doklady Akademii Nauk*, volume 248, pages 1037-1040. Russian Academy of Sciences, 1979.

[27] J. W. Both, N. A. Barnafi, F. A. Radu, P. Zunino, and A. Quarteroni. Iterative splitting schemes for a soft material poromechanics model. *Computer Methods in Applied Mechanics and Engineering*, page 29, 2022.

[28] J. W. Both, M. Borregales, J. M. Nordbotten, K. Kumar, and F. A. Radu. Robust fixed stress splitting for Biot's equations in heterogeneous media. *Applied Mathematics Letters*, 68:101-108, June 2017.

[29] J. W. Both, I. S. Pop, and I. Yotov. Global existence of weak solutions to unsaturated poroelasticity. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55(6):2849-2897, 2021.

[30] E. Bourgeois. Mécanique des milieux poreux en transformation finie: formulation des problèmes et méthodes de résolution. *Thèse de doctorat, École Nationale des Ponts et Chaussées*, 1997.

[31] R. M. Bowen. Incompressible porous media models by use of the theory of mixtures. *International Journal of Engineering Science*, 18(9):1129-1148, 1980.

[32] F. Boyer and P. Fabrie. *Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models*, volume 183. Springer Science & Business Media, 2012.

[33] M. Bukač, I. Yotov, R. Zakerzadeh, and P. Zunino. Partitioning strategies for the interaction of a fluid with a poroelastic material based on a Nitsche's coupling approach. *Comput. Methods Appl. Mech. Engrg.*, page 33, 2015.

[34] M. Bukač, I. Yotov, and P. Zunino. An operator splitting approach for the interaction between a fluid and a multilayered poroelastic structure. *Numerical Methods for Partial Differential Equations*, 31(4):1054-1100, July 2015.

[35] M. Bukač, I. Yotov, and P. Zunino. Dimensional model reduction for flow through fractures in poroelastic media. *ESAIM: Mathematical Modelling and Numerical Analysis*, Nov. 2016.
[36] N. Burq and P. Gérard. *Contrôle optimal des équations aux dérivées partielles*. École Polytechnique, Département de mathématiques, 2002.

[37] B. Burtschell, D. Chapelle, and P. Moireau. Effective and energy-preserving time discretization for a general nonlinear poromechanical formulation. *Computers & Structures*, 182:313–324, 2017.

[38] B. Burtschell, P. Moireau, and D. Chapelle. Numerical analysis for an energy-stable total discretization of a poromechanics model with inf-sup stability. *Acta Mathematicae Applicatae Sinica, English Series*, 35(1):28–53, 2019.

[39] V. Calo, N. Brasher, Y. Bazilevs, and T. Hughes. Multiphysics model for blood flow and drug transport with application to patient-specific coronary artery flow. *Computational Mechanics*, 43(1):161–177, 2008.

[40] S. Čanić, Y. Wang, and M. Bukač. A Next-Generation Mathematical Model for Drug-Eluting Stents. *SIAM Journal on Applied Mathematics*, 81(4):1503–1529, Jan. 2021.

[41] Y. Cao, S. Chen, and A. Meir. Analysis and numerical approximations of equations of nonlinear poroelasticity. *Discrete & Continuous Dynamical Systems-B*, 18(5):1253, 2013.

[42] P. Causin, G. Guidoboni, A. Harris, D. Prada, R. Sacco, and S. Terragni. A poroelastic model for the perfusion of the lamina cribrosa in the optic nerve head. *Mathematical Biosciences*, 257:33–41, Nov. 2014.

[43] A. Cesmelioglu. Analysis of the coupled Navier-Stokes/Biot problem. *Journal of Mathematical Analysis and Applications*, 456(2):970–991, Dec. 2017.

[44] R. Chabiniok, V. Y. Wang, M. Hadjicharalambous, L. Asner, J. Lee, M. Sermesant, E. Kuhl, A. A. Young, P. Moireau, M. P. Nash, et al. Multiphysics and multiscale modelling, data-model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics. *Interface focus*, 6(2):20150083, 2016.

[45] D. Chapelle, J.-F. Gerbeau, J. Sainte-Marie, and I. Vignon-Clementel. A poroelastic model valid in large strains with applications to perfusion in cardiac modeling. *Computational Mechanics*, 46(1):91–101, 2010.

[46] D. Chapelle and P. Moireau. General coupling of porous flows and hyperelastic formulations—from thermodynamics principles to energy balance and compatible time schemes. *European Journal of Mechanics-B/Fluids*, 46:82–96, 2014.

[47] L. Chesnel and P. Ciarlet. T-coercivity and continuous Galerkin methods: application to transmission problems with sign changing coefficients. *Numerische Mathematik*, 124(1):1–29, 2013.

[48] D. Chou, J. C. Vardakis, L. Guo, B. J. Tully, and Y. Ventikos. A fully dynamic multi-compartmental poroelastic system: Application to aqueductal stenosis. *Journal of Biomechanics*, 49(11):2306–2312, July 2016.

[49] P. G. Ciarlet. *Mathematical Elasticity: Volume I: three-dimensional elasticity*. North-Holland, 1988.

[50] P. Ciarlet Jr. T-coercivity: Application to the discretization of Helmholtz-like problems. *Computers & Mathematics with Applications*, 64(1):22–34, 2012.

[51] O. Coussy. *Poromechanics*. John Wiley & Sons, 2004.

[52] C. M. Dafermos. On the existence and the asymptotic stability of solutions to the equations of linear thermoelasticity. *Archive For Rational Mechanics And Analysis*, 29(4):241–271, Jan. 1968.