# Incremental variational approach to gradient damage coupled with poroelasticity of saturated media

Xiao-Dong Zhang, Long Cheng, Djimédo Kondo, Albert Giraud

► To cite this version:

Xiao-Dong Zhang, Long Cheng, Djimédo Kondo, Albert Giraud. Incremental variational approach to gradient damage coupled with poroelasticity of saturated media. Journal of the Mechanics and Physics of Solids, 2024, 187, pp.105614. 10.1016/j.jmps.2024.105614. hal-05512111

HAL Id: hal-05512111

https://hal.science/hal-05512111v1

Submitted on 27 Feb 2026

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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Incremental variational approach to gradient damage coupled with
poroelasticity of saturated media

Xiao-Dong Zhang<sup>a</sup>, Long Cheng<sup>a,*</sup>, Djimédo Kondo<sup>b</sup>, Albert Giraud<sup>a</sup>

<sup>a</sup>GeoRessources Laboratory, Université de Lorraine (ENSG), CNRS UMR7359, F54518, Vandoeuvre-lès-Nancy, France

<sup>b</sup>Institut Jean Le Rond d'Alembert, CNRS UMR7190, Sorbonne Université, 4 place Jussieu, F75005, Paris, France

**Abstract**

In this study, we aim at investigating the coupling between poroelasticity (including the fluid flow) and gradient damage phenomena in saturated porous media. To this end, we first extend the thermodynamics-based Biot-Coussy theory of poroelasticity in order to incorporate gradient damage processes. Taking advantage of this framework, we establish a variational formulation for the proposed model, expressed as a four-field incremental minimization problem of a total energy functional which includes poroelastic energy as well as dissipated energy related to damage growth together with Darcy-like fluid flow. Moreover, after a careful analysis of the so-called kinetic porosity in each increment, the incremental minimization is reduced to that of three fields, namely, the displacement and damage fields of the skeleton phase of the porous media as well as the pore fluid pressure. Subsequently, we proceed to the numerical implementation of the variational problem via a semi-staggered optimization algorithm and apply it to a first benchmark modeling for which corresponding solution as well as numerical results for hydraulic fracturing are available. Finally, we perform an application of the proposed model for the evaluation of an Excavation Damage Zone (EDZ) around an underground waste storage gallery. The model is shown to be able to deliver sufficiently reliable predictions.

**Keywords:** Gradient damage model, Poroelasticity, Thermodynamics formulation, Incremental variational approach, Darcy flow, Phase field method, Simulation by FEM

**1. Introduction**

During the last decades, a large number of studies have been devoted to Continuum Damage Mechanics (CDM) and its ability to predict a transition to the failure of various dissipative structures. Following the pioneering work of Kachanov (1958) and the Lemaitre and Chaboche (1978), CDM has been reformulated (see for instance Marigo (1981)) in the framework of Generalized Standard Materials (GSM) that has been previously introduced by Halphen and Nguyen (1975) in the line of Ziegler (1963); Coleman and Gurtin (1967); Rice (1971); Moreau (1974). It was at the end of 1980s, and more markedly in the 1990s, that the question of regularized damage modeling received particular attentions via various non-local approaches. The latter aims at leading

<footnote>*Corresponding author</footnote>
to robust numerical tools for a better description of the transition from damage to fracture. The
most probably important developments may be found in those accounting for the non-local integral
formulation (Bazant and Pijaudier-Cabot, 1988), the higher-order deformation gradients approach
(Peerlings et al., 1996) as well as gradient damage models (Frémond and Nedjar, 1993; Frémond
and Nedjar, 1996; Lorentz and Andrieux, 1999; Lorentz and Benallal, 2005)<fnref n="1" />, etc.. The latter one,
which proceeds to incorporate the gradient of damage in the modeling, is adopted for the present
study, not only for its simplicity, but also by its potential links with recent variational approaches
of fracture mechanics (see for instance Pham et al. (2011); Marigo et al. (2016)). Notably, these
studies are related to the work of Bourdin et al. (2000) in which the variational regularization ap-
proach results in a minimization of a two-fields (displacement field and a scalar phase field) energy
functional. This approach is known nowadays as the variational Phase Field method which allows
to efficiently predict fracture nucleation and growth phenomena. Another notable contribution
is the consideration of unilateral contact effect of damage under compressive loadings by Amor
et al. (2009), who appropriately split the strain energy in a spherical and a deviatoric parts. The
number of applications of the above variational phase-field approach to the mechanical responses
of different engineering structures are numerous and will not be referred to here. Some of them
can be found in Bourdin et al. (2008); Wu et al. (2020) and in Martínez-Pañeda et al. (2018) (for
hydrogen assisted cracking).

The above-mentioned studies mainly focus on the coupling between elasticity and damage
phenomena. Despite their interest, they are not directly applicable to several civil engineering or
petroleum engineering problems that primarily involve hydromechanical couplings. This can be
addressed by relying on poroelasticity theory which goes back to the pioneering fundamental work
of Biot (1941) (see also Biot (1955)). The Biot's poroelasticity theory has been later developed
by Rice and Cleary (1976) and reformulated in Coussy (1995) (see also Coussy (2004)) by taking
advantage of thermodynamics of irreversible processes. One can also refer to the textbooks of Cheng
(2016) or for microporomechanics to Dormieux et al. (2006), and also to Dormieux et al. (2002).
Among others, recent applications to hydraulic fracturing can be found in Gao and Detournay
(2020).

The extension of variational phase field method to hydraulic fractures problem in permeable
rocks has been done by Yoshioka and Bourdin (2016) and applied in several studies (see for instance
Chukwudozie et al. (2019); Mollaali et al. (2019)). In these researches, a decoupling has been
assumed between the poroelasticity and fluid transport. This consists first in the determination
of fluid pressure field by solving the mass balance equation together with that of Darcy law for
fluid flow. Then, knowing the pressure field, the variational hydraulic fracturing problem becomes
that of the minimization of two fields (i.e. displacement and damage) functional. Concerning the
coupled poroelasticity-damage, another interesting study was that proposed by Miehe et al. (2015)
(see also Mauthe and Miehe (2017)). This delivered a related variational principle for continuous
evolution problem as well as its incremental form. The approach proposed by these authors was
based on the Darcy-Biot-type fluid transport, in which the fluid mass content constitutes one of the
state variables. However, in these studies, the poroelastic parameters (e.g. Biot coefficient and Biot
modulus) are not affected by the damage. Yet, all the above mentioned works have concentrated
on the objectives of hydraulic fracturing modeling but unfortunately have not demonstrated their
efficiency in the context of poromechanical damage, known as of great importance in geomechanics

<footnote n="1">¹It is interesting to note that the work of Lorentz and Andrieux (1999) is one of the first which proceeds by an incremental variational approach of the coupling between elasticity and damage.</footnote>
and geotechnics. By defining the so-called normalized fluid content variation, Ulloa et al. (2022) recently extended the work of Miehe et al. (2015) and proposed a rate form of saddle point principle to incorporate some microporomechanics considerations and non-associative frictional plasticity law. Moreover, their treatment of the coupling between the poromechanics and the fluid transport was carried out with a staggered procedure by separately solving the elastoplastic evolution, fluid balance and damage evolution equations, where the solution of the former one is obtained by linearizing the mechanical balance equation with plastic strains updated via local return-mapping.

We aim in the present work at developing an incremental variational approach to the gradient damage model for saturated poroelastic media. By doing so, the solution of the fully coupled problem will be obtained from the minimization of an incremental energy functional. More specifically, after a brief recall of some elements of thermodynamics for saturated porous media, we propose in the context of generalized standard materials an extension which allows to couple poroelasticity with gradient damage (see Section 2). Based on such framework, an incremental variational principle is then developed in Section 3 as dependent on four fields, which are the displacement and damage fields of the skeleton, the fluid pressure and the so-called kinetic porosity. The latter one is then carefully approximated in each increment that allows to reduce and express the incremental variational principle as based on functional of the three first fields. This is further applied in Section 4 via a benchmark modeling by assuming the fluid transport obeying the linear Darcy flow. Next, by making use of a semi-staggered optimization algorithm, we assess and validate the proposed model in Section 5 by comparing its predictions with some available results of hydraulic fracturing problem. This allows to demonstrate the efficiency and accuracy of the poroelastic gradient damage model. Finally, a large scale structural application has been performed in order to investigate the occurrence and evolution of an excavation damage zone (EDZ) around an underground opening. Pore pressure field as well as fracture pattern are also shown.

## 2. Thermodynamics of poroelasticity coupled with gradient damage

Thermodynamics provides a theoretical framework allowing to appropriately describe the evolution of a mechanical system depending on all forms of energy. In this section, based on the basic elements of continuous saturated media established by Coussy (2004), we formulate a suitable thermodynamic description of poroelasticity with fluid flow coupled to gradient damage. To this end, let us first consider a saturated porous medium, as shown in Fig.1, denoted by $\Omega$, in which each material point (i.e. infinitesimal elementary volume $d\Omega$) is considered to be made of a superposition of solid and fluid phases. Specifically, the skeleton is made up of an infinitesimally deformable solid matrix and connected pore space, whose volume fraction constitutes the material porosity in the present study<fnref n="2" />. This connected porous space is assumed to be saturated by an incompressible fluid continuum.

### 2.1. Brief recall of thermodynamic principles for poromechanics

As well known, the first principle of the thermodynamics states the conservation of energy and, the second one corresponds to the positivity of energy dissipation (this allows to distinguish

<footnote n="2">²If there exist occluded pores in the skeleton, saturated or not, they are considered as totally segregated from the connected ones without any fluid filtration. The corresponding porosity, i.e. the so-called occluded porosity, and its effects on the skeleton will be taken into account in the behaviors of the solid matrix by default.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Saturated porous media and its infinitesimal volume</figcaption></figure>

between the reversible and irreversible transformations and those of the irreversible dissipation of the concerned system). As classically, they are respectively expressed as:

$$ \dot{\mathcal{E}} = -\mathcal{P}_{int} + \mathcal{P}_{cal} \qquad (1) $$

for the first principle and

$$ \dot{\mathcal{S}} \geq \int_{\Omega} \frac{r}{T} d\Omega - \int_{\partial\Omega} \frac{\mathbf{q} \cdot \mathbf{n}}{T} dS \qquad (2) $$

for the second one.

$\mathcal{E}$ is the global internal energy of the concerned system; $\mathcal{P}_{int}$ gives its power of internal work including that of possible microforces; $\mathcal{P}_{cal}$ symbolizes the calorific power received by the system, which is furnished by the external heat supplies per unit volume, denoted by $r$, and per unit surface, denoted by $\mathbf{q}$.

$S$ denotes the global entropy; $\mathbf{n}$ is the unit outward normal on the boundary of $\partial\Omega$ and $T$ the local absolute temperature.

Coming to poromechanics, let us suppose, as aforementioned, the infinitesimal transformation for the skeleton, and neglect the gravity as well as inertia forces of $d\Omega$. The two principles can be written into their local forms, which read (see for instance Coussy (2004)):

$$ \dot{\mathcal{e}} = \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} + r - \operatorname{div}(\mathbf{q}) - \operatorname{div} \left[ \rho^f \mathbf{v} \left( e^f + \frac{p}{\rho^f} \right) \right] \qquad (3) $$

and

$$ \mathcal{D} = \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} - \dot{\mathcal{e}} + T\dot{s} + T\operatorname{div}(s^f \rho^f \mathbf{v}) - \operatorname{div} \left[ \rho^f \mathbf{v} \left( e^f + \frac{p}{\rho^f} \right) \right] - \frac{\nabla T}{T} \cdot \mathbf{q} \geq 0 \qquad (4) $$

Specifically, in Eq.(3), $e$ is the overall density of internal energy per unit volume (and not per unit mass) of $d\Omega$, which can be calculated through $e = \rho^s (1 - \phi) e^s + \rho^f \phi e^f$ with $\phi$ the Lagrangian porosity, $e^s$ and $e^f$ denoting the specific internal energy of the solid and fluid parts, $\rho^s$ and $\rho^f$ their mass density, respectively; $p$ is the pore fluid pressure; $\boldsymbol{\sigma}$ is the Cauchy stress tensor; $\boldsymbol{\varepsilon} = \nabla_s \mathbf{u}$ gives the linearized strain tensor of the skeleton with $\mathbf{u}$ the displacement vector and $\nabla_s$ the symmetric gradient operator; $\mathbf{v} = \phi(\mathbf{v}^f - \mathbf{v}^s)$ furnishes the filtration velocity of the fluid with respect to the skeleton with $\mathbf{v}^f$ and $\mathbf{v}^s$ their velocities. It is worthy to note here that the fluid continuum is assumed to obey to a laminar and sufficiently slow flow. Eq.(4) characterizes the inequality of local
entropy (i.e. Clausius - Duhem inequality) with $\mathcal{D}$ the total dissipation, $s$ standing for the overall entropy per unit volume of $d\Omega$ calculated by $s = \rho^s (1-\phi) s_s^s + \rho^f \phi s^f$ with $s_s^s$ and $s^f$ the specific entropy of the skeleton and that of the fluid part per unit mass<fnref n="3" />, respectively. Indeed, these two equations constitute an approximation, under small perturbations assumption, of the Lagrangian local form of the energy balance equation and that of the entropy inequality formulated by Coussy (2004) (see for instance his Eqs.(3.21) and Eq.(3.23) without the consideration of $r$ but with that of the gravity and inertia forces). Readers could be referred to Appendix A for more detail of their development.

For the sake of a better analysis in next sections, we indicate at this stage that the skeleton constituent of the overall internal energy, i.e. $e^s$, is a function of its infinitesimal strain $\varepsilon$, its entropy $s^s$ and other state variables (e.g. internal variables), whereas the fluid part, i.e. $e^f$, depends only on its specific entropy $s^f$ due to the incompressibility of fluid.

## 2.2. Coupling poroelasticity with gradient damage

This section deals with the gradient damage model for saturated porous media by incorporating the internal variables, which are the damage of the skeleton, denoted by $\alpha$, and its gradient $\nabla\alpha$ which characterizes its non-local effect. Moreover, in order to carry out the corresponding thermodynamics-based formulation, we mainly rely on an approach introduced by Frémond and Nedjar (1996) and recently followed for gradient plasticity theory by Nguyen (2016) both under the assumption of infinitesimal transformations. In this context, the power of interior forces requires not only the classical term $\sigma : \dot{\varepsilon}$, the pore fluid pressure $p$ together with the filtration velocity $v$ but also the microscopic forces associated to damage $\alpha$ and its gradient $\nabla\alpha$, respectively denoted by $Y_\alpha$ and $Y_{\nabla\alpha}$. Relying then on Coussy (1995) and Frémond and Nedjar (1996) (see also Frémond (2002)), one has:

$$ \mathcal{P}_{int} = - \int_{\Omega} [\sigma : \dot{\varepsilon} - \operatorname{div}(p\mathbf{v})] \, d\Omega - \int_{\Omega} (Y_{\alpha}\dot{\alpha} + \mathbf{Y}_{\nabla\alpha} \cdot \nabla\dot{\alpha}) \, d\Omega \quad (5) $$

By considering the above equation and starting from Eqs.(1) and (2), similar calculations to those which lead to (3) and (4) (see Appendix A) gives the local energy balance equation as well as the Clausius-Duhem inequality for the saturated porous media involving gradient damage, which respectively read:

$$ \dot{e} = \sigma : \dot{\varepsilon} + r - \operatorname{div}(q) - \operatorname{div} \left[ \rho^f \mathbf{v} \left( e^f + \frac{p}{\rho^f} \right) \right] + Y_\alpha \dot{\alpha} + \mathbf{Y}_{\nabla\alpha} \cdot \nabla \dot{\alpha} \quad (6) $$

$$ \mathcal{D} = \sigma : \dot{\varepsilon} - \dot{e} + T\dot{s} + T\operatorname{div}(\rho^f s^f \mathbf{v}) - \operatorname{div} \left[ \rho^f \mathbf{v} \left( e^f + \frac{p}{\rho^f} \right) \right] - \frac{\nabla T}{T} \cdot \mathbf{q} + Y_\alpha \dot{\alpha} + \mathbf{Y}_{\nabla\alpha} \cdot \nabla \dot{\alpha} \geq 0 \quad (7) $$

This is an extension of available result in Coussy (2004) by incorporation of damage terms. However, as the damage process is associated to the skeleton part, it would be more convenient to make use of its own energy potentials. For this purpose, we follow a similar procedure that was carried out in Coussy (2004) by introducing the local Lagrangian Helmholtz free energy per unit volume of $d\Omega$, denoted by $w$, satisfying the following partial Legendre transform:

$$ -w = Ts - e \quad (8) $$

<footnote n="3">³To distinguish from the local entropy of the skeleton per unit volume of $d\Omega$, which will be introduced in the following and denoted as $s^s$, we use the notation $s_s^s$ here.</footnote>
It can be readily understood that each of $w$ and $s$ are composed of a contribution from the skeleton and another one associated to the fluid part, namely:

$$w = w^s + \rho^f \phi w^f; \quad s = s^s + \rho^f \phi s^f \qquad (9)$$

where $w^s$ and $s^s$ denote the Helmholtz free energy and entropy of the skeleton per unit volume of $d\Omega$, respectively, while $w^f$ is the specific (per unit mass) Helmholtz free energy of the fluid part. By introducing Eqs.(8) and (9) into (7) and after some calculation, one gets:

$$\mathcal{D} = \underbrace{\boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} - \dot{w}^s - s^s \dot{T} + p \dot{\phi} + Y_{\alpha} \dot{\alpha} + \mathbf{Y}_{\nabla \alpha} \cdot \nabla \dot{\alpha}}_{\mathcal{D}_{\text{int}}} - \underbrace{-q \cdot \frac{\nabla T}{T}}_{\mathcal{D}_{\text{th}}} - \underbrace{-v \cdot \nabla p}_{\mathcal{D}_{\text{f}}} \geq 0 \qquad (10)$$

where $\mathcal{D}_{\text{int}}$ is the intrinsic dissipation of the skeleton, $\mathcal{D}_{\text{th}}$ denotes the dissipation of the thermal conduction and $\mathcal{D}_{\text{f}}$ that of the isothermal flow of the fluid.

Moreover, in order to keep the present work focused on the coupling between poroelasticity and gradient damage, we assume henceforward an isothermal evolution for all of the saturated porous system, i.e. $\dot{T} = 0$ and $\nabla T = 0$. Consequently, considering the above definition of $\mathcal{D}_{\text{int}}$, the rate of Helmholtz free energy of the skeleton reads:

$$\dot{w}^s = [(\boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} + Y_{\alpha} \dot{\alpha} + \mathbf{Y}_{\nabla \alpha} \cdot \nabla \dot{\alpha}) - \mathcal{D}_{\text{int}}] + p \dot{\phi} \qquad (11)$$

This equation reveals that the variation of $w^s$ can be divided into two parts: the first one in the bracket corresponds to non-dissipative mechanical work in the skeleton with fixed porosity; the other, i.e. $p\dot{\phi}$, accounts for variations of $w^s$ when all other parameters hold constant and is obviously associated with the influence of fluid pore pressure$^4$. Hence, one can define $w^s$ as a function of $\boldsymbol{\varepsilon}$, $\alpha$, $\nabla \alpha$ and $\phi$, which are its four natural variables, namely:

$$w^s := w^s(\boldsymbol{\varepsilon}, \phi, \alpha, \nabla \alpha) \qquad (12)$$

### 2.3. Thermodynamic forces and state laws

Since the state variables are independent, Eq.(11) can be recast into:

$$\dot{w}^s = \frac{\partial w^s}{\partial \boldsymbol{\varepsilon}} : \dot{\boldsymbol{\varepsilon}} + \frac{\partial w^s}{\partial \alpha} \dot{\alpha} + \frac{\partial w^s}{\partial \nabla \alpha} \cdot \nabla \dot{\alpha} + \frac{\partial w^s}{\partial \phi} \dot{\phi} = [(\boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} + Y_{\alpha} \dot{\alpha} + \mathbf{Y}_{\nabla \alpha} \cdot \nabla \dot{\alpha}) - \mathcal{D}_{\text{int}}] + p \dot{\phi} \qquad (13)$$

which allows to identify the following non-dissipative thermodynamics forces (denoted with superscript “nd”) as:

$$p^{nd} = \frac{\partial w^s}{\partial \phi} \quad \text{and} \quad \boldsymbol{\sigma}^{nd} = \frac{\partial w^s}{\partial \boldsymbol{\varepsilon}}; \quad Y_{\alpha}^{nd} = \frac{\partial w^s}{\partial \alpha}; \quad \mathbf{Y}_{\nabla \alpha}^{nd} = \frac{\partial w^s}{\partial \nabla \alpha}; \qquad (14)$$

Moreover, as the strain field $\boldsymbol{\varepsilon}$ and the Lagrangian porosity $\phi$ are both non-dissipative, one has $\boldsymbol{\sigma}^d = \boldsymbol{\sigma} - \boldsymbol{\sigma}^{nd} = 0$ and $p^d = p - p^{nd} = 0$, where superscript “d” represent the dissipative thermodynamics force of related variable. Eqs.(14) can then be recast into:

$$p = p^{nd} = \frac{\partial w^s}{\partial \phi}; \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{nd} = \frac{\partial w^s}{\partial \boldsymbol{\varepsilon}}; \quad Y_{\alpha}^{nd} = \frac{\partial w^s}{\partial \alpha}; \quad \mathbf{Y}_{\nabla \alpha}^{nd} = \frac{\partial w^s}{\partial \nabla \alpha}; \qquad (15)$$

<footnote>$^4$This reasoning is inspired by Maitournam (2013), who considered the case without poroelasticity and gradient damage.</footnote>
which are the state laws.

Yet, since the porosity $\phi$ is less “controllable”, it would be more appropriate to adopt its conjugate variable, i.e. the fluid pressure $p$, as an external state variable. In this context, it can be preferable to consider, instead of the density of the Helmholtz free energy $w^s$, the Lagrangian thermodynamics potential of the skeleton per unit volume, denoted by $\psi^s$, which can be obtained by the following partial Legendre transform (with respect to $\phi$) of $w^s$:

$$-\psi^s = p\phi - w^s \qquad (16)$$

It follows that:

$$\psi^s = \psi^s(\varepsilon, p, \alpha, \nabla\alpha) \qquad (17)$$

Note that as same as the Lagrangian porosity $\phi$, the pore fluid pressure $p$ is also a non-dissipative state variable. Particular attention must also be payed to the fact that $w^s$ is convex with respect to all of its arguments. It follows by duality that $\psi^s$ is a concave function of the pressure field $p$, while its convexity with respect to other variables still holds.

Considering Eqs.(10) and (16), the intrinsic dissipation function can be rewritten as:

$$D_{\text{int}} = \sigma : \dot{\varepsilon} - \dot{\psi}^s - \dot{p}\phi + Y_{\alpha}\dot{\alpha} + Y_{\nabla\alpha}\cdot\nabla\dot{\alpha} \geq 0 \qquad (18)$$

From similar calculations that lead to Eq.(15)<fnref n="5" />, the non-dissipative forces associated with the state variables $\varepsilon, p, \alpha$ and $\nabla\alpha$ can obtained as:

$$\sigma = \frac{\partial \psi^s}{\partial \varepsilon}; \quad \phi = -\frac{\partial \psi^s}{\partial p}; \quad Y_{\alpha}^{nd} = \frac{\partial \psi^s}{\partial \alpha}; \quad Y_{\nabla\alpha}^{nd} = \frac{\partial \psi^s}{\partial \nabla\alpha}; \qquad (19)$$

Inserting Eqs.(19) into (18) yields:

$$D_{\text{int}} = Y_{\alpha}^{d}\dot{\alpha} + Y_{\nabla\alpha}^{d}\cdot\nabla\dot{\alpha} \geq 0 \qquad (20)$$

where $Y_{\alpha}^d$ and $Y_{\nabla\alpha}^d$ are respectively the dissipative (i.e. irreversible) thermodynamic forces of $\alpha$ and $\nabla\alpha$, satisfying:

$$Y_{\alpha}^{d} = Y_{\alpha} - Y_{\alpha}^{nd}; \quad Y_{\nabla\alpha}^{d} = Y_{\nabla\alpha} - Y_{\nabla\alpha}^{nd} \qquad (21)$$

**Remark 1.** It should be noted at this stage that the skeleton should be thermodynamically equilibrated at any given moment. Due to this reason, the porosity calculated in the second equation of (19) must be recognized as that at these "equilibrium states". Moreover, by referring to the calculation from Eq.(7) to (18), it can be readily understood that the fluid transport is coupled with the deformation of the skeleton. In other words, the latter should tend to make room for entering the fluid mass content (i.e. $m = \rho^f\phi$) into the infinitesimal skeleton volume, which corresponds to a "kinetic behavior" of the porosity at the same moment.

In addition, as the pore fluid pressure $p$ and the porosity $\phi$ are both non-dissipative, the reversible process of the skeleton is described by the power $\psi^s + \dot{p}\phi$, whose stationarity allows to retrieve the state equation of $p$ (i.e. second equation in (19)). Readers are referred to Appendix B for more details. Note also that a similar situation was mentioned in Heuzé and Stainier (2022) (see also Yang et al. (2006); Stainier (2013)) in the context of the thermomechanics coupled with heat conduction.

<footnote n="5">⁵A more straightforward calculation of Eq.(19) consist in rewriting (15) by taking into account the partial Legendre transform (16).</footnote>
## 2.4. Evolution laws for the coupled poroelasticity-gradient damage model - the dissipation potentials

We are looking now for a dissipation potential, denoted by $\varphi$, which will allow to describe the irreversible evolution of quasi-brittle saturated porous media. $\varphi$ of the rate of the internal state variables of the skeleton (i.e., $\dot{\alpha}$ and $\nabla\dot{\alpha}$) and parametrically depending on its current state variables. Moreover, owing to the fluid viscosity and the permeability of the porous media, the fluid presents some spatial loss of pressure during its transport. Thus, the gradient of fluid pressure $\nabla p$ is considered as an additional variable of $\varphi$. Again, by considering the non-dissipative character of the strain field $\varepsilon$ and the pressure $p$, $\varphi$ is defined as:

$$ \varphi := \varphi(\dot{\alpha}, \nabla\dot{\alpha}, \nabla p; \alpha, \nabla\alpha) \qquad (22) $$

which can be treated as a jointed potential, such that:

$$ \varphi = \varphi^s + \varphi^f \qquad (23) $$

where the skeleton part is defined by

$$ \varphi^s := \varphi^s(\dot{\alpha}, \nabla\dot{\alpha}; \alpha, \nabla\alpha) \qquad (24) $$

Concerning the fluid part $\varphi^f$, and in link to the definition of the fluid dissipation $\mathcal{D}_f$ (in Eq.(10)), the dissipation potential depends on $\nabla p$, which is defined as:

$$ \varphi^f := \varphi^f(\nabla p); \quad \text{such that} \quad \mathbf{v} = \frac{\partial \varphi^f}{\partial \nabla p} \qquad (25) $$

Note that $\varphi^s$ must be convex, while $\varphi^f$ is concave with respect to their own arguments. Moreover, at this stage, it is emphasized that an assumption of the classical linear Darcy flow will be made in Section 4 for $\varphi^f$.

**Remark 2.** Alternatively, the fluid dissipation could be also interpreted as induced by the momentum of fluid particles. In this point of view, it is expressed by a kinetic energy function of the filtration velocity $\mathbf{v}$ that is obtained through the Legendre-Fenchel transform of $\varphi^f$ with respect to $\nabla p$, namely:

$$ -\varphi^{f*}(\mathbf{v}) = \sup\{\nabla p \cdot \mathbf{v} - \varphi^f(\nabla p)\} \quad \text{such that} \quad \nabla p = -\frac{\partial}{\partial \mathbf{v}}\varphi^{f*}(\mathbf{v}) \qquad (26) $$

where $\varphi^{f*}$ is convex with respect to $\mathbf{v}$.

In the context of Generalized Standard Materials (GSM) (Halphen and Nguyen, 1975), $\varphi^s$ and $\varphi^f$ are indeed the potentials describing the dissipative processes, respectively corresponding to the gradient damage and the fluid transport, being convex with respect to $\dot{\alpha}$, $\nabla\dot{\alpha}$ but concave to $\nabla p$. Having in hand Eq.(22) - (23), one can get the non-zero dissipative thermodynamics forces$^6$:

$$ Y_{\alpha}^{d} = \frac{\partial \varphi}{\partial \dot{\alpha}}; \quad \mathbf{Y}_{\nabla \alpha}^{d} = \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \qquad (27) $$

<footnote>$^6$In the case where the potential $\varphi$ is not differentiable, one should refer to the notion of subdifferential $\partial\varphi$, e.g. the evolution law for $\alpha$ then taking the form: $A_{\alpha} \in \partial_{\dot{\alpha}}\varphi$. The irreversible force $A_{\alpha}$ is said to belong to the sub-gradient of $\varphi$ at the considered point.</footnote>
with respect to $\alpha$ and $\nabla\alpha$, respectively.

It readily follows that:

$$Y_{\alpha} = \frac{\partial \psi^s}{\partial \alpha} + \frac{\partial \varphi}{\partial \dot{\alpha}}; \quad Y_{\nabla\alpha} = \frac{\partial \psi^s}{\partial \nabla\alpha} + \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \quad (28)$$

Let us summarize here the main elements of the thermodynamics-based gradient damage model of saturated porous media that has been established above:

* The starting point is an appropriate choice of the state variable of the skeleton i.e. the infinitesimal strain $\varepsilon$, the fluid pressure $p$, the damage $\alpha$ and its gradient $\nabla\alpha$.

* With the above choice, the reversible evolution of the saturated porous media with evolving damage is given by the thermodynamics potential of the skeleton $\psi^s$ (see Eq.(17)), while the irreversible process is described by means of the dissipation potential $\varphi = \varphi^s + \varphi^f$, which is assumed to be the superposition of the skeleton dissipation $\varphi^s$ and that of fluid $\varphi^f$ (see Section 2.4).

These two potentiels, $\psi^s$ and $\varphi$, allow to obtain the whole set of equations defining the searched constitutive laws. They will be very useful for the variational characterization of the model (sections 3 and 4) in view of structural analysis that will be presented in sections 5 and 6.

### 3. General formulation of incremental variational principle

In this section, we aim at proposing an incremental variational framework for the gradient damage model of saturated porous media. To this end, based on the thermodynamics framework developed in Section 2, we first propose a rate form of total energy functional, whose integral over the time increment allows to construct an incremental energy functional. The stationarity of the latter one will be demonstrated as being able to describe the governing equations of the fully coupled problem.

#### 3.1. Incremental variational principle

We consider that the whole system is subjected to some mechanical and/or external actions due to fluid. Specifically, the mechanical loads could be composed of surface forces $T$ applied on the Neumann boundary $\partial_T\Omega$ and an external displacement $\bar{u}$ on the Dirichlet part $\partial_u\Omega$, where $\partial_T\Omega \cup \partial_u\Omega = \partial\Omega$ and $\partial_T\Omega \cap \partial_u\Omega = \emptyset$. Concerning the external fluid charge, the boundary of the porous media could be subjected to a fluid injection $Q$ on its Neumann surface $\partial_Q\Omega$ and a fluid pressure $\bar{p}$ on the Dirichlet border $\partial_p\Omega$, for which we have $\partial_Q\Omega \cup \partial_p\Omega = \partial\Omega$ and $\partial_Q\Omega \cap \partial_p\Omega = \emptyset$. Moreover, the porous media might also be subjected to a body force per unit volume, denoted by $F$. We hence introduce the following external power functional:

$$P_{\text{ext}} := \int_{\Omega} \boldsymbol{F} \cdot \dot{\boldsymbol{u}} \, \mathrm{d}\Omega + \int_{\partial_T \Omega} \boldsymbol{T} \cdot \dot{\boldsymbol{u}} \, \mathrm{d}S + \int_{\partial_Q \Omega} Q p \, \mathrm{d}S \quad (29)$$

Inspired from some available studies (see for instance Ortiz and Stainier (1999) for elasto(visco)-plasticity, Yang et al. (2006); Stainier (2013) for thermo-mechanics, Nguyen (2016) for gradient plasticity, etc.), let us introduce the rate form of total energy functional for the poroelastic material with gradient damage:

$$\Pi(\dot{\boldsymbol{u}}, \dot{p}, \phi, \dot{\alpha}; p) = \int_{\Omega} \psi^s(\dot{\boldsymbol{\epsilon}}, \dot{p}, \dot{\alpha}, \nabla \dot{\alpha}) + \dot{p}\phi + \varphi(\dot{\alpha}, \nabla \dot{\alpha}, \nabla p; \alpha, \nabla \alpha) \, \mathrm{d}\Omega - P_{\text{ext}} \quad (30)$$
Recall that in the above equation, the thermodynamics potential $\dot{\psi}^s$ along with $\dot{p}\phi$ characterize the reversible power per unit volume of $d\Omega$. It is important to emphasize again that the porosity simultaneously presents some equilibrium and kinetic behaviors at any given time $t$ (see Remark 1). For the sake of a better comprehension, we denote here the so-called “kinetic porosity” by $\hat{\phi}$ at the pseudo non-equilibrium state $t^-$. The kinetic porosity $\hat{\phi}$ must be carefully distinguished from the so-called “equilibrium porosity” $\phi$ at the equilibrium state $t^+$, which is calculated via the corresponding state law (see Eq.(19)). The characterization of this kinetic porosity will be done in the following by considering Eqs.(36) and (37) in the variational point of view. In this perspective, we rewrite Eq.(30) as:

$$ \Pi(\dot{\mathbf{u}}, \dot{p}, \hat{\phi}, \dot{\alpha}; p) = \int_{\Omega} \dot{\psi}^s(\dot{\varepsilon}, \dot{p}, \dot{\alpha}, \nabla \dot{\alpha}) + \dot{p}\hat{\phi} + \varphi(\dot{\alpha}, \nabla \dot{\alpha}, \nabla p; \alpha, \nabla \alpha) \, d\Omega - P_{\text{ext}} \quad (31) $$

We now proceed to establish an incremental variational problem in order to reduce time-dependent problem such that the solution is expected to be given by the optimization of a total energy functional in the incremental form. Specifically, we discretize the total time path $[t_0, t_f]$ into a sequence of $t_0 = 0, ..., t_n, t_{n+1}, ..., t_N = t_f$, and seek to characterize the state at each of these time steps. This will be performed by consistently optimizing in any time interval $[t_n, t_{n+1}]$ the state $(\varepsilon_{n+1}, p_{n+1}, \hat{\phi}_{n+1}, \alpha_{n+1})$ at $t_{n+1}$ with that of $(\varepsilon_n, p_n, \hat{\phi}_n, \alpha_n)$ at $t_n$ being known. In this perspective, we define the following functional:

$$ I_n(\varepsilon_{n+1}, p_{n+1}, \hat{\phi}_{n+1}, \alpha_{n+1}) = \int_{t_n}^{t_{n+1}} \Pi(\dot{\mathbf{u}}, \dot{p}, \hat{\phi}, \dot{\alpha}; p) \, dt \quad (32) $$

By considering the integral of the exact differential $\dot{\psi}^s$, substitution of Eqs.(17), (22) and (29) into (32) yields:

$$ \begin{aligned} I_n(\mathbf{u}_{n+1}, p_{n+1}, \hat{\phi}_{n+1}, \alpha_{n+1}) = & \int_{\Omega} \psi^s_{n+1}(\mathbf{u}_{n+1}, p_{n+1}, \alpha_{n+1}, \nabla \alpha_{n+1}) - \psi^s_n(\mathbf{u}_n, p_n, \alpha_n, \nabla \alpha_n) \, d\Omega \\ & + \int_{t_n}^{t_{n+1}} \int_{\Omega} \dot{p}\hat{\phi} + \varphi(\dot{\alpha}, \nabla \dot{\alpha}, \nabla p; \alpha, \nabla \alpha) \, d\Omega \, dt - \int_{t_n}^{t_{n+1}} \int_{\partial Q \cap \Omega} Qp \, dS \, dt \\ & - \int_{\Omega} \mathbf{F}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \, d\Omega - \int_{\partial T \cap \Omega} \mathbf{T}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \, dS \, dt \end{aligned} \quad (33) $$

Taking variation of Eq.(33) with respect to the fluid pressure $p$ gives the related weak form:

$$ \int_{\Omega} \frac{\partial \psi^s_{n+1}}{\partial p_{n+1}} \delta p_{n+1} d\Omega + \int_{t_n}^{t_{n+1}} \int_{\Omega} \left( \hat{\phi} \delta \dot{p} + \frac{\partial \varphi}{\partial \nabla p} \cdot \nabla \delta p \right) d\Omega \, dt - \int_{t_n}^{t_{n+1}} \int_{\partial Q \cap \Omega} Q \delta p \, dS \, dt = 0 \quad (34) $$

whose integration by parts with respect to time and space yields:

$$ \begin{aligned} & \int_{\Omega} \left( \frac{\partial \psi^s_{n+1}}{\partial p_{n+1}} + \hat{\phi}_{n+1} \right) \delta p_{n+1} d\Omega - \int_{t_n}^{t_{n+1}} \int_{\Omega} \left[ \dot{\hat{\phi}} + \text{div} \left( \frac{\partial \varphi}{\partial \nabla p} \right) \right] \delta p \, d\Omega \, dt \\ & + \int_{t_n}^{t_{n+1}} \int_{\partial Q \cap \Omega} \left( \frac{\partial \varphi}{\partial \nabla p} \cdot \mathbf{n} - Q \right) \delta p \, dS \, dt = 0 \end{aligned} \quad (35) $$
Ensuring stationarity with respect to $p$ furnishes at $t = t_{n+1}$:

$$ \hat{\phi}_{n+1} = -\frac{\partial \psi^{s}_{n+1}}{\partial p_{n+1}}; \quad \forall \mathbf{x} \in \Omega \qquad (36) $$

and in the time interval $[t_n, t_{n+1}]$

$$ \dot{\hat{\phi}} + \operatorname{div}\left(\frac{\partial \varphi}{\partial \nabla p}\right) = 0; \quad \forall \mathbf{x} \in \Omega \qquad (37) $$

$$ \frac{\partial \varphi}{\partial \nabla p} \cdot \mathbf{n} = Q; \quad \forall \mathbf{x} \in \partial Q \Omega \qquad (38) $$

By considering (19) and (36), it can be readily shown that $\hat{\phi}_{n+1} = \phi_{n+1}$ at $t_{n+1}$. Also, from (25) and (23), together with (37), the rate of $\hat{\phi}$ satisfies the mass conservation of the fluid in $[t_n, t_{n+1}]$, i.e. $\dot{\hat{\phi}} + \operatorname{div}(\mathbf{v}) = 0$ with the corresponding Neumann boundary condition on $\partial Q \Omega$, such that $\mathbf{v} \cdot \mathbf{n} = Q$.

By analogy, taking variations of Eq. (33) with respect to the other variables yields the following weak forms respectively for $\mathbf{u}$, $\alpha$ and $\hat{\phi}$:

$$ -\int_{\Omega} \left[ \operatorname{div} \left( \frac{\partial \psi^{s}_{n+1}}{\partial \varepsilon_{n+1}} \right) + \mathbf{F}_{n+1} \right] \cdot \delta \mathbf{u}_{n+1} \, d\Omega + \int_{\partial T \Omega} \left( \frac{\partial \psi^{s}_{n+1}}{\partial \varepsilon_{n+1}} \cdot \mathbf{n} - \mathbf{T}_{n+1} \right) \cdot \delta \mathbf{u}_{n+1} \, dS \qquad (39) $$

$$ \begin{aligned} & \int_{\Omega} \left[ \left( \frac{\partial \psi^{s}_{n+1}}{\partial \alpha_{n+1}} + \frac{\partial \varphi}{\partial \dot{\alpha}} \right)_{|n+1} - \operatorname{div} \left( \frac{\partial \psi^{s}_{n+1}}{\partial \nabla \alpha_{n+1}} + \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \right)_{|n+1} \right] \delta \alpha_{n+1} \, d\Omega \\ & + \int_{\partial \Omega} \left[ \left( \frac{\partial \psi^{s}_{n+1}}{\partial \nabla \alpha_{n+1}} + \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \right)_{|n+1} \cdot \mathbf{n} \right] \delta \alpha_{n+1} \, dS \\ & + \int_{t_n}^{t_{n+1}} \int_{\Omega} \left\{ \left[ \frac{\partial \varphi}{\partial \alpha} - \frac{d}{dt} \left( \frac{\partial \varphi}{\partial \dot{\alpha}} \right) \right] - \operatorname{div} \left[ \frac{\partial \varphi}{\partial \nabla \alpha} - \frac{d}{dt} \left( \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \right) \right] \right\} \delta \alpha \, d\Omega \, dt \\ & + \int_{t_n}^{t_{n+1}} \int_{\partial \Omega} \left\{ \left[ \frac{\partial \varphi}{\partial \nabla \alpha} - \frac{d}{dt} \left( \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \right) \right] \cdot \mathbf{n} \right\} \delta \alpha \, dS \, dt = 0 \end{aligned} \qquad (40) $$

$$ \int_{t_n}^{t_{n+1}} \int_{\Omega} \dot{p} \, \delta \hat{\phi} \, d\Omega \, dt = 0 \qquad (41) $$

where $|_{n+1}$ symbolizes the estimate of considered quantity at time $t_{n+1}$.
Enforcing the stationarities of Eqs.(39) - (41) yields at $t_{n+1}$:

$$ egin{cases} 
\operatorname{div}\left(\frac{\partial \psi^{s}_{n+1}}{\partial \varepsilon_{n+1}}\right) + \mathbf{F}_{n+1} = \mathbf{0}; & \forall \mathbf{x} \in \Omega \\[2ex] 
\frac{\partial \psi^{s}_{n+1}}{\partial \varepsilon_{n+1}} \cdot \mathbf{n} = \mathbf{T}_{n+1}; & \forall \mathbf{x} \in \partial_{T} \Omega \\[2ex] 
\left(\frac{\partial \psi^{s}_{n+1}}{\partial \alpha_{n+1}} + \frac{\partial \varphi}{\partial \dot{\alpha}}|_{n+1}\right) - \operatorname{div}\left(\frac{\partial \psi^{s}_{n+1}}{\partial \nabla \alpha_{n+1}} + \frac{\partial \varphi}{\partial \nabla \dot{\alpha}}|_{n+1}\right) = 0; & \forall \mathbf{x} \in \Omega \\[2ex] 
\left(\frac{\partial \psi^{s}_{n+1}}{\partial \nabla \alpha_{n+1}} + \frac{\partial \varphi}{\partial \nabla \dot{\alpha}}|_{n+1}\right) \cdot \mathbf{n} = 0; & \forall \mathbf{x} \in \partial \Omega 
\end{cases} 
\quad (42) $$

and in the interval $[t_n, t_{n+1}]$:

$$ egin{cases} 
\left[ \frac{\partial \varphi}{\partial \alpha} - \frac{d}{dt} \left( \frac{\partial \varphi}{\partial \dot{\alpha}} \right) \right] - \operatorname{div} \left[ \frac{\partial \varphi}{\partial \nabla \alpha} - \frac{d}{dt} \left( \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \right) \right] = 0 & \forall \mathbf{x} \in \Omega \\[2ex] 
\left[ \frac{\partial \varphi}{\partial \nabla \alpha} - \frac{d}{dt} \left( \frac{\partial \varphi}{\partial \nabla \dot{\alpha}} \right) \right] \cdot \mathbf{n} = 0 & \forall \mathbf{x} \in \partial \Omega \\[2ex] 
\dot{p} = 0 & \forall \mathbf{x} \in \Omega 
\end{cases} 
\quad (43) $$

Eqs.(36) and (42) furnish the Euler-Lagrange equations of $I_n$ and the Neumann boundary conditions at $t_{n+1}$. Eq.(43) along with (37) and (38) allow to determine the optimization path<fnref n="7" /> of the incremental energy functional $I_n$. All of them are expressed thanks to the thermodynamics potentials allowing to take into account the non-local damage effects. Essentially, these weak forms make $I_n$ to be stationary in $[t_n, t_{n+1}]$. Let us suppose the thermodynamics potentials as classically being quadratic, and recall their convexity as discussed in Sections 2.2 and 2.4. It follows that $I_n$ must be convex with respect to $u_{n+1}$ and $\alpha_{n+1}$ but concave with respect to $p_{n+1}$. We can hence establish the following incremental variational principle:

$$ (\mathbf{u}_{n+1}, p_{n+1}, \hat{\phi}_{n+1}, \alpha_{n+1}) = \operatorname{Arg} \inf_{\mathbf{u}'_{n+1}, \alpha'_{n+1}, p'_{n+1}, \hat{\phi}'_{n+1}} \operatorname{stat} \left\{ I_n \left( \mathbf{u}'_{n+1}, p'_{n+1}, \hat{\phi}'_{n+1}, \alpha'_{n+1} \right) \right\} \quad (44) $$

with $(\mathbf{u}_{n+1}, p_{n+1}, \hat{\phi}_{n+1}, \alpha_{n+1})$ henceforth the estimate of the expected solution and $(\mathbf{u}'_{n+1}, p'_{n+1}, \hat{\phi}'_{n+1}, \alpha'_{n+1})$ denoting the admissible fields of the displacement, the fluid pressure, the kinetic porosity as well as that of the damage, respectively. We recall that $I_n$ has been expressed in Eq.(33).

### 3.2. A reduced form of the variational scheme

It can be understood that the application of the incremental variational principle (44) requires a calculation of the optimization path, whose numerical implementation might be performed via a fully implicit algorithm. However, as discussed in Ortiz and Martin (1989) (see also Yang et al. (2006)), this path calculation could be treated as avoidable except for some particular case and

<footnote n="7">⁷See for Ortiz and Repetto (1999) an account for minimization path.</footnote>
hence, a weak quantification of $I_n$ compatible with field equations at $t_{n+1}$ could be sufficient. In this light, we first adopt here the backward Euler approximations for the following rate fields:

$$ \dot{\mathbf{u}} \approx \frac{\mathbf{u}_{n+1} - \mathbf{u}_n}{\Delta t}; \quad \dot{\alpha} \approx \frac{\alpha_{n+1} - \alpha_n}{\Delta t}; \quad \dot{p} \approx \frac{p_{n+1} - p_n}{\Delta t}; \qquad (45) $$

where $\Delta t = t_{n+1} - t_n$.

Moreover, as $\hat{\phi}$ being kinetic, which is, as aforementioned, consistent with the fluid filtration at the pseudo non-equilibrium state, one of a family of its appropriate numerical identifications in $[t_n, t_{n+1}]$ could be the porosity at the beginning of the current increment, namely, $\hat{\phi} = \phi_n$. In this spirit and still from an incremental point of view, the continuity of $\hat{\phi}$ in time dimension makes:

$$ \hat{\phi} = \phi_n = - \frac{\partial \psi^s}{\partial p} (\mathbf{u}_n, p_n, \alpha_n, \nabla \alpha_n) \qquad (46) $$

which is indeed its optimization at the last increment. It should be emphasized here that this treatment is not unique and expressed at this stage by Eq.(46) just as a modeling choice purpose. Since the kinetic characters of $\hat{\phi}$ can be applied at any instant between $t_n$ and $t_{n+1}$ except for $t_{n+1}$, the latter is considered as the stationary state to be calculated. In short, each of these adequate choices is of course an approximation, but as will be shown later, is able to suitably describe the coupled constitutive relations and to deliver sufficiently reliable results when $\Delta t \mapsto 0$.

By substituting Eqs. (45) and (46) into (33) and making use of the rectangular rule for the time integration of $\phi$ as well as that of the term $Qp$, one gets the incremental energy to be minimized:

$$ \begin{aligned} I_n(\mathbf{u}_{n+1}, p_{n+1}, \alpha_{n+1}) \approx & \int_{\Omega} (\psi_{n+1}^s - \psi_n^s) + \phi_n (p_{n+1} - p_n) \, d\Omega + \Delta t \int_{\Omega} (\varphi_{n+1}^s + \varphi_{n+1}^f) \, d\Omega \\ & - \int_{\Omega} \mathbf{F}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \, d\Omega - \int_{\partial T_{\Omega}} \mathbf{T}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \, dS \\ & - \Delta t \int_{\partial Q_{\Omega}} Q_{n+1} p_{n+1} \, dS \end{aligned} \qquad (47) $$

Let us compute the variation of Eq.(47) with respect to $\mathbf{u}_{n+1}, p_{n+1}, \alpha_{n+1}$; the stationarity gives:

$$ \left\{ \begin{aligned} & \operatorname{div} \left( \frac{\partial \psi_{n+1}^s}{\partial \varepsilon_{n+1}} \right) + \mathbf{F}_{n+1} = 0 && \forall \mathbf{x} \in \partial \Omega \\ & \frac{\partial \psi_{n+1}^s}{\partial \varepsilon_{n+1}} \cdot \mathbf{n} - \mathbf{T}_{n+1} = 0; && \forall \mathbf{x} \in \partial T_{\Omega} \\ & \frac{-\frac{\partial \psi_{n+1}^s}{\partial p_{n+1}} - \phi_n}{\Delta t} + \operatorname{div} \left( \frac{\partial \varphi_{n+1}^f}{\partial \nabla p_{n+1}} \right) = 0; && \forall \mathbf{x} \in \Omega \\ & \frac{\partial \varphi_{n+1}^f}{\partial \nabla p_{n+1}} \cdot \mathbf{n} - Q_{n+1} = 0; && \forall \mathbf{x} \in \partial Q_{\Omega} \\ & \left( \frac{\partial \psi_{n+1}^s}{\partial \alpha_{n+1}} + \frac{\partial \varphi_{n+1}^s}{\partial \dot{\alpha}_{n+1}} \right) - \operatorname{div} \left( \frac{\partial \psi_{n+1}^s}{\partial \nabla \alpha_{n+1}} + \frac{\partial \varphi_{n+1}^s}{\partial \nabla \dot{\alpha}_{n+1}} \right) = \mathcal{O}(\Delta t); && \forall \mathbf{x} \in \Omega \\ & \left( \frac{\partial \psi_{n+1}^s}{\partial \nabla \alpha_{n+1}} + \frac{\partial \varphi_{n+1}^s}{\partial \nabla \dot{\alpha}_{n+1}} \right) \cdot \mathbf{n} = \mathcal{O}(\Delta t); && \forall \mathbf{x} \in \partial \Omega \end{aligned} \right. \qquad (48) $$
Furthermore, by considering (19) and (28), Eq.(48) can be readily recast into:

$$
\left\{
\begin{array}{ll}
\operatorname{div}(\boldsymbol{\sigma}_{n+1}) + \boldsymbol{F}_{n+1} = 0 & \forall \boldsymbol{x} \in \Omega \\
\boldsymbol{\sigma}_{n+1} \cdot \boldsymbol{n} - \boldsymbol{T}_{n+1} = 0; & \forall \boldsymbol{x} \in \partial_T \Omega \\
\frac{\phi_{n+1} - \phi_n}{\Delta t} + \operatorname{div}(\boldsymbol{v}_{n+1}) = 0; & \forall \boldsymbol{x} \in \Omega \\
\boldsymbol{v}_{n+1} \cdot \boldsymbol{n} - Q_{n+1} = 0; & \forall \boldsymbol{x} \in \partial_Q \Omega \\
Y_{\alpha}|_{n+1} - \operatorname{div}(\boldsymbol{Y}_{\nabla\alpha}|_{n+1}) = \mathcal{O}(\Delta t); & \forall \boldsymbol{x} \in \Omega \\
\boldsymbol{Y}_{\nabla\alpha}|_{n+1} \cdot \boldsymbol{n} = \mathcal{O}(\Delta t); & \forall \boldsymbol{x} \in \partial\Omega
\end{array}
\right.
\qquad (49)
$$

which are the governing equations and Neumann boundary conditions associated to $u_{n+1}$, $p_{n+1}$, $\alpha_{n+1}$, respectively.

As the functional $I_n(u_{n+1}, p_{n+1}, \alpha_{n+1})$ defined by Eq.(47) must also be separately convex with respect to $u_{n+1}$ and $\alpha_{n+1}$ but concave with $p_{n+1}$, the incremental variational principle of Eq.(44) could be reduced into the following form as dependent on three fields, such that:

$$ (u_{n+1}, p_{n+1}, \alpha_{n+1}) = \operatorname{Arg} \inf_{u_{n+1}', \alpha_{n+1}'} \sup_{p_{n+1}'} \left\{ I_n \left( u'_{n+1}, p'_{n+1}, \alpha'_{n+1}; \hat{\phi} \right) \right\} \quad (50) $$

with $\hat{\phi}$ being beforehand identified in the current increment.

It must be emphasized here that the incremental variational principle described by (50) together with (47) constitutes one of the most important contributions of the present study.

## 4. A variational model for poroelasticity coupled with gradient damage

The objective of this section is to apply the established incremental variational principle to a specific gradient damage model of saturated porous media. As an example, we will make use of some appropriate thermodynamic potentials by assuming that the skeleton of the porous media is poroelastic with an evolving damage process. The fluid transport in the pore space is supposed to be sufficiently slow and laminar obeying the linear Darcy's law.

### 4.1. Thermodynamics and dissipation potentials

We propose as:

* Thermodynamic potential for the Biot-type poroelastic damaged skeleton:

$$ \psi^s(\boldsymbol{\varepsilon}, p, \alpha) = \frac{1}{2} \boldsymbol{\varepsilon} : \mathbb{C}(\alpha) : \boldsymbol{\varepsilon} - p\phi_0 - (p - p_0)\boldsymbol{b}(\alpha) : \boldsymbol{\varepsilon} - \frac{(p - p_0)^2}{2N(\alpha)} \quad (51) $$

where $\mathbb{C}(\alpha)$, $\boldsymbol{b}(\alpha)$, $1/N(\alpha)$ are respectively the stiffness tensor, the Biot tensor and the Biot modulus of the damaged skeleton. Moreover, as usual for phase field models, we assume that $\mathbb{C}(\alpha) = g(\alpha)\mathbb{C}^s$ with $g(\alpha)$ the corresponding degradation function that is taken as $g(\alpha) = (1 - \alpha)^2$. $\mathbb{C}^s$ denotes the elasticity tensor for the sound skeleton depending on the
initial porosity $\phi_0$. $\mathbb{C}^s = (1 - \phi_0)\mathbb{C}^0$, the tensor $\mathbb{C}^0$ being the stiffness tensor of the solid phase. The following approximation is then adopted:

$$ \mathbb{C}(\alpha) = (1 - \alpha)^2(1 - \phi_0)\mathbb{C}^0 \quad (52) $$

Having Eq.(52) in hand, the Biot tensor $\mathbf{b}(\alpha) = b(\alpha)\mathbf{1}$ and the Biot modulus $1/N(\alpha)$ can be readily computed in the isotropic case as:

$$ \mathbf{b}(\alpha) = b(\alpha)\mathbf{1} \quad \text{and} \quad \frac{1}{N(\alpha)} = \frac{b(\alpha) - \phi_0}{\kappa^0}, \quad \text{with} \quad b(\alpha) = 1 - \frac{\kappa(\alpha)}{\kappa^0} \quad (53) $$

where $b(\alpha)$ is the Biot coefficient; $\kappa(\alpha)$ and $\kappa^0$ denote the compressibility modulus of the skeleton and that of its solid matrix, respectively. It can be readily seen that when $\alpha = 0$, $b$ is equal to the initial porosity $\phi_0$. This is in fact the lower bound of its estimate (Dormieux et al., 2006), but, seems not so adequate in practical and engineering point of view. In order to figure out this limitation, we calibrate it in an ad-hoc manner by assuming that

$$ b(\alpha) = 1 - (1 - \alpha)^2(1 - b_0) \quad (54) $$

with $b_0$ being the initial Biot coefficient that could be identified from experimental measurements or micromechanical formulation (see for instance Dormieux et al. (2006)). The above equation indicates a monotonic increasing variation of $b$ with respect to $\alpha$ from the initial value $b_0$ until the totally damage one, namely $b = 1$. Without resorting to micromechanical arguments, this can be primarily interpreted as the consequence of the skeleton stiffness degradation which is known to affect the poroelastic properties.

Next, the reversible thermodynamic forces obtained from the state laws reads:

$$ \begin{gathered} \sigma = \mathbb{C}(\alpha) : \varepsilon - (p - p_0)\mathbf{b}(\alpha); \quad \phi - \phi_0 = \mathbf{b}(\alpha) : \varepsilon + \frac{p - p_0}{N(\alpha)}; \\ Y_{\alpha}^{nd} = \frac{\partial \psi^s}{\partial \alpha}; \quad \mathbf{Y}_{\nabla\alpha}^{nd} = \mathbf{0} \end{gathered} \quad (55) $$

* Dissipation potential:

$$ \varphi(\dot{\alpha}, \nabla\dot{\alpha}, \nabla p; \alpha, \nabla\alpha) = \varphi^s(\dot{\alpha}, \nabla\dot{\alpha}; \alpha, \nabla\alpha) + \varphi^f(\nabla p) \quad (56) $$

For the purpose of a guiding example, the dissipation potential for the skeleton damage is supposed to be in the form (see also Kpotufe (2021)):

$$ \varphi^s(\dot{\alpha}, \nabla\dot{\alpha}, \nabla p; \alpha, \nabla\alpha) = Y_c(\alpha)\dot{\alpha} + 2l_0^2 w_1 \nabla\alpha \cdot \nabla\dot{\alpha} \quad (57) $$

where $Y_c(\alpha)$ denotes the critical damage energy at current state, $l_0$ is the material internal length and $w_1$ the first threshold in energy release rate that can be identified from different damage descriptions (e.g. Pham et al. (2011) for AT1 model and Bourdin et al. (2000) for AT2 model.). It can be readily observed that Eq.(57) is positively homogeneous of degree one with decoupling of the contribution of the internal variables $\alpha$ and that of $\nabla\alpha$, which leads the related optimization path to be vanished (see the first two equations of (43) which are automatically satisfied). In the present case, this corresponds to the simple dissipative
systems introduced by Ehrlacher and Fedelich (1989), for which the total dissipated energy until the current time is function only of the current value of internal variables. For this reason, it is convenient to make use of the current dissipated energy of local damage:

$$ \Phi^s = \int_0^t \varphi^s(\dot{\alpha}) \, dt \qquad (58) $$

Moreover, by introducing $\omega(\alpha) = \frac{1}{w_1} \int_0^\alpha Y_c(\beta)d\beta$ that describes the local contribution of the energy dissipated during the evolution of the damage until the current time, Eq.(58) can be recast into:

$$ \Phi^s = w_1 l_0 \left( \frac{\omega(\alpha)}{l_0} + l_0 \| \nabla \alpha \|^2 \right) \qquad (59) $$

where $\| \cdot \|$ symbolizes the $L^2$ norm.

The corresponding incremental dissipation over $[t_n, t_{n+1}]$ can then be readily calculated as:

$$ \Phi_n^s = \Phi^s(\alpha_{n+1}) - \Phi^s(\alpha_n) = w_1 l_0 \left[ \frac{\omega(\alpha_{n+1}) - \omega(\alpha_n)}{l_0} + l_0 (\| \nabla \alpha_{n+1} \|^2 - \| \nabla \alpha_n \|^2) \right] \qquad (60) $$

By looking at Eq.(47), this incremental formulation yields:

$$ \Delta t \int_{\Omega} \varphi_{n+1}^s d\Omega \approx \int_{\Omega} \Phi_n^s d\Omega; \quad \text{when } \Delta t \mapsto 0 \qquad (61) $$

It is worthy to note here that Eq.(61) allows to retrieve the closed-form of the related incremental dissipation for such a simple dissipative process without the approximation that was made in Eq.(47). Nevertheless, we indicate that the incremental variational approach developed in Section 3 can also be applied to a more general situation, such as those corresponding to dissipation potentials in which $\dot{\alpha}$ and $\nabla\alpha$ are coupled.

Coming now to the description of fluid diffusion, as aforementioned, we make use of the linear Darcy law, namely,

$$ v = -\frac{\mathbf{k}}{\mu^f} \cdot \nabla p \qquad (62) $$

with $\mu^f$ the dynamic viscosity of the fluid and $\mathbf{k}$ the permeability tensor. The corresponding dissipation potential $\varphi^f$ is:

$$ \varphi^f(\nabla p) = -\frac{1}{2} \nabla p \cdot \frac{\mathbf{k}}{\mu^f} \cdot \nabla p \qquad (63) $$

It is worthy to emphasize that in the case of moderate (diffuse) damage, there may be typically an isotropic increase of permeability. On the other hand, once the structure presents a highly localized damage band, assimilated to regularized cracks, one may have induced anisotropic change of permeability. Interestingly, these two aspects have been explored in the literature in the context of variational damage-based phase field modeling:

1. For instance, from a phenomenological point of view, Zhu et al. (2013) (see also Pillai et al. (2018); Mollaali et al. (2019)) proposed an exponential function for the definition of an isotropic permeability evolution, which reads $\mathbf{k} = k^s (\phi/\phi_0)^3 \exp(c_\alpha \alpha) \mathbf{1}$, with $k^s$ the initial isotropic permeability of the porous material and $c_\alpha$ the so-called damage-permeability effect coefficient.
2. Alternatively, still adopting a phenomenological point of view, Miehe et al. (2015) proposed a so-called modification of Darcy's law by adding an additional anisotropic permeability to the isotropic one of the sound skeleton. The former corresponds to the description of a Poiseuille-type flow in opened regularized crack, which is represented by a highly localized damage band (with $\alpha \mapsto 1$). This type of characterization for the permeability change was also recently adopted by some other works on the phase field modeling applied to poromechanics (see for instance Miehe and Mauthe (2016); Mauthe and Miehe (2017); Xia et al. (2017); Noii et al. (2021); Ulloa et al. (2022)).

In the present study, we follow and make use of the second approach with the anisotropic permeability tensor in the form (see for instance Mauthe and Miehe (2017)):

$$k(\varepsilon, \alpha, \nabla\alpha) = k^s \mathbf{1} + \alpha^\epsilon k^\alpha(\varepsilon, \nabla\alpha) \quad \text{with} \quad k^\alpha(\varepsilon, \nabla\alpha) = \left(\frac{\eta^2}{12} - k^s\right) (\mathbf{1} - \tilde{\mathbf{n}} \otimes \tilde{\mathbf{n}}) \quad (64)$$

where $\alpha^\epsilon k^\alpha$ represents the additional permeability with $\epsilon$ an interpolation parameter that allows to incorporate the transition from Darcy-type flow to Poiseuille-type flow. The vector $\tilde{\mathbf{n}} = \nabla\alpha/|\nabla\alpha|$ gives the normal orientation to the localized damage zone representing the regularized crack surface. $\eta = (\tilde{\mathbf{n}} \cdot \varepsilon \cdot \tilde{\mathbf{n}})L_\perp$ is a crack opening variable with $L_\perp$ a numerical parameter approximated by the element scale.

As already mentioned by Ulloa et al. (2022), this characterization of permeability tensor is heuristic for the representation of the regularized crack opening $\eta$ and the normal direction $\tilde{\mathbf{n}}$. However, it is worth noting that it can allow a rather complete description of the change of permeability from the moderate damage state to that of full damage, such that:

- by choosing a relatively high value for $\epsilon$ (e.g. $\epsilon = 50$, which will be adopted in the results sections), the additional term becomes negligible for $\alpha < 0.9$. In this regime, the permeability remains quite isotropic and approximately equal to that of the intact skeleton, i.e. $k \simeq k^s \mathbf{1}$.

- once $\alpha > 0.9$, $k$ significantly increases in an anisotropy manner. This idea has been also expressed in Dormieux and Kondo (2005) (see also Dormieux and Kondo (2004)), who have identified in their micromechanical analysis (of the coupling between permeability and damage) a critical threshold of the micro-crack density parameter (damage), beyond which there is an important change of the permeability.

- in the limit case where $\alpha = 1$ (regularized crack), one has from Eq. (64) that:

$$k(\varepsilon; \alpha = 1) = k^s \tilde{\mathbf{n}} \otimes \tilde{\mathbf{n}} + \frac{\eta^2}{12} (\mathbf{1} - \tilde{\mathbf{n}} \otimes \tilde{\mathbf{n}}) \quad (65)$$

which is precisely the expression considered for fracture anisotropic permeability by Dormieux and Kondo (2004) (based on Guéguen and Palciauskas (1994)).

For completeness, it must be mentioned here that $k$ is considered as a modeling parameter in the present work, which is evaluated in the previous state, namely, $k \simeq k(\epsilon_n, \alpha_n, \nabla\alpha_n)$, and will be updated at the end of the current increment. This explicit approach has also advantage to preserve the requisite variational structure of the model.

Substitution of Eqs. (56), (57) and (63) into (27) gives the dissipative thermodynamic forces of the proposed model:

$$\boldsymbol{\sigma}^d = \mathbf{0}; \quad Y_p^d = 0; \quad Y_\alpha^d = Y_c(\alpha); \quad \mathbf{Y}_{\nabla\alpha}^d = 2l_0^2 w_1 \nabla\alpha \quad (66)$$
4.2. Final expression of the incremental functional $I_n$

We now look at the so-called kinetic porosity $\hat{\phi}$, which is indispensable, as already analyzed in Section 3, for the thermodynamics description of the concerned system. Recall that one of the relevant choices in the current increment $[t_n, t_{n+1}]$ could be $\hat{\phi} = \phi_n$ (see Section 3.2), which can be computed as:

$$ \hat{\phi} = -\frac{\partial \psi^s}{\partial p}(\boldsymbol{\varepsilon}_n, p_n, \alpha_n) = \phi_0 + \mathbf{b}(\alpha_n) : \boldsymbol{\varepsilon}_n + \frac{p_n - p_0}{N(\alpha_n)} $$

Moreover, in order to well incorporate the effect of the damage on the degradation of poroelasticity during the variational calculation, we suppose that the dependence of $\mathbf{b}(\alpha)$ and $1/N(\alpha)$ on the damage field $\alpha$ in the current increment would be satisfied through $\alpha_{n+1}$, such that:

$$ \hat{\phi} \approx \phi_0 + \mathbf{b}(\alpha_{n+1}) : \boldsymbol{\varepsilon}_n + \frac{p_n - p_0}{N(\alpha_{n+1})} \quad (67) $$

and also:

$$ \psi_n^s \approx \frac{1}{2} \boldsymbol{\varepsilon}_n : \mathcal{C}(\alpha_n) : \boldsymbol{\varepsilon}_n - p_n \phi_0 - (p_n - p_0) \mathbf{b}(\alpha_{n+1}) : \boldsymbol{\varepsilon} - \frac{(p_n - p_0)^2}{2N(\alpha_{n+1})} \quad (68) $$

These two equations lead to:

$$ \dot{p}\hat{\phi} \approx \frac{p_{n+1} - p_n}{\Delta t} \hat{\phi} = \frac{p_{n+1} - p_n}{\Delta t} \left[ \phi_0 + \mathbf{b}(\alpha_{n+1}) : \boldsymbol{\varepsilon}_n + \frac{p_n - p_0}{N(\alpha_{n+1})} \right] \quad (69) $$

and

$$ \begin{aligned} \psi_{n+1}^s - \psi_n^s \approx & \frac{1}{2} \boldsymbol{\varepsilon}_{n+1} : \mathcal{C}(\alpha_{n+1}) : \boldsymbol{\varepsilon}_{n+1} - \frac{1}{2} \boldsymbol{\varepsilon}_n : \mathcal{C}(\alpha_n) : \boldsymbol{\varepsilon}_n - (p_{n+1} - p_n) \phi_0 \\ & - \mathbf{b}(\alpha_{n+1}) : [(p_{n+1} - p_0) \boldsymbol{\varepsilon}_{n+1} - (p_n - p_0) \boldsymbol{\varepsilon}_n] - \frac{(p_{n+1} - p_n)(p_{n+1} + p_n - 2p_0)}{2N(\alpha_{n+1})} \end{aligned} \quad (70) $$

Gathering Eqs.(56), (57), (61), (63), (69) and (70), the incremental energy functional given by (47) can be recast in the form:

$$ \begin{aligned} I_n(\boldsymbol{\varepsilon}_{n+1}, \alpha_{n+1}, p_{n+1}) = & \int_{\Omega} \left[ \frac{1}{2} \boldsymbol{\varepsilon}_{n+1} : \mathcal{C}(\alpha_{n+1}) : \boldsymbol{\varepsilon}_{n+1} - \frac{1}{2} \boldsymbol{\varepsilon}_n : \mathcal{C}(\alpha_n) : \boldsymbol{\varepsilon}_n \right. \\ & \left. - (p_{n+1} - p_0) \mathbf{b}(\alpha_{n+1}) : (\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n) - \frac{(p_{n+1} - p_n)^2}{2N(\alpha_{n+1})} \right] d\Omega \\ & + \int_{\Omega} w_1 l_0 \left[ \frac{\omega(\alpha_{n+1}) - \omega(\alpha_n)}{l_0} + l_0 (\| \nabla \boldsymbol{\varepsilon}_{n+1} \|^2 - \| \nabla \boldsymbol{\varepsilon}_n \|^2) \right] d\Omega \quad (71) \\ & - \Delta t \int_{\Omega} \frac{1}{2} \nabla p_{n+1} \cdot \frac{\mathbf{k}}{\mu_f} \cdot \nabla p_{n+1} \, d\Omega - \int_{\Omega} \mathbf{F}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \, d\Omega \\ & - \int_{\partial_T \Omega} \mathbf{T}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \, dS - \Delta t \int_{\partial_Q \Omega} Q_{n+1} p_{n+1} \, dS \end{aligned} $$
In the context of variational optimization, for which the quantities that depend only on the state at $t_n$ can be considered constant, the above equation makes the same sense with respect to:

$$ 
\begin{aligned}
I_n(\varepsilon_{n+1}, \alpha_{n+1}, p_{n+1}) = & \int_{\Omega} \left[ \frac{1}{2} \varepsilon_{n+1} : \mathbb{C}(\alpha_{n+1}) : \varepsilon_{n+1} - (p_{n+1} - p_0) \mathbf{b}(\alpha_{n+1}) : (\varepsilon_{n+1} - \varepsilon_n) \right. \\
& \left. - \frac{(p_{n+1} - p_n)^2}{2N(\alpha_{n+1})} \right] d\Omega + \int_{\Omega} w_1 l_0 \left( \frac{\omega(\alpha_{n+1})}{l_0} + l_0 \| \nabla \alpha_{n+1} \|^2 \right) d\Omega \\
& - \Delta t \int_{\Omega} \frac{1}{2} \nabla p_{n+1} \cdot \frac{\mathbf{k}}{\mu_f} \cdot \nabla p_{n+1} d\Omega - \int_{\Omega} \mathbf{F}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) d\Omega \\
& - \int_{\partial T \Omega} \mathbf{T}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) dS - \Delta t \int_{\partial Q \Omega} Q_{n+1} p_{n+1} dS
\end{aligned} 
\quad (72) $$

Formulation (72) reveals that the incremental energy functional is quadratic of each of its arguments and separately convex with respect to $\mathbf{u}_{n+1}$ and $\alpha_{n+1}$, but concave with respect to $p_{n+1}$. Hence, the resolution of the concerned problem can be achieved by applying the incremental variational principle (50), which is based on the three fields dependent functional (72).

For completeness, we provide below the governing equations and corresponding boundary conditions, which are obtained based on the weak form of Eq. (72):

$$ 
\left\{
\begin{aligned}
& \operatorname{div} [\mathbb{C}(\alpha_{n+1}) : \boldsymbol{\varepsilon}_{n+1} - (p_{n+1} - p_0)\mathbf{b}(\alpha_{n+1})] + \mathbf{F}_{n+1} = \mathbf{0} \\
& [\mathbb{C}(\alpha_{n+1}) : \boldsymbol{\varepsilon}_{n+1} - (p_{n+1} - p_0)\mathbf{b}(\alpha_{n+1})] \cdot \mathbf{n} - \mathbf{T}_{n+1} = \mathbf{0} \\
& \frac{\partial \psi^s_{n+1}}{\partial \alpha_{n+1}} + w_1 \frac{\partial \omega(\alpha_{n+1})}{\partial \alpha_{n+1}} - 2w_1 l_0^2 \nabla^2 \alpha_{n+1} = 0 \\
& (2w_1 l_0^2 \nabla \boldsymbol{\alpha}_{n+1}) \cdot \mathbf{n} = 0 \\
& \frac{1}{\Delta t} \left[ \mathbf{b}(\alpha_{n+1}) : (\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n) + \frac{(p_{n+1} - p_n)}{N(\alpha_{n+1})} \right] - \operatorname{div} \left( \frac{\mathbf{k}}{\mu_f} \cdot \nabla p_{n+1} \right) = 0 \\
& \left( -\frac{\mathbf{k}}{\mu_f} \cdot \nabla p_{n+1} \right) \cdot \mathbf{n} + Q_{n+1} = 0
\end{aligned}
\right. 
\quad (73) $$

with $\nabla^2$ the Laplacian operator.

**Remark 3.** Following the methodology of Amor et al. (2009), the unilateral contact effect of damage under compressive loading is also implemented by calibrating the thermodynamic potential (i.e. Eq. (51)) as:

$$ 
\begin{aligned}
\psi^s(\boldsymbol{\varepsilon}, \alpha, p) = & (1 - \alpha)^2 (1 - \phi_0) \left[ \frac{\kappa^0}{2} \langle \operatorname{tr}(\boldsymbol{\varepsilon}) \rangle_+^2 + \mu \boldsymbol{\varepsilon}_{dev} : \boldsymbol{\varepsilon}_{dev} \right] + (1 - \phi_0) \frac{\kappa^0}{2} \langle \operatorname{tr}(\boldsymbol{\varepsilon}) \rangle_-^2 \\
& - p\phi_0 - (p - p_0) \mathbf{b}(\alpha) : \boldsymbol{\varepsilon} - \frac{(p - p_0)^2}{2N(\alpha)}
\end{aligned} 
\quad (74) $$
where $\kappa^0$ is the bulk modulus of the solid matrix and $\mu^0$ its shear modulus; $\varepsilon_{dev}$ is the deviatoric strain tensor; $\langle a \rangle_+ := (a + |a|)/2$, $\langle a \rangle_- := (a - |a|)/2$ are the positive and negative ramp functions, respectively. Note here that this consideration does not change the incremental variational principle (50).

# 5. Numerical implementation and assessments

This section deals with the numerical implementation procedure for the gradient damage model proposed in Section 4. Since its total energy functional, i.e. Eq.(72), includes a quadratic term of the product between $\nabla u$ and $\alpha$, $I_n$ is hence separately but not globally convex with respect to them (see Bourdin et al. (2000) for more details). This makes the corresponding global minimization unachievable. Nevertheless, this fact does not exist between the displacement $u$ and the pressure $p$. Due to this reason, the variational optimization could be numerically carried out through a semi-staggered algorithm, which is achieved in the present work by using the open-source computing platform FEniCS (Alnaes et al., 2015). The whole procedure will be next applied to the hydraulic fracturing process, for which the model predictions will be compared to some available results for the purpose of assessment and validation.

## 5.1. Variational updates and numerical algorithm

We assume that the optimization of Eq.(44) in the current increment, i.e. $[t_n, t_{n+1}]$, can be sequentially achieved by first switching off the damage variation, such that $\alpha := \alpha_n$, and solving the following saddle-point problem where the poroelasticity and the fluid diffusion process are coupled:

$$ (\mathbf{u}_{n+1}, p_{n+1}) = \underset{\mathbf{u}'_{n+1}, p'_{n+1}}{\operatorname{Arg\,stat}} \left\{ I_n^{\text{poro}} (\mathbf{u}'_{n+1}, p'_{n+1}; \alpha_n) \right\} \quad (75) $$

By considering Eq.(72) and neglecting the constant terms, $I_n^{\text{poro}}$ can be explicitly expressed as:

$$ \begin{aligned} I_n^{\text{poro}}(\mathbf{u}_{n+1}, p_{n+1}; \alpha_n) = & \int_{\Omega} \left[ \frac{1}{2} \varepsilon_{n+1} : \mathbb{C}(\alpha_n) : \varepsilon_{n+1} - (p_{n+1} - p_0) \mathbf{b}(\alpha_n) : (\varepsilon_{n+1} - \varepsilon_n) \right. \\ & \left. - \frac{(p_{n+1} - p_n)^2}{2N(\alpha_n)} \right] d\Omega - \Delta t \int_{\Omega} \frac{1}{2} \nabla p_{n+1} \cdot \frac{\mathbf{k}}{\mu_f} \cdot \nabla p_{n+1} d\Omega \\ & - \int_{\Omega} \mathbf{F}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) d\Omega - \int_{\partial T \cap \Omega} \mathbf{T}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) dS \\ & - \Delta t \int_{\partial Q \cap \Omega} Q_{n+1} p_{n+1} dS \end{aligned} \quad (76) $$

We next proceed to estimate the damage solution $\alpha_{n+1}$ from the following minimizer:

$$ \alpha_{n+1} = \underset{\alpha'_{n+1}}{\operatorname{Arg\,min}} \left\{ I_n^{\text{dam}} (\alpha'_{n+1}) \right\} \quad (77) $$

where

$$ I_n^{\text{dam}}(\alpha_{n+1}) = \int_{\Omega} \left[ \frac{1}{2} \bar{\varepsilon}_{n+1} : \mathbb{C}(\alpha_{n+1}) : \bar{\varepsilon}_{n+1} - (\bar{p}_{n+1} - p_0) \mathbf{b}(\alpha_{n+1}) : (\bar{\varepsilon}_{n+1} - \varepsilon_n) \right. \\ \left. - \frac{(\bar{p}_{n+1} - p_n)^2}{2N(\alpha_{n+1})} + w_1 l_0 \left( \frac{\omega(\alpha_{n+1})}{l_0} + l_0 \| \nabla \alpha_{n+1} \|^2 \right) \right] d\Omega \quad (78) $$
where the bar superscript indicates the optimization of $u_{n+1}$ and $p_{n+1}$ that have been obtained from (75).

Recall that the permeability tensor $k$ (Eq.(64)) only parametrically depends on $\varepsilon$ and $\alpha$, its updates will be made at the end of each increment for the use of next one. Eq.(75) could be numerically achieved by making use of the *Nonlinear Variational Solver* module of the open-source DOLFIN (Logg et al., 2012) that allows an automatic iterations in each time increment with user-defined tolerances of absolute and relative errors as well as the maximal number of iterations, etc. Moreover, we employ the *mixed finite elements* which are particularly useful for the coupling problem of poroelasticity. For instance, the linear and quadratic Lagrange finite elements (i.e. CG1 and CG2 elements) are adopted to discretize the corresponding mixed function space that includes the fluid pore pressure and the displacement fields, respectively. The corresponding discussion is beyond the objective of this work, for which readers could be referred to Miehe and Mauthe (2016) for more details. In turn, we make use of the *PETSc TAOSolver* of DOLFIN for the damage optimization (see Eq.(77)) by primarily defining a corresponding error tolerance for the sake of computation convenience. Linear Lagrange finite elements (i.e. CG1 elements) are utilized for this minimization. It is also worthy to point out that this module also allows to numerically enforce the damage irreversibility condition, i.e. $\dot{\alpha} \ge 0$, by uploading its lower bound at each endpoint of increment, such that $\alpha \in [\alpha_{min}, 1]$ with $\alpha_{min} \leftarrow \alpha_{n+1}$ after its computation. Constant increment are adopted, such that $\Delta t = (t_f - f_0)/N$. Recall that $N$ is the total number of increments. The whole numerical procedure is detailed as below:
**Algorithm 1:** Semi-staggered algorithm proposed and implemented for the developed variational model in $[t_0, t_f]$

Input: {$\mathbf{u}_0, \alpha_0, p_0$} at $t_0$

Output: {$\mathbf{u}_n, \alpha_n, p_n$} at $t = t_1, \dots, t_n, t_{n+1}, \dots, t_N = t_f$

Set initial and boundary conditions;

Initialize $t = t_0$; set final time $t_f$, absolute tolerance $err_{tol}$ and calculate the constant time increment $\Delta t = \frac{t_f - t_0}{N}$;

Initialize the increment counter with $n = 0$;

while $t \le t_f$ do

$t = t_{n+1} = t_0 + (n + 1)\Delta t;$

Initialize an iteration counter with $i := 0$; set ($\mathbf{u}_{n+1}^0, p_{n+1}^0$) := ($\mathbf{u}_n, p_n$) and $\alpha = \alpha_n$;

Solve the saddle-point problem of poroelasticity coupled with fluid flow, i.e. Eq.(75), via NonlinearVariationalSolver from the following iterations in $[t_n, t_{n+1}]$;

repeat

Compute ($\mathbf{u}_{n+1}^{i+1}, p_{n+1}^{i+1}$) with ($\mathbf{u}_{n+1}^i, p_{n+1}^i$);

$i = i + 1$;

until $||\mathbf{u}_{n+1}^{i+1} - \mathbf{u}_{n+1}^i||_2 < err_{tol}$ and $||\mathbf{p}_{n+1}^{i+1} - \mathbf{p}_{n+1}^i||_2 < err_{tol}$;

Update $\mathbf{u}_n \leftarrow \mathbf{u}_{n+1}$ and $p_n \leftarrow p_{n+1}$;

Initialize another iteration counter with $j := 0$; set $\alpha_{n+1}^0 = \alpha_n$ and ($\mathbf{u}, p) := (\mathbf{u}_n, p_n)$;

Solve the minimization problem of gradient damage, i.e. Eq.(77), via

PETScTAOSolver from the following iterations in $[t_n, t_{n+1}]$;

repeat

Compute $\alpha_{n+1}^{i+1}$ with $\alpha_{n+1}^i$;

$j = j + 1$;

until $||\alpha_{n+1}^{i+1} - \alpha_{n+1}^i||_2 < err_{tol}$;

Update $\alpha_n \leftarrow \alpha_{n+1}$;

Update the lower bound of phase field $\alpha_{min} = \alpha_{n+1}$;

Update the permeability tensor $\mathbf{k}$ from Eq. (64).

end

## 5.2. Assessment and validation via application to hydraulic fracturing

In this section, the predictions of the proposed model in Section 4 is assessed by comparison to some available results which concern hydraulic fracturing. More specifically, we make use of the well-known AT2 model (Bourdin et al., 2000) for the description of the damage phenomena in the skeleton. By considering the equivalence between the energy dissipated by damage at failure and the material toughness (denoted by $G_c$), one has:

$$Y_0(\alpha) = 2w_1\alpha; \quad l_0 = \frac{G_c}{2w_1}; \quad \omega(\alpha) = \frac{1}{w_1} \int_0^\alpha 2w_1\beta \, d\beta \quad (79)$$

We consider a square plate (see Fig.2) that is embedded with a void at its center. Two preexisting cracks are considered at initial state on the top and bottom boundary of the void with the same length $L_{c0}$. The external surfaces of the plate are set with drained condition (i.e. $p = 0$). For
comparison purpose, we follow Mollaali et al. (2019) by adopting the same simulation procedure as well as the geometrical and material parameters in the context of CO₂ fracturing. Specifically, we make use of the so-called equivalent prescribed displacement that is numerically computed by subjecting the in-situ stress on the outside boundary of the specimen. The obtained displacements are next adopted as corresponding Dirichlet boundary conditions for the simulation. Another Neumman boundary condition, i.e. gas flux, denoted by $\bar{Q}_g$, is applied to the internal boundary of the central hole. It is worthy to note here that in Mollaali et al. (2019), the degradation of the Biot coefficient was described by following Chukwudozie et al. (2013) through $(1 - \alpha)b_0$, which might guarantee a good convergence especially when the fluid flux boundary conditions are applied. Due to this reason, we adopt in this test a similar treatment, but by calibrating the Biot coefficient through $[1 - (\phi_0 + \phi)]b(\alpha)$ with $b(\alpha)$ given in Eq. (54) and $\phi$ updated at the end of each increment. This would favour a good computation convergence without sacrificing the damage controlled degradation of the Biot coefficient.

The specimen dimension, the diameter of the central hole as well as the initial length of cracks are defined as $L = 170$ mm, $D = 20$ mm and $L_{c0} = 10$ mm, respectively. Note also that the characteristic size of elements is set as 0.5 mm. We investigate two cases by separately defining an isotropic stress loading $\sigma_1 = \sigma_3 = 1$ MPa and an anisotropic one $\sigma_1 = 3$ MPa and $\sigma_3 = 2$ MPa. The initial pressure in the cavity is set as $p_0 = 0.1$ MPa. The material parameters are reported in Table 1.

<figure><img src="image_3.png" /><figcaption>Figure 2: Geometry of the square plate with central void and predefined default cracks</figcaption></figure>

In Fig.3, we show the contours of the damage field $\alpha$, the horizontal displacement $u_x$ as well as the pressure distribution $p$ at a quasi-initial state when $t_0 \simeq 10^{-3}t_f$, two intermediate instants, e.g. $0.6 t_f$, $0.8 t_f$, and the final time $t_f$, respectively. Note here that this illustration is made only in the case of isotropic in-situ stresses because of its qualitatively similar contours with respect to the anisotropic ones. It can be observed in Fig.3(a) - 3(d) that the gas-driven fracture vertically propagate with a detectable horizontal gradient that could be related to the crack opening displacements as qualitatively captured in Fig.3(i) - 3(l). This fact is described in the proposed model through the parameter $\eta$ that also significantly affects the permeability in the damage zone (see Eq.(64)). Also, the gas pressure horizontally displays an important gradient, as shown in Fig.3(e) - 3(h), that significantly changes in return the crack opening.
<figure><img src="image_4.png" /><figcaption>Figure 3: Hydraulic fracturing process of the square plate with a central void and default cracks under isotropic in-situ stress $\sigma_1 = \sigma_3 = 1MPa$. (a)-(d): damage field; (e)-(h): pressure field; (i)-(l): horizontal displacement.</figcaption></figure>
<table><caption>Table 1: Material parameters for the simulation of CO₂ fracturing reported by Wang et al. (2018) and Mollaali et al. (2019)</caption><thead><tr><th>Parameter</th><th>Name</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td><i>E</i></td><td>Young's modulus</td><td>6 &times; 10<sup>3</sup></td><td>MPa</td></tr><tr><td><i>&nu;</i></td><td>Poisson's ratio</td><td>0.34</td><td>-</td></tr><tr><td><i>G</i><sub>c</sub></td><td>Critical fracture energy</td><td>0.30</td><td>MPa &middot; mm</td></tr><tr><td><i>b</i><sub>0</sub></td><td>Initial Biot coefficient</td><td>0.85</td><td>-</td></tr><tr><td><i>k</i><sup>s</sup></td><td>Permeability of healthy rock</td><td>1 &times; 10<sup>-12</sup></td><td>mm<sup>2</sup></td></tr><tr><td><i>&mu;</i><sup>f</sup></td><td>Dynamic viscosity</td><td>4.0 &times; 10<sup>-11</sup></td><td>MPa &middot; s</td></tr><tr><td><i>l</i><sub>0</sub></td><td>Internal length of damage</td><td>1.60</td><td>mm</td></tr><tr><td><i>&sigma;</i><sub>T</sub></td><td>Rock's tensile strength</td><td>11</td><td>MPa</td></tr></tbody></table>

Furthermore, in order to carry out a quantitative assessment, we first pay attention to the breakdown pressure, i.e. pressure in the void at the onset of crack propagation, denoted by $p_{br}$. The related analytical solution was derived by Haimson et al. (1967) and reads:

$$p_{br} = \frac{3\sigma_3 - \sigma_1 + \sigma_T + p_0}{1 + \frac{v}{1-v}b} \qquad (80)$$

It is worth noting here that for the adopted AT2 model adopted in the present study, the tensile strength $\sigma_T$ and $w_1$ are related by $\sigma_T^2 = (\frac{8\sqrt{2}}{3\sqrt{3}})^2 E w_1$ (see for instance Marigo et al. (2016)). Simultaneously considering Eq. (79) yields $l_0 = 3G_c E/(16\sigma_T/3)^2$, where the parameters $E$, $G_c$, $l_0$ and $\sigma_T$ provided in Table 1 are in perfect agreement.

We display in Table 2 the corresponding numerical estimation obtained from the proposed model, the analytical prediction of (see Eq. (80)), some reference computational results obtained by Mollaali et al. (2019) as well as the experimental measurement of Ishida et al. (2012). It can be observed that the model predictions are in good agreement with those of Mollaali et al. (2019), but both of them overestimate the analytical and experimental results.

For completeness, let us introduce, following Miehe et al. (2010a) (see also Miehe et al. (2010b)) the so-called *regularized crack functional*<fnref n="8" />:

$$\mathcal{L}_c = \frac{1}{2l_0} \int_{\Omega} \alpha^2 + l_0 \| \nabla \alpha \|^2 d\Omega \qquad (81)$$

We show in Fig. 4 the evolution of its normalized value, i.e. $\mathcal{L}_c/L$, as well as the pressure in the void (denoted by $p_v$) before and during the damage propagation. These predictions are also compared to those obtained by Mollaali et al. (2019). A good agreement can be found for the

<footnote n="8">⁸Practically, Eq. (81) is called the normalized dissipation energy of the AT2 model since it is equivalent to $\mathcal{L}_c = \Phi^s/G_c$ by considering (79).</footnote>
<table><thead><tr><th>In situ stress</th><th>Numer.<sup>a</sup></th><th>Theoretical<br>Haimson et al. (1967)</th><th>Ref. Simulation<sup>b</sup></th><th>Experiments<br>Ishida et al. (2012)</th></tr></thead><tbody><tr><td>Isotropic case</td><td>10.66</td><td>9.10</td><td>10.90</td><td>8.44</td></tr><tr><td>Anisotropic case</td><td>11.63</td><td>9.80</td><td>11.90</td><td>-</td></tr></tbody></table>

<sup>a</sup> Numerical predictions of the proposed model

<sup>b</sup> Numerical results of Mollaali et al. (2019) available in literature

<caption>Table 2: Numerical estimations of the breakdown pressure $p_{tb}$ (MPa) and the comparison to the analytical results of Haimson et al. (1967), the numerical computations of Mollaali et al. (2019) as well as to the experimental observation of Ishida et al. (2012) in the case of the isotropic in-situ stresses $\sigma_1 = \sigma_3 = 1$ MPa and that of anisotropic ones $\sigma_1 = 3$ MPa and $\sigma_3 = 2$ MPa, respectively.</caption>

entire loading history. Another interesting observation is related to the pressure value when these curves simultaneously change their slope. This indicates the onset of the fracture propagation when $p_v = p_{br}$.

## 6. Predictions for excavation damaged zone

This section is devoted to a numerical study of the so-called excavation damaged zone (EDZ) that is an important issue in geomechanical engineering. This will be achieved by applying the proposed variational model to the structural computation of an underground gallery with the host Callovo-Oxfordian clay rocks. The lattices are selected by the French authorities for the study of the feasibility of storage of high-level, long-lived radioactive waste. In practice, the depth of these galleries could reach about 500 m, yet, the EDZ is usually generated just in a relatively limited zone (Armand et al., 2013, 2017). We hence make use here of a square model, as illustrated in Fig.5, with 50 m of size. Moreover, the galleries are approximated as in a circular form that, owing to its symmetry, could be schematized from a quarter of its geometry. The characteristic size of elements is set as $0.04m$. According to the experimental works of Armand et al. (2013, 2014), the whole structure is assumed to be subjected to a water pressure $p_0 = 4.7$ MPa and to a vertical and horizontal compression in-situ stresses, denoted by $\sigma_v$ and $\sigma_h$, respectively. We hence make use of the model incorporating the unilateral effect, by means of the thermodynamic potential expressed in Eq.(74). In addition, we proposed here to carry out a computation in the case of isotropic loading ($\sigma_v = \sigma_h = 12.7$MPa). The material parameters are summarized in Table 3 which can also be found in Armand et al. (2013, 2014), except for that of the fracture energy $G_c$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Numerical predictions of the normalized crack length $L_c/L$ and the pressure in the void $p_v$ as well as their comparison to the computations realized by Mollaali et al. (2019). (a): Isotropic case; (b): Anisotropic case.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Geometric model and boundary conditions of the underground gallery</figcaption></figure>

Figs.6 show the contour of the damage evolution. Attentions must be payed to its colorbar when $t \le 500$ days, whose values are limited in an interval, such that $0 \le \alpha \le 0.5$ in order to appropriately illustrate the onset and propagation of the damage field$^9$. We observe that at the beginning, e.g. $t = 6$ days, an instantaneous but slight EDZ could be already produced with the maximal value of $\alpha$ being almost 0.2. As time progresses, e.g. until $t = 500$ or $600$ days, the damage field shows localization zone near the gallery wall. Two fractures with $\alpha \approx 1$ are then nucleated (see Fig.6(e)), which are symmetrically located near the top and the bottom of the gallery. Each of

<footnote>$^9$Note here that at $t = 500$ days, the maximal value of $\alpha$ is computed as being almost equal to 0.49.</footnote>
<table><caption>Table 3: Material parameters of the Callovo-Oxfordian clay rocks for the storage of radioactive wastes by Armand et al. (2013, 2014) except for that of $G_c$.</caption><thead><tr><th>Parameter</th><th>Name</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td><i>E</i></td><td>Young's modulus</td><td>5 × 10<sup>3</sup></td><td>MPa</td></tr><tr><td><i>ν</i></td><td>Poisson's ratio</td><td>0.295</td><td>-</td></tr><tr><td><i>G</i><sub>c</sub></td><td>Critical fracture energy</td><td>0.275</td><td>MPa·mm</td></tr><tr><td><i>b</i><sub>0</sub></td><td>Initial Biot coefficient</td><td>0.6</td><td>-</td></tr><tr><td><i>k</i><sup>s</sup></td><td>Permeability of healthy rock</td><td>6 × 10<sup>-20</sup></td><td>m<sup>2</sup></td></tr><tr><td><i>μ</i><sup>f</sup></td><td>Dynamic viscosity</td><td>1.0 × 10<sup>-10</sup></td><td>MPa·s</td></tr><tr><td><i>l</i><sub>0</sub></td><td>Internal length of damage</td><td>10</td><td>mm</td></tr><tr><td><i>φ</i><sub>0</sub></td><td>Initial porosity</td><td>0.16</td><td>-</td></tr></tbody></table>

them, as shown in Fig. 6(f), propagates towards another one. Note here that very slight asymmetry between these two fractures could be noted. The reason may probably be searched in the fact that the numerical solution issued from the minimization procedure can bifurcate from the fundamental branch; this can lead to some loss of symmetry.
<figure><img src="image_7.png" /><figcaption>Figure 6: Contour of the damage field $\alpha$ during the evolution of EDZ</figcaption></figure>

In Fig.7, we display the profile of pressure field during the evolution of EDZ. By simultaneously examining Fig.6, it can be observed that even the damage field is less important at initial times, e.g. $t = 6$ days, the water pressure significantly decreases near the edge of the gallery in a non-uniform manner. This should be due to the permeability description that is affected by the damage state (see Eq.(64)). During the evolution of EDZ, the pressure profile is shown to have a gradient distribution near and in the fracture zone.
<figure><img src="image_8.png" /><figcaption>Figure 7: Contour of the pressure field $p$ during the evolution of EDZ</figcaption></figure>

For a more comprehensive visualization, we show in Fig.8 the displacements evolution and that of the water pressure at some targeted points on the bottom boundary of the geometric model. As illustrated in Fig.5, these points are sequentially denoted as $Pt_i$ with $i \in [1, 2, 3, 4, 5]$, being 0.05m, 1.10m, 1.90m, 4.80m and 10.00m away from the drifts. It can be observed from Figs.8(a) that at the beginning, e.g. $t \mapsto 0$ day, the horizontal displacements $u_x$ temporally increase until about $t = 50 - 100$ days. This should be due to the slight damage evolution at this initial period on the drift edge that yet allows a relatively noticeable modification of the elasticity modulus (c.f. Eq.(52)) as well as that of the permeability. The latter change can also be captured, as shown in Figs.8(b), by a remarkable decrease of the water pressure during this period. By contrast, their evolution then becomes somewhat negligible during a long time, even as shown in Figs.8(b), the decrease of the water pressure is relatively a little more important. After about $t = 700$ days (as shown in Fig.6(e) the damage already reaches $\alpha = 1$), a significant change of the slope of the displacement evolution at $Pt_1$ is observed. This is probably due to an effect of crack nucleation
and growth.

<figure><img src="image_9.png" /><figcaption>Figure 8: Evolution of displacements and pressure on the bottom boundary of the geometry model, represented by the points $Pt_i$ with $i \in [1, 2, 3, 4, 5]$, respectively. (a): horizontal displacement $u_x$; (b): water pressure $p$.</figcaption></figure>

## 7. Conclusions

In this work, we have developed an incremental variational approach for saturated porous media that allows to incorporate the full coupling between the poroelasticity, the fluid flow as well as the damage and its gradient. The developed non-local formulation relies on the Generalized Standard Materials framework extended here to poroelasticity together with gradient damage. In this context, three points must be pointed out.

* The proposed incremental variational principle consists in the minimization of a four-fields incremental energy functional. The consideration of a *kinetic porosity*, required for an appropriate description of the coupling between the fluid transport (via Darcy law) and the deformation of the skeleton, has allowed to reduce the above minimization to that of a three-fields functional.

* By taking advantage of this variational approach, we have formulated a complete poroelastic-gradient damage model, whose numerical implementation has been carried out through a semi-staggered algorithm by using the open-source computing platform FEniCS.

* For assessment and application purpose of the proposed model, two main tests have been presented: the first one concerned the propagation of hydraulic fractures, while the second one is devoted to the investigation of an Excavation Damage Zone around a waste storage underground gallery.

In detail, the irreversible process is quantified via a dissipation potential that are expressed in a joint form between the damage dissipation and the fluid diffusion. This permits a relevant description of the coupled hydromechanical aspects of quasi brittle saturated media. The thermodynamics descriptions lead to the construction of a total energy functional, whose incremental
form is dependent on the displacement and damage fields of the skeleton phase, the fluid pressure as well as on the kinetic porosity. It was then demonstrated that the corresponding weak forms lead to the Euler-Lagrange equations (i.e. governing equations) and to the boundary conditions for the fully coupled system. An incremental variational principle was hence proposed. Moreover, as already mentioned, this variational principle has been reduced as three variables dependent based on an adequate consideration of the kinetic porosity in the time increment.

Practical application of this incremental variational approach has been done by considering the linear Darcy flow for the time dependent fluid transport. Accordingly, after numerically implementing the incremental variational approach, we proposed its assessment through a first application ting in an hydraulic fracturing problem, for which the predictions are well compared to some available results. This numerical procedure is then applied to the study of an EDZ growth until occurrence and growth of fractures around an underground excavation involved in a radioactive waste storage project in France.

In addition to all the foregoing, it should be noted that the proposed variational approach can be extended to other types of coupled behaviors involving non-local damage process. For instance, the consideration of the thermoporoelasticity and the poroplasticity coupled with gradient damage may also constitute important extensions, which could pave the way for better and more realistic predictions in different engineering problems.

## Data availability

Data will be made available on request.

## Acknowledgment

This work was funded by the China Scholarship Council (No. 201906450075). L.C. and D.K. would also like to acknowledge partial supports from the French National Research Agency through the ANR project BigBen (grant ANR-22-CE51-0034-02) and the NEEDS (Nucléaire: Energie, Environnement, Déchets, Société) through the exploratory project MENUFrac-THM. L.C. thanks also Géry de Saxcé and Fabrice Golfier for the discussions and comments that allowed to improve the paper. Finally, X-D.Z. would like to express his thanks to Vahid Ziaei-Rad for fruitful exchanges.

## Appendix A. Local form of thermodynamics principles for saturated porous media

Let us first look at the energy balance equation. In Eq.(1), the rate of global internal energy can be expressed as (see for instance Eq.(3.15) in Coussy (2004)):

$$ \dot{\mathcal{E}} = \int_{\Omega} \left[ \frac{d^s e}{dt} + e \operatorname{div}(\mathbf{v}^s) + \operatorname{div}(e^f \rho^f \mathbf{v}) \right] d\Omega \quad (A.1) $$

where the operator $\frac{d^s a}{dt} = \frac{\partial a}{\partial t} + \nabla a \cdot \mathbf{v}^s$ determines the particle derivative of a scalar field a with respect to the skeleton particle.

The calorific power can be classically written as:

$$ \mathcal{P}_{cal} = \int_{\Omega} r \, d\Omega - \int_{\partial\Omega} \mathbf{q} \cdot \mathbf{n} \, dS \quad (A.2) $$
with $n$ being the unit outward normal on the boundary of $d\Omega$.

We neglect the gravity and inertial forces both of the skeleton and fluid parts, the power of internal work could be given by (see also Eq.(2.45) in Coussy (2004) for the so-called total strain work rate):

$$P_{int} = - \int_{\Omega} [\boldsymbol{\sigma} : \boldsymbol{d}^s - \mathrm{div}(p\boldsymbol{v})] \, d\Omega \quad (A.3)$$

where $d^s$ symbolizes the local strain rate of the skeleton and $p$ the fluid pressure.

Inserting Eqs.(A.1) - (A.3) into (1), one gets the following local form of the first principle of poroelasticity:

$$\frac{d^s e}{dt} + e\mathrm{div}(v^s) = \boldsymbol{\sigma} : d^s + r - \mathrm{div}(q) - \mathrm{div}\left[\rho^f \boldsymbol{v} \left(e^f + \frac{p}{\rho^f}\right)\right] \quad (A.4)$$

for which readers could also be referred to Eq.(3.19) of Coussy (2004) for a similar expression without the consideration of $r$ but with that of the gravity and inertia forces. Taking into account the infinitesimal transformation of the skeleton, Eq.(A.4) could be reduced to Eq.(3).

Coming now to the second principle of the concerned media and following a similar calculation of (A.1), the global entropy can be calculated as:

$$\dot{S} = \int_{\Omega} \left[ \frac{d^s s}{dt} + s\mathrm{div}(v^s) + \mathrm{div}(\rho^f s^f v) \right] d\Omega \quad (A.5)$$

Putting Eq.(A.5) into (2) leads to the local form of entropy inequality, which reads:

$$\dot{s} + s\mathrm{div}(v^s) \geq \frac{r}{T} - \mathrm{div}\left(\frac{q}{T}\right) - \mathrm{div}\left(\rho^f s^f v\right) \quad (A.6)$$

Again, in the case of infinitesimal transformations, Eq.(A.6) takes the following form:

$$\dot{s} \geq \frac{r}{T} - \mathrm{div}\left(\frac{q}{T}\right) - \mathrm{div}\left(\rho^f s^f v\right) \quad (A.7)$$

which is a counterpart of Eq.(3.23) in Coussy (2004) without the consideration of $r$.

Finally, the combination of Eqs.(3) and (A.7) yields the Clausius - Duhem inequality, i.e. Eq.(4).

## Appendix B. Variational analysis of reversible process in poroelasticity

This section is devoted to a variational analysis on the reversible process in poroelasticity. Inspired from Heuzé and Stainier (2022), we make use of the so-called residual of the rate of Legendre transform (see Eq.(16)), which allows to define the following Lagrangian function:

$$\chi(\dot{\phi}, \dot{p}, \lambda) = \dot{w}^s - p\dot{\phi} + \lambda \frac{d}{dt}(p\phi + \psi^s - w^s) \quad (B.1)$$

where $\lambda$ is the Lagrange multiplier. Note that as we would like to search for the stationary action that is mainly described by $\dot{\phi}$ and $\dot{p}$, the rate of Helmholtz free energy $\dot{w}^s$, the displacement rate $\dot{u}$ as well as the non-reversible variables (e.g. those defining the fluid diffusion and gradient damage) are considered as known here.
The stationarity of Eq.(B.1) with respect to the non-dissipative arguments (i.e. $\phi$ and $p$) and the Lagrangian multiplier $\lambda$ gives:

$$ \left\{ \begin{array}{l} -p + \lambda p = 0 \\ \lambda \left( \phi + \frac{\partial \psi^s}{\partial p} \right) = 0 \\ \frac{d}{dt} (p\phi + \psi^s - w^s) = 0 \end{array} \right. \qquad (B.2) $$

which yields $\lambda = 1$, $\phi = -\frac{\partial \psi^s}{\partial p}$ (see also Eq.(19)) as well as the rate form of Legendre transform. It can be readily calculated that under these conditions, Eq.(B.1) can be recast into:

$$ \underset{\phi, \dot{p}, \lambda}{\text{stat}} \chi = \chi(\phi, p; \lambda = 1) = \dot{\psi}^s + \dot{p}\phi \qquad (B.3) $$

Likewise, we could also suppose the rate of potential function $\dot{\psi}^s$ to be known, which allows to the definition of the following Lagrangian function by introducing the Lagrange multiplier $\lambda'$:

$$ \chi'(\dot{\phi}, \dot{p}, \lambda') = \dot{\psi}^s + \dot{p}\phi + \lambda' \frac{d}{dt} (p\phi + \psi^s - w^s) \qquad (B.4) $$

In the same way, the stationarity with respect to $\dot{p}$ and $\phi$ as well as to $\lambda'$ leads to

$$ \left\{ \begin{array}{l} \phi + \lambda' \phi = 0 \\ \lambda' \left( p - \frac{\partial w^s}{\partial \phi} \right) = 0 \\ \frac{d}{dt} (p\phi + \psi^s - w^s) = 0 \end{array} \right. \qquad (B.5) $$

allowing to recast Eq.(B.4) into:

$$ \underset{\phi, \dot{p}, \lambda'}{\text{stat}} \chi' = \chi'(\phi, p; \lambda' = -1) = \dot{w}^s - p\dot{\phi} \qquad (B.6) $$

Apparently, both of Eqs.(B.1) and (B.4) are able to describe the reversible power in the context of poroelasticity. In the context of thermodynamics, both of them could be eventually coupled with some dissipative process, such as the fluid diffusion and the gradient damage. This is carried out in Section 3 by making use of an incremental variational analysis for saturated porous media.

## References

* Alnaes, M.S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N., 2015. The FEniCS project version 1.5. Archive of Numerical Software 3.

* Amor, H., Marigo, J.J., Maurini, C., 2009. Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. Journal of the Mechanics and Physics of Solids 57, 1209–1229.

* Armand, G., Conil, N., Talandier, J., Seyed, D.M., 2017. Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: From experimental studies to model calibration and validation. Computers and Geotechnics 85, 277–286.

* Armand, G., Leveau, F., Nussbaum, C., de La Vaissiere, R., Noiret, A., Jaeggi, D., Landrein, P., Righini, C., 2014. Geometry and properties of the excavation-induced fractures at the meuse/haute-marne url drifts. Rock Mech Rock Eng 47, 21–41.

* Armand, G., Noiret, A., Zghondi, J., Seyed, D., 2013. Short- and long-term behaviors of drifts in the callovo-oxfordian claystone at the meuse/haute-marne underground research laboratory. Journal of Rock Mechanics and Geotechnical Engineering 5, 221–230.
Bazant, Z.P., Pijaudier-Cabot, G., 1988. Nonlocal continuum damage, localization instability and convergence. Journal of Applied Mechanics 55, 287-293.

Biot, M., 1955. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics 26, 182-185.

Biot, M.A., 1941. General theory of three-dimensional consolidation. Journal of applied physics 12, 155-164.

Bourdin, B., Francfort, G.A., Marigo, J.J., 2000. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48, 797-826.

Bourdin, B., Francfort, G.A., Marigo, J.J., 2008. The variational approach to fracture. Journal of elasticity 91, 5-148.

Cheng, A.H.D., 2016. Poroelasticity. Springer.

Chukwudozie, C., Bourdin, B., Yoshioka, K., 2013. A variational approach to the modeling and numerical simulation of hydraulic fracturing under in-situ stresses. Proceedings of the 38th Workshop on Geothermal Reservoir Engineering.

Chukwudozie, C., Bourdin, B., Yoshioka, K., 2019. A variational phase-field model for hydraulic fracturing in porous media. Comput. Methods Appl. Mech. Engrg 347, 957-982.

Coleman, B.D., Gurtin, M.E., 1967. Thermodynamics with internal state variables. J. Chem. Phys. 47, 597-613.

Coussy, O., 1995. Mechanics of porous media. Chichester: Wiley.

Coussy, O., 2004. Poromechanics. Wiley, New Jersey.

Dormieux, L., Kondo, D., 2004). Approche micromécanique du couplage perméabilité-endommagement. C. R. Mecanique 332, 135-140.

Dormieux, L., Kondo, D., 2005. Poroelasticity and damage theory for saturated cracked media. In: Dormieux, L., Ulm, F.-J. (Eds.), Applied Micromechanics of Porous Materials. CISM Courses and Lectures, Springer, Wien, New York 480, 207-288.

Dormieux, L., Kondo, D., Ulm, F.J., 2006. Microporomechanics. John Wiley and Sons, Ltd, West Sussex, England.

Dormieux, L., Molinari, A., Kondo, D., 2002. Micromechanical approach to the behavior of poroelastic materials. Journal of the Mechanics and Physics of Solids 50(10), 2203-2231.

Ehrlacher, A., Fedelich, B., 1989. Stability and bifurcation of simple dissipative systems ; application to brutal damage. cracking and damage : Strain localization and size effect. edited by J. Mazars and Z.P.Bazand, Taylor & Francis, 221-227.

Frémond, M., Nedjar, B., 1996. Damage, gradient of damage and principle of virtual power. International journal of solids and structures 33, 1083-1103.

Frémond, M., 2002. Non-Smooth Thermomechanics. Springer-Verlag Berlin Heidelberg New York.

Frémond, M., Nedjar, B., 1993. Endommagement et principe des puissances virtuelles. Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre 317(7), 57-864.

Frémond, M., Nedjar, B., 1996. Damage, gradient of damage and principle of virtual power. International Journal of Solids and Structures 33(8), 1083-1103.

Gao, Y., Detournay, E., 2020. A poroelastic model for laboratory hydraulic fracturing of weak permeable rock. Journal of the Mechanics and Physics of Solids 143, 104090.

Guéguen, Y., Palciauskas, V., 1994. Introduction to Physics of Rocks. Princeton Univ. Press.

Haimson, B., Fairhurst, C., et al., 1967. Initiation and extension of hydraulic fractures in rocks. Society of Petroleum Engineers Journal 7, 310-318.

Halphen, B., Nguyen, Q.C., 1975. Sur les matériaux standards généralisés. Journal de Mécanique 14(1), 39-63.

Heuzé, T., Stainier, L., 2022. A variational formulation of thermomechanical constitutive update for hyperbolic conservation laws. Comput. Methods Appl. Mech. Engrg. 394, 114893.

Ishida, T., Aoyagi, K., Niwa, T., Chen, Y., Murata, S., Chen, Q., Nakayama, Y., 2012. Acoustic emission monitoring of hydraulic fracturing laboratory experiment with supercritical and liquid co2. Geophysical research letters 39, L16309.

Kachanov, L.M., 1958. Time of the rupture process under creep conditions. Izv Akad. Nauk USSR Otd Tekh. 8, 26-31.

Kpotufe, K., 2021. Modèles à gradient d'endommagement : cadre thermodynamique, formulation variationnelle et applications. Ph.D. thesis. Sorbonne Université et Université de Lomé.

Lemaitre, J., Chaboche, J., 1978. Aspect phénoménologique de la rupture par endommagement. J. Mécanique Appliquée 2, 3, 317-365.

Logg, A., Mardal, K., Wells, G.N.e.a., 2012. Automated Solution of Differential Equations by the Finite Element Method. Springer.

Lorentz, E., Andrieux, S., 1999. A variational formulation for nonlocal damage models. International Journal of
Plasticity 15(2), 119–138.

Lorentz, E., Benallal, A., 2005. Gradient constitutive relations: numerical aspects and application to gradient damage. Computer Methods in Applied Mechanics and Engineering 194, 5191–5220.

Maitournam, H., 2013. Mécanique des structures anélastique. Editions de l'Ecole Polytechnique.

Marigo, J.J., 1981. Formulation d'une loi d'endommagement d'un matériau élastique. Comptes Rendus Académie des Sciences, Serie II 292, 1309 – 1312.

Marigo, J.J., Maurini, C., Pham, K., 2016. An overview of the modelling of fracture by gradient damage models. Meccanica 51, 3107–3128.

Martínez-Pañeda, E., Golahmarb, A., Niordsonb, C.F., 2018. A phase field formulation for hydrogen assisted cracking. Computer Methods in Applied Mechanics and Engineering 342, 742 – 761.

Mauthe, S., Miehe, C., 2017. Hydraulic fracture in poro-hydro-elastic media. Mechanics Research Communications 80, 69–83.

Miehe, C., Hofacker, M., Welschinger, F., 2010a. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199, 2765–2778.

Miehe, C., Mauthe, S., 2016. Phase field modeling of fracture in multi-physics problems. part iii. crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media. Computer Methods in Applied Mechanics and Engineering 304, 619–655.

Miehe, C., Mauthe, S., Teichtmeister, S., 2015. Minimization principles for the coupled problem of darcy–biot-type fluid transport in porous media linked to phase field modeling of fracture. Journal of the Mechanics and Physics of Solids 82, 186–217.

Miehe, C., Welschinger, F., Hofacker, M., 2010b. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. International journal for numerical methods in engineering 83, 1273–1311.

Mollaali, M., Ziaei-Rad, V., Shen, Y., 2019. Numerical modeling of co2 fracturing by the phase field approach. Journal of Natural Gas Science and Engineering 70, 102905.

Moreau, J.J., 1974. On Unilateral Constraints, Friction and Plasticity. In New Variational Techniques in Mathematical Physics, ed. G. Capriz and G. Stampacchia, Centro Internazionale Matematico Estivo, II Ciclo, Edizioni Cremonese, Roma.

Nguyen, Q.S., 2016. Quasi-static responses and variational principles in gradient plasticity. Journal of the Mechanics and Physics of Solids 97, 156–167.

Noii, N., Khodadadian, A., Wick, T., 2021. Bayesian inversion for anisotropic hydraulic phase-field fracture. Comput. Methods Appl. Mech. Engrg 386, 114118.

Ortiz, M., Martin, J., 1989. Symmetry - preserving return mapping algorithms and incrementally external paths: a unification of concepts. Int. J. Numer. Methods Eng. 28, 1839–1853.

Ortiz, M., Repetto, E., 1999. Nonconvex energy minimization and dislocation structures in ductile single crystals. J. Mech. Phys. Solids 47 (7-8), 379–462.

Ortiz, M., Stainier, L., 1999. The variational formulation of viscoplastic constitutive updates. Advances in Applied Mechanics 171, 419–444.

Peerlings, R., de Borst, R., Brekelmans, W., Vree, J.D., 1996. Gradient enhanced damage for quasi-brittle materials. International Journal for Numerical Methods in Engineering 39, 3391 – 3403.

Pham, K., Amor, H., Marigo, J.J., Maurini, C., 2011. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics 20, 618–652.

Pillai, U., Heider, Y., Markert, B., 2018). A diffusive dynamic brittle fracture model for heterogeneous solids and porous materials with implementation using a user-element subroutine. Comput. Mater. Sci. 153, 36–47.

Rice, J.R., 1971. Inelastic constitutive relations for solids: An internal-variable theory and its application to metal plasticity. Journal of the Mechanics and Physics of Solids 19, 433 – 455.

Rice, J.R., Cleary, M.P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Reviews of Geophysics 14(2), 227–241.

Stainier, L., 2013. Chapter two - a variational approach to modeling coupled thermo-mechanical nonlinear dissipative behaviors. Advances in Applied Mechanics 46, 69–126.

Ulloa, J., Noii, N., Alessi, R., Aldakheel, F., Degrande, G., François, S., 2022. Variational modeling of hydromechanical fracture in saturated porous media: A micromechanics-based phase-field approach. Comput. Methods Appl. Mech. Engrg 396, 115084.

Wang, J., Elsworth, D., Wu, Y., Liu, J., Zhu, W., Liu, Y., 2018. The influence of fracturing fluids on fracturing processes: a comparison between water, oil and sc-co2. Rock Mechanics and Rock Engineering 51, 299–313.
Wu, J., Nguyen, V., Nguyen, C., Sutula, D. adn Sina, S., P.A.B, S., 2020. Chapter one - phase-field modeling of fracture. Advances in Applied Mechanics 53, 1-183.

Xia, L., Yvonnet, J., Ghabezloo, S., 2017. Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media. Engineering Fracture Mechanics 186, 158 – 180.

Yang, Q., Stainier, L., Ortiz, M., 2006. A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids. Journal of the Mechanics and Physics of Solids 54, 401 – 424.

Yoshioka, K., Bourdin, B., 2016. A variational hydraulic fracturing model coupled to a reservoir simulator. International Journal of Rock Mechanics and Mining Sciences 88, 137–150.

Zhu, W., Wei, C., Li, S., Wei, J., Zhang, M., 2013. Numerical modeling on destress blasting in coal seam for enhancing gas drainage. International Journal of Rock Mechanics and Mining Sciences 59, 179–190.

Ziegler, H., 1963. Some extremum principles in irreversible thermodynamics with applications to continuum mechanics. in: I.N. Sneddon and R. Hill, eds., Progress in Solid Mechanics 4 (North-Holland, Amsterdam).