# Phase transition in saturated porous media: Pore-fluid segregation in consolidation

Emilio N.M. Cirillo, Nicoletta Ianiro, Giulio Sciarra

► To cite this version:

Emilio N.M. Cirillo, Nicoletta Ianiro, Giulio Sciarra. Phase transition in saturated porous media: Pore-fluid segregation in consolidation. Physica D: Nonlinear Phenomena, 2011, 240 (17), 10.1016/j.physd.2011.05.015. hal-01362173

HAL Id: hal-01362173

https://hal.science/hal-01362173

Submitted on 3 Jun 2024

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 Attribution 4.0 International License
# Phase transition in saturated porous media: Pore-fluid segregation in consolidation

Emilio N.M. Cirillo, Nicoletta Ianiro, Giulio Sciarra

► To cite this version:

Emilio N.M. Cirillo, Nicoletta Ianiro, Giulio Sciarra. Phase transition in saturated porous media: Pore-fluid segregation in consolidation. Physica D: Nonlinear Phenomena, 2011, 240 (17), 10.1016/j.physd.2011.05.015. hal-01362173

HAL Id: hal-01362173

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

Submitted on 3 Jun 2024

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 Attribution 4.0 International License
# Phase transition in saturated porous media: pore-fluid segregation in consolidation

Emilio N.M. Cirillo¹, Nicoletta Ianiro¹ and Giulio Sciarra²,*

¹ Dipartimento Me. Mo. Mat., Università di Roma “La Sapienza”, via A. Scarpa 16, 00161 Rome, Italy

² Dipartimento Ingegneria Chimica Materiali Ambiente, Università di Roma “La Sapienza”, via Eudossiana 18, 00184 Rome, Italy

April 27, 2022

## Abstract

Consider the consolidation process typical of soils, this phenomenon is expected not to exhibit a unique state of equilibrium, depending on the *external loading* and the constitutive parameters. Beyond the standard solution, also pore-fluid segregation, which is typically associated with fluidization of the granular material, can arise. Pore-fluid segregation has been recognized as a phenomenon typical of the short time behavior of a saturated porous slab or a saturated porous sphere, during consolidation. In both circumstances Biot's three dimensional model provides time increasing values of the water pressure (and fluid mass density) at the center of the slab (or of the sphere), at early times, if the Lamé constant $\mu$ of the skeleton is different from zero. This localized pore-fluid segregation is known in the literature as Mandel-Cryer effect. In this paper a non linear poromechanical model is formulated. The model is able to

<footnote>*Corresponding author: e-mail giulio.sciarra@uniroma1.it Tel. +390644585230 Fax +390644585618</footnote>
# Phase transition in saturated porous media: pore-fluid segregation in consolidation

Emilio N.M. Cirillo¹, Nicoletta Ianiro¹ and Giulio Sciarra²,*

¹ Dipartimento Me. Mo. Mat., Università di Roma “La Sapienza”, via A. Scarpa 16, 00161 Rome, Italy

² Dipartimento Ingegneria Chimica Materiali Ambiente, Università di Roma “La Sapienza”, via Eudossiana 18, 00184 Rome, Italy

April 27, 2022

## Abstract

Consider the consolidation process typical of soils, this phenomenon is expected not to exhibit a unique state of equilibrium, depending on the *external loading* and the constitutive parameters. Beyond the standard solution, also pore-fluid segregation, which is typically associated with fluidization of the granular material, can arise. Pore-fluid segregation has been recognized as a phenomenon typical of the short time behavior of a saturated porous slab or a saturated porous sphere, during consolidation. In both circumstances Biot's three dimensional model provides time increasing values of the water pressure (and fluid mass density) at the center of the slab (or of the sphere), at early times, if the Lamé constant $\mu$ of the skeleton is different from zero. This localized pore-fluid segregation is known in the literature as Mandel-Cryer effect. In this paper a non linear poromechanical model is formulated. The model is able to

<footnote>*Corresponding author: e-mail giulio.sciarra@uniroma1.it Tel. +390644585230 Fax +390644585618</footnote>
describe the occurrence of two states of equilibrium and the switching from one to the other by considering a kind of *phase transition*. Extending classical Biot's theory a more than quadratic strain energy potential is postulated, depending on the strain of the porous material and the variation of the fluid mass density (measured with respect to the skeleton reference volume). When the consolidating pressure is strong enough the existence of two distinct minima is proven.

**Keywords:** Granular materials, Porous media, Phase transition, Bifurcation theory

# 1 Introduction

Porous media confined into a fluid infinite reservoir, or suffering consolidation loading, can exhibit solid-solid and solid-fluid phase transitions. These last can be observed because of different phenomena, as gravity driven solid-fluid separation of suspensions (see Bürger et al., 2000) or solid-fluid segregation of soils (see e.g. Nichols et al., 1994; Vardoulakis, 2004a,b). It will be the purpose of this paper to investigate this last transition.

The Biot (1941, 1955) three dimensional linear model describes short time pore-fluid segregation, in case of consolidation, only for special geometries of the porous material, a slab (Mandel, 1953) or a sphere (Cryer, 1963). The existence of a fluid richer stationary state of the porous medium can not be proven in the context of linear Biot theory endowed with Darcy solid-fluid viscous coupling. Recently, in-situ observations and experimental studies have pointed out the formation of compaction bands in rocks and soils (Mollema and Antonellini, 1996). This phenomenon is typically connected with the occurrence of pore-fluid segregation (Holcomb and Olsson, 2003; Holcomb et al., 2007) in consolidation processes and, eventually, soil fluidization (Kolymbas, 1998): under the effect of an applied external pressure, or because of gravity, solid internal remodeling can induce the formation of non-connected fluid-filled cavities. Thus increasing the external loading causes fluid to remain trapped and therefore the fluid mass in the trapping chambers to increase with respect to that of the fluid flowing out of the solid matrix. The pore-fluid pressure is now
describe the occurrence of two states of equilibrium and the switching from one to the other by considering a kind of *phase transition*. Extending classical Biot's theory a more than quadratic strain energy potential is postulated, depending on the strain of the porous material and the variation of the fluid mass density (measured with respect to the skeleton reference volume). When the consolidating pressure is strong enough the existence of two distinct minima is proven.

**Keywords:** Granular materials, Porous media, Phase transition, Bifurcation theory

# 1 Introduction

Porous media confined into a fluid infinite reservoir, or suffering consolidation loading, can exhibit solid-solid and solid-fluid phase transitions. These last can be observed because of different phenomena, as gravity driven solid-fluid separation of suspensions (see Bürger et al., 2000) or solid-fluid segregation of soils (see e.g. Nichols et al., 1994; Vardoulakis, 2004a,b). It will be the purpose of this paper to investigate this last transition.

The Biot (1941, 1955) three dimensional linear model describes short time pore-fluid segregation, in case of consolidation, only for special geometries of the porous material, a slab (Mandel, 1953) or a sphere (Cryer, 1963). The existence of a fluid richer stationary state of the porous medium can not be proven in the context of linear Biot theory endowed with Darcy solid-fluid viscous coupling. Recently, in-situ observations and experimental studies have pointed out the formation of compaction bands in rocks and soils (Mollema and Antonellini, 1996). This phenomenon is typically connected with the occurrence of pore-fluid segregation (Holcomb and Olsson, 2003; Holcomb et al., 2007) in consolidation processes and, eventually, soil fluidization (Kolymbas, 1998): under the effect of an applied external pressure, or because of gravity, solid internal remodeling can induce the formation of non-connected fluid-filled cavities. Thus increasing the external loading causes fluid to remain trapped and therefore the fluid mass in the trapping chambers to increase with respect to that of the fluid flowing out of the solid matrix. The pore-fluid pressure is now
capable to induce unbalance of forces acting on the soil grains, so causing fluidization. Nichols

<figure><img src="image_2.png" /><figcaption>Figure 1: Qualitative picture of the behavior of the fluid mass, $m_f$, with respect to the consolidating pressure $p$.</figcaption></figure>

et al. (1994) have experimentally demonstrated the existence of two phases in granular materials saturated by a fluid. A cylindrical perspex test vessel is filled by a granular test sample and water is injected from the bottom through the granular layer. Tuning water pressure increases the flow velocity and therefore the upward drag action on the grains. At low flow velocity, when the drag force is smaller than gravity, the standard phase is observed: the fluid flows through the solid which remains undeformed. When the velocity is increased another phase shows up: the drag force balances gravity and fluidization of the grains occurs.

The aim of this paper is to model solid-fluid phase transition describing in particular the showing up of a pore-fluid segregating state, when a consolidating external pressure is applied to the porous material. For this purpose we use a simple one-dimensional model, generalizing the Biot theory, where the phase transition is achieved modifying the standard Biot internal energy functional. A suitable potential energy including the effects of the external pressure $p$ will be considered.

Also in this model the free energy depends on two fields, the deformation $\varepsilon$ of the porous matrix and the density of the fluid $m_f$, measured with respect to the solid reference volume. We shall describe the following phase transition: there exists $p_c > 0$ such that for $0 \le p \le p_c$, that is low
capable to induce unbalance of forces acting on the soil grains, so causing fluidization. Nichols

<figure><img src="image_2.png" /><figcaption>Figure 1: Qualitative picture of the behavior of the fluid mass, $m_f$, with respect to the consolidating pressure $p$.</figcaption></figure>

et al. (1994) have experimentally demonstrated the existence of two phases in granular materials saturated by a fluid. A cylindrical perspex test vessel is filled by a granular test sample and water is injected from the bottom through the granular layer. Tuning water pressure increases the flow velocity and therefore the upward drag action on the grains. At low flow velocity, when the drag force is smaller than gravity, the standard phase is observed: the fluid flows through the solid which remains undeformed. When the velocity is increased another phase shows up: the drag force balances gravity and fluidization of the grains occurs.

The aim of this paper is to model solid-fluid phase transition describing in particular the showing up of a pore-fluid segregating state, when a consolidating external pressure is applied to the porous material. For this purpose we use a simple one-dimensional model, generalizing the Biot theory, where the phase transition is achieved modifying the standard Biot internal energy functional. A suitable potential energy including the effects of the external pressure $p$ will be considered.

Also in this model the free energy depends on two fields, the deformation $\varepsilon$ of the porous matrix and the density of the fluid $m_f$, measured with respect to the solid reference volume. We shall describe the following phase transition: there exists $p_c > 0$ such that for $0 \le p \le p_c$, that is low
pressure, there exists a single stationary state with fluid mass $m_{f1}(p)$ and solid deformation $\varepsilon_1(p)$; while for $p > p_c$, that is high pressure, a second phase ($m_{f2}(p)$, $\varepsilon_2(p)$) appears. Fluid density $m_{f2}$ is greater than $m_{f1}$ and increases with the pressure $p$, see Fig 1.

The first solution corresponds to the case when the fluid is not confined inside the matrix and can flow freely back and forth from the inside to the outside of the solid. At equilibrium the density of the internal fluid equals that of the external infinite reservoir or, in other words the pore of the solid matrix are connected. The second solution, on the other hand, corresponds to the case when the mass density of the fluid is not that of the fluid in the external reservoir. This is due to the fact that the porous material starts to behave as a closed system rather than an open system, as the pore connecting ducts become thinner and thinner.

This point can be supported by means of a thermodynamic argument: the flow of the liquid back and forth from the inside to the outside of the porous medium is a thermodynamic transformation at constant temperature, pressure, and volume, typical of open systems. The equilibrium is achieved when the internal Gibbs free energy $G$ equals the external one; if the infinitesimal mass $dm$ of fluid exits the matrix, the mass outside will vary of the amount $-dm$. Since the infinitesimal variation of the internal and the external Gibbs free energy are given by $dG_i = \mu_i dm_i = \mu_i dm$ and $dG_e = \mu_e dm_e = -\mu_e dm$ respectively, where $\mu$ is the chemical potential of the fluid, we have that at equilibrium $dG_i = dG_e$ implies $\mu_i = \mu_e$. Recalling that $\mu = \partial U/\partial m$, with $U$ the internal energy of the fluid, for reasonable choices of the function $U$, we have that the equality of the internal and external chemical potential reflects into the equality of the internal and external fluid density.

It is quite natural that, supposed to limit our discussion to small values of the external pressure, the deformation of the solid is grossly proportional to the external pressure $p$, namely, $\varepsilon(p) \approx -Bp$ for some positive constant $B$ (recall that for a compressed solid matrix the deformation is negative). In the unique low pressure phase, since as a consequence of the pressure $p$ some liquid will exit the solid, we suppose that $m_{f1}(p)$ decreases proportionally to $p$. Concerning the second phase $m_{f2}(p)$, the only constraint will be $m_{f2}(p) > m_{f1}(p)$. Indeed in this phase we guess the solid structure is modified and room for some liquid is made.
pressure, there exists a single stationary state with fluid mass $m_{f1}(p)$ and solid deformation $\varepsilon_1(p)$; while for $p > p_c$, that is high pressure, a second phase ($m_{f2}(p)$, $\varepsilon_2(p)$) appears. Fluid density $m_{f2}$ is greater than $m_{f1}$ and increases with the pressure $p$, see Fig 1.

The first solution corresponds to the case when the fluid is not confined inside the matrix and can flow freely back and forth from the inside to the outside of the solid. At equilibrium the density of the internal fluid equals that of the external infinite reservoir or, in other words the pore of the solid matrix are connected. The second solution, on the other hand, corresponds to the case when the mass density of the fluid is not that of the fluid in the external reservoir. This is due to the fact that the porous material starts to behave as a closed system rather than an open system, as the pore connecting ducts become thinner and thinner.

This point can be supported by means of a thermodynamic argument: the flow of the liquid back and forth from the inside to the outside of the porous medium is a thermodynamic transformation at constant temperature, pressure, and volume, typical of open systems. The equilibrium is achieved when the internal Gibbs free energy $G$ equals the external one; if the infinitesimal mass $dm$ of fluid exits the matrix, the mass outside will vary of the amount $-dm$. Since the infinitesimal variation of the internal and the external Gibbs free energy are given by $dG_i = \mu_i dm_i = \mu_i dm$ and $dG_e = \mu_e dm_e = -\mu_e dm$ respectively, where $\mu$ is the chemical potential of the fluid, we have that at equilibrium $dG_i = dG_e$ implies $\mu_i = \mu_e$. Recalling that $\mu = \partial U/\partial m$, with $U$ the internal energy of the fluid, for reasonable choices of the function $U$, we have that the equality of the internal and external chemical potential reflects into the equality of the internal and external fluid density.

It is quite natural that, supposed to limit our discussion to small values of the external pressure, the deformation of the solid is grossly proportional to the external pressure $p$, namely, $\varepsilon(p) \approx -Bp$ for some positive constant $B$ (recall that for a compressed solid matrix the deformation is negative). In the unique low pressure phase, since as a consequence of the pressure $p$ some liquid will exit the solid, we suppose that $m_{f1}(p)$ decreases proportionally to $p$. Concerning the second phase $m_{f2}(p)$, the only constraint will be $m_{f2}(p) > m_{f1}(p)$. Indeed in this phase we guess the solid structure is modified and room for some liquid is made.
In spite of the appealing simplicity of these arguments, it turns out that the actual situation is complex to be studied. The first goal is to understand basic physical phenomena and to determine how phase transition in question is affected by the external pressure.

The present paper is organized as follows. In §2 we describe the generalized Biot model and introduce the free energy functional. In §3, via an analytical minimization of the functional, we study the stationary points and their character. Finally in §4 we discuss our results.

## 2 The model

*Kinematics.* Let $\mathcal{B}_S, \mathcal{B}_f$ be the reference configurations of the solid and fluid components; and $\mathcal{E}$ the Euclidean space of positions. To specify the current configuration of the system, that is the configuration at time $t \in \mathbb{R}$, two families of diffeomorphisms, $\{\chi_{S,t} : \mathcal{B}_S \to \mathcal{E}, t \in \mathbb{R}\}$ and $\{\phi_{f,t} : \mathcal{B}_S \to \mathcal{B}_f, t \in \mathbb{R}\}$ are introduced.

The current solid configuration is given by the solid placement map $\chi_{S,t}$; for any $X_S \in \mathcal{B}_S$, $x = \chi_{S,t}(X_S)$ is the position occupied, at time $t$ in the Euclidean space $\mathcal{E}$, by the solid material particle $X_S$. The map $\phi_{f,t}$, on the other hand, identifies the fluid material particle $X_f$ in $\mathcal{B}_f$ which, at time $t$, occupies the same current place $x$ as the solid particle $X_S$. This description of the kinematics of the fluid is completely consistent with the Eulerean point of view adopted in standard fluid mechanics: the focus is not on the placement of the fluid particles, but on the particle, which at time $t$, occupies the current place $x$. As a consequence the reference configuration of the solid $\mathcal{B}_S$, which in the following will be the reference configuration of the system, is assumed to be a known subdomain of $\mathcal{E}$, while the one of the fluid, $\mathcal{B}_f$, is unknown, to be determined by the map $\phi_{f,t}$. As usual also the current configuration of the solid is unknown. Bearing in mind the definition of the map $\phi_{f,t}$ we shall call $\mathcal{B} := \chi_{S,t}(\mathcal{B}_S)$ be the current configuration of the system. The so called fluid placement map $\chi_{f,t} : \mathcal{B}_f \to \mathcal{E}$ can be constructed starting from $\chi_{S,t}$ and $\phi_{f,t}$. Indeed once we set $\chi_{f,t} := \chi_{S,t} \circ \phi_{f,t}^{-1}$, for any $X_f \in \mathcal{B}_f$, $\chi_{f,t}(X_f)$ represents the position occupied in the current configuration by the fluid particle $X_f$. For more details we refer to Sciarra et al. (2008).
In spite of the appealing simplicity of these arguments, it turns out that the actual situation is
complex to be studied. The first goal is to understand basic physical phenomena and to determine
how phase transition in question is affected by the external pressure.

The present paper is organized as follows. In §2 we describe the generalized Biot model and
introduce the free energy functional. In §3, via an analytical minimization of the functional, we
study the stationary points and their character. Finally in §4 we discuss our results.

## 2 The model

*Kinematics.* Let $\mathcal{B}_S, \mathcal{B}_f$ be the reference configurations of the solid and fluid components; and $\mathcal{E}$ the Euclidean space of positions. To specify the current configuration of the system, that is the configuration at time $t \in \mathbb{R}$, two families of diffeomorphisms, $\{\chi_{S,t} : \mathcal{B}_S \to \mathcal{E}, t \in \mathbb{R}\}$ and $\{\phi_{f,t} : \mathcal{B}_S \to \mathcal{B}_f, t \in \mathbb{R}\}$ are introduced.

The current solid configuration is given by the solid placement map $\chi_{S,t}$; for any $X_S \in \mathcal{B}_S$, $x = \chi_{S,t}(X_S)$ is the position occupied, at time $t$ in the Euclidean space $\mathcal{E}$, by the solid material particle $X_S$. The map $\phi_{f,t}$, on the other hand, identifies the fluid material particle $X_f$ in $\mathcal{B}_f$ which, at time $t$, occupies the same current place $x$ as the solid particle $X_S$. This description of the kinematics of the fluid is completely consistent with the Eulerean point of view adopted in standard fluid mechanics: the focus is not on the placement of the fluid particles, but on the particle, which at time $t$, occupies the current place $x$. As a consequence the reference configuration of the solid $\mathcal{B}_S$, which in the following will be the *the reference configuration of the system*, is assumed to be a known subdomain of $\mathcal{E}$, while the one of the fluid, $\mathcal{B}_f$, is unknown, to be determined by the map $\phi_{f,t}$. As usual also the current configuration of the solid is unknown. Bearing in mind the definition of the map $\phi_{f,t}$ we shall call $\mathcal{B} := \chi_{S,t}(\mathcal{B}_S)$ be the *current configuration of the system*. The so called *fluid placement map* $\chi_{f,t} : \mathcal{B}_f \to \mathcal{E}$ can be constructed starting from $\chi_{S,t}$ and $\phi_{f,t}$. Indeed once we set $\chi_{f,t} := \chi_{S,t} \circ \phi_{f,t}^{-1}$, for any $X_f \in \mathcal{B}_f$, $\chi_{f,t}(X_f)$ represents the position occupied in the current configuration by the fluid particle $X_f$. For more details we refer to Sciarra et al. (2008).
**Strain.** Let $F_{S,t} := \nabla \chi_{S,t}$ and $\Phi_{f,t} := \nabla \phi_{f,t}$ be the gradients of the maps $\chi_{S,t}$ and $\phi_{f,t}$, respectively (to clarify notations we remark that $\nabla$ indicates, in this context, spatial derivative independently of the domain of the map on which it operates). $F_{S,t}$ is typically named deformation of the solid. Since derivatives are taken with respect to $X_S \in \mathcal{B}_S$, those gradients are usually called Lagrangean gradients. We also define the gradient $F_{f,t} := \nabla \chi_{f,t}$; the chain rule yields $F_{f,t}(X_f) = \nabla \chi_{f,t}(X_f) = \nabla(\chi_{S,t}(\phi_{f,t}^{-1}(X_f))) = F_{S,t}(\phi_{f,t}^{-1}(X_f))\nabla\phi_{f,t}^{-1}(X_f) = F_{S,t}(X_S)\Phi_{f,t}(X_S)^{-1}$, where by definition $X_f = \phi_{f,t}(X_S)$.

Let $J_{\alpha,t} := |F_{\alpha,t}|$, with $\alpha = s,f$, be the Jacobian of the transformation $\chi_{\alpha,t}$ measuring the ratio between current and reference volumes; we define the Green–Lagrange strain of the solid $\varepsilon := (F_{S,t}^\top F_{S,t} - I)/2$, where $I$ is the second order identity tensor.

From now on we shall restrict our attention just to one–dimensional space of positions, having in mind to consider applications of the present model to the so–called consolidation problem (see Biot, 1941; Terzaghi, 1946; Cryer, 1963). In this framework tensorial quantities restrict to scalars; in particular the Green–Lagrange strain reduces to $\varepsilon := (J_{S,t}^2 - 1)/2$.

**Mass balance.** Let $\varrho_{0,\alpha} : \mathcal{B}_\alpha \to \mathbb{R}$, with $\alpha = s,f$, be the solid and fluid reference densities. The total mass of each of the two components,

$$ M_{\alpha} := \int_{\mathcal{B}_{\alpha}} \varrho_{0,\alpha}(X_{\alpha}) \, dX_{\alpha} \qquad (1) $$

with $\alpha = s,f$, is supposed to be constant with respect to time $t$, which means that the mass is conserved when passing from the reference to the current configuration. Let $\varrho_{\alpha,t}$, with $\alpha = s,f$, the solid and fluid current densities, mass conservation reads

$$ \int_{\mathcal{B}_{\alpha}} \varrho_{0,\alpha}(X_{\alpha}) \, dX_{\alpha} = \int_{\mathcal{B}} \varrho_{\alpha,t}(x) \, dx = \int_{\mathcal{B}_{\alpha}} \varrho_{\alpha,t}(\chi_{\alpha,t}(X_{\alpha})) J_{\alpha,t}(X_{\alpha}) \, dX_{\alpha} \qquad (2) $$

which in the local form, i.e. $\forall X_\alpha \in \mathcal{B}_\alpha$, becomes $\varrho_{\alpha,t}(\chi_{\alpha,t}(X_\alpha))J_{\alpha,t}(X_\alpha) = \varrho_{0,\alpha}(X_\alpha)$. Note that the initial value $\varrho_{\alpha,0}$ of the current density of the $\alpha$-th constituent can be assumed equal to the reference density $\varrho_{0,\alpha}$. Using the map $\phi_{f,t}$ allows for introducing the solid Lagrangean mass density
**Strain.** Let $F_{S,t} := \nabla \chi_{S,t}$ and $\Phi_{f,t} := \nabla \phi_{f,t}$ be the gradients of the maps $\chi_{S,t}$ and $\phi_{f,t}$, respectively (to clarify notations we remark that $\nabla$ indicates, in this context, spatial derivative independently of the domain of the map on which it operates). $F_{S,t}$ is typically named deformation of the solid. Since derivatives are taken with respect to $X_S \in \mathcal{B}_S$, those gradients are usually called Lagrangean gradients. We also define the gradient $F_{f,t} := \nabla \chi_{f,t}$; the chain rule yields $F_{f,t}(X_f) = \nabla \chi_{f,t}(X_f) = \nabla(\chi_{S,t}(\phi_{f,t}^{-1}(X_f))) = F_{S,t}(\phi_{f,t}^{-1}(X_f))\nabla\phi_{f,t}^{-1}(X_f) = F_{S,t}(X_S)\Phi_{f,t}(X_S)^{-1}$, where by definition $X_f = \phi_{f,t}(X_S)$.

Let $J_{\alpha,t} := |F_{\alpha,t}|$, with $\alpha = s,f$, be the Jacobian of the transformation $\chi_{\alpha,t}$ measuring the ratio between current and reference volumes; we define the Green–Lagrange strain of the solid $\varepsilon := (F_{S,t}^\top F_{S,t} - I)/2$, where $I$ is the second order identity tensor.

From now on we shall restrict our attention just to one–dimensional space of positions, having in mind to consider applications of the present model to the so–called consolidation problem (see Biot, 1941; Terzaghi, 1946; Cryer, 1963). In this framework tensorial quantities restrict to scalars; in particular the Green–Lagrange strain reduces to $\varepsilon := (J_{S,t}^2 - 1)/2$.

**Mass balance.** Let $\varrho_{0,\alpha} : \mathcal{B}_\alpha \to \mathbb{R}$, with $\alpha = s,f$, be the solid and fluid reference densities. The total mass of each of the two components,

$$M_{\alpha} := \int_{\mathcal{B}_{\alpha}} \varrho_{0,\alpha}(X_{\alpha}) \, dX_{\alpha} \qquad (1)$$

with $\alpha = s,f$, is supposed to be constant with respect to time $t$, which means that the mass is conserved when passing from the reference to the current configuration. Let $\varrho_{\alpha,t}$, with $\alpha = s,f$, the solid and fluid current densities, mass conservation reads

$$\int_{\mathcal{B}_{\alpha}} \varrho_{0,\alpha}(X_{\alpha}) \, dX_{\alpha} = \int_{\mathcal{B}} \varrho_{\alpha,t}(x) \, dx = \int_{\mathcal{B}_{\alpha}} \varrho_{\alpha,t}(\chi_{\alpha,t}(X_{\alpha})) J_{\alpha,t}(X_{\alpha}) \, dX_{\alpha} \qquad (2)$$

which in the local form, i.e. $\forall X_\alpha \in \mathcal{B}_\alpha$, becomes $\varrho_{\alpha,t}(\chi_{\alpha,t}(X_\alpha))J_{\alpha,t}(X_\alpha) = \varrho_{0,\alpha}(X_\alpha)$. Note that the initial value $\varrho_{\alpha,0}$ of the current density of the $\alpha$-th constituent can be assumed equal to the reference density $\varrho_{0,\alpha}$. Using the map $\phi_{f,t}$ allows for introducing the solid Lagrangean mass density
of the fluid constituent:

$$m_{f,t}(X_s) := \varrho_{0,f}(\phi_{f,t}(X_s))\det\Phi_{f,t}(X_s) \quad (3)$$

The admissible deformations of the porous continuum are therefore completely known once the Green–Lagrange strain and the solid Lagrangean mass density of the fluid are determined.

*Overall potential.* To study the equilibrium property of the system at constant temperature a suitable overall potential energy $\Phi$, per unit volume, given by the sum of the *Helmoltz free energy* $\Psi$ and the potential of external forces, can be introduced. Since our goal is that of modeling soil consolidation, external loading will only consist of a pure pressure acting on the solid skeleton, which implies the overall potential to be $\Phi = \Psi + pJ_s$.

As already noticed in a one–dimensional space of positions $\varepsilon = (J_s^2 - 1)/2$, if we restrict the discussion to the regime of small deformations, namely $J_s \approx 1$ and $\varepsilon \approx 0$, we can expand around $J_s = 1$ and get $\varepsilon \approx J_s$ (*geometrical linearization*).

For the sake of simplicity, from now on we shall denote with $m$ the increment of fluid mass density $m_f$, with respect to a suitable reference value $m_{0,f}$. For small deformations a reasonable expression for the dimensionless potential density of isotropic porous materials is given, in the framework of the Biot (1941, 1955) theory, by the following quadratic form

$$\Phi^B(m, \varepsilon) = p\varepsilon + \frac{1}{2}\varepsilon^2 + \frac{1}{2}a(m - b\varepsilon)^2 \quad (4)$$

where $a > 0$ is the ratio between the Biot modulus $M$ (see Coussy, 2004), and the solid bulk modulus, while $b > 0$ is the so–called Biot coefficient and measures the coupling between the solid and the fluid components; $p$ is made dimensionless with respect to the solid bulk modulus. It is immediate to show that the only stationary state of (4) is $m_B(p) = -bp$ and $\varepsilon_B(p) = -p$. Our hope is to describe the phase transition, driven by the pressure $p$, by considering additional third and fourth order terms in the overall potential. In this perspective we introduce a fourth order potential which reduces to the Biot one for small values of the parameters $m$ and $\varepsilon$, namely, in our
of the fluid constituent:

$$m_{f,t}(X_s) := \varrho_{0,f}(\phi_{f,t}(X_s))\det\Phi_{f,t}(X_s) \quad (3)$$

The admissible deformations of the porous continuum are therefore completely known once the Green–Lagrange strain and the solid Lagrangean mass density of the fluid are determined.

*Overall potential.* To study the equilibrium property of the system at constant temperature a suitable overall potential energy $\Phi$, per unit volume, given by the sum of the *Helmoltz free energy* $\Psi$ and the potential of external forces, can be introduced. Since our goal is that of modeling soil consolidation, external loading will only consist of a pure pressure acting on the solid skeleton, which implies the overall potential to be $\Phi = \Psi + pJ_s$.

As already noticed in a one-dimensional space of positions $\varepsilon = (J_s^2 - 1)/2$, if we restrict the discussion to the regime of small deformations, namely $J_s \approx 1$ and $\varepsilon \approx 0$, we can expand around $J_s = 1$ and get $\varepsilon \approx J_s$ (geometrical linearization).

For the sake of simplicity, from now on we shall denote with $m$ the increment of fluid mass density $m_f$, with respect to a suitable reference value $m_{0,f}$. For small deformations a reasonable expression for the dimensionless potential density of isotropic porous materials is given, in the framework of the Biot (1941, 1955) theory, by the following quadratic form

$$\Phi^B(m, \varepsilon) = p\varepsilon + \frac{1}{2}\varepsilon^2 + \frac{1}{2}a(m - b\varepsilon)^2 \quad (4)$$

where $a > 0$ is the ratio between the Biot modulus $M$ (see Coussy, 2004), and the solid bulk modulus, while $b > 0$ is the so-called Biot coefficient and measures the coupling between the solid and the fluid components; $p$ is made dimensionless with respect to the solid bulk modulus. It is immediate to show that the only stationary state of (4) is $m_B(p) = -bp$ and $\varepsilon_B(p) = -p$. Our hope is to describe the phase transition, driven by the pressure $p$, by considering additional third and fourth order terms in the overall potential. In this perspective we introduce a fourth order potential which reduces to the Biot one for small values of the parameters $m$ and $\varepsilon$, namely, in our
expression, the second order terms will be precisely the Biot ones. We set

$$ \Phi(m, \varepsilon) = \frac{1}{12}\alpha m^2(3m^2 - 8b\varepsilon m + 6b^2\varepsilon^2) + \Phi^B(m, \varepsilon) \quad (5) $$

where $\alpha = \alpha(a, b)$ is a positive real function of the physical parameters $a$ and $b$. While the function $\alpha$ can be chosen freely to tune the results with reasonable physical behaviors, the coefficients of the second order trinomial have been chosen so that (5) admits the local minimum $(m_1(p), \varepsilon_1(p))$, with $m_1 = b\varepsilon_1$, describing the *trivial* phase of the system similar to the one obtained in the framework of the Biot theory.

# 3 Phase transition

We discuss the stationary states of the system which are identified with the minima of the two variables function $\Phi(\varepsilon, m)$. We shall show that a point similar to the Biot one always exists, while, depending on the pressure, more precisely for a sufficiently large pressure, a second stationary state shows up. This section is devoted to the mathematical discussion of the phenomenon, its physical interpretation is postponed to the Section 4.

## 3.1 Stationary points of the overall potential

In order to find the stationary points of the overall potential we compute the first order partial derivatives of the potential energy $\Phi$

$$ \Phi_m = (m - b\varepsilon)(\alpha m^2 - \alpha b\varepsilon m + a) \quad (6) $$

and

$$ \Phi_{\varepsilon} = \frac{1}{3}\alpha b m^2 (3b\varepsilon - 2m) - ab(m - b\varepsilon) + \varepsilon + p \quad (7) $$
expression, the second order terms will be precisely the Biot ones. We set

$$ \Phi(m, \varepsilon) = \frac{1}{12} \alpha m^2 (3m^2 - 8b\varepsilon m + 6b^2\varepsilon^2) + \Phi^B(m, \varepsilon) \quad (5) $$

where $\alpha = \alpha(a, b)$ is a positive real function of the physical parameters $a$ and $b$. While the function $\alpha$ can be chosen freely to tune the results with reasonable physical behaviors, the coefficients of the second order trinomial have been chosen so that (5) admits the local minimum $(m_1(p), \varepsilon_1(p))$, with $m_1 = b\varepsilon_1$, describing the *trivial* phase of the system similar to the one obtained in the framework of the Biot theory.

# 3 Phase transition

We discuss the stationary states of the system which are identified with the minima of the two variables function $\Phi(\varepsilon, m)$. We shall show that a point similar to the Biot one always exists, while, depending on the pressure, more precisely for a sufficiently large pressure, a second stationary state shows up. This section is devoted to the mathematical discussion of the phenomenon, its physical interpretation is postponed to the Section 4.

## 3.1 Stationary points of the overall potential

In order to find the stationary points of the overall potential we compute the first order partial derivatives of the potential energy $\Phi$

$$ \Phi_m = (m - b\varepsilon)(\alpha m^2 - \alpha b\varepsilon m + a) \quad (6) $$

and

$$ \Phi_{\varepsilon} = \frac{1}{3}\alpha b m^2 (3b\varepsilon - 2m) - ab(m - b\varepsilon) + \varepsilon + p \quad (7) $$
Letting $\Phi_m = 0$ we get the following three solutions: $m_1 = b\epsilon$ and

$$m_{\pm} = \frac{b}{2} \left[ \epsilon \pm \sqrt{\epsilon^2 - \frac{4a}{\alpha b^2}} \right] \qquad (8)$$

To get the corresponding values of $\epsilon$, the equation $\Phi_{\epsilon} = 0$ must be solved with $m = m_1, m_{\pm}$. We first note that by inserting $m = m_1$ in $\Phi_{\epsilon} = 0$, we get

$$p = -\epsilon - \frac{1}{3}\alpha b^4 \epsilon^3 =: f_1(\epsilon) \qquad (9)$$

Studying the function $f_1$ it is immediate to deduce that (9) has a single real solution for any $p > 0$; we let $\epsilon_1$ such a unique solution and remark that the stationary point $(m_1, \epsilon_1)$ of the overall potential does exist for any choice of the parameter of the model.

To solve (7) with $m = m_{\pm}$, we bound the discussion to the physically relevant region $\epsilon < 0$. Note that, if $|\epsilon| \ge 2/(b\sqrt{\alpha/a})$ then

$$m_{\pm} = \frac{1}{2} b \epsilon \left[ 1 \mp \sqrt{1 - \frac{4a}{\alpha b^2 \epsilon^2}} \right] < 0 \qquad (10)$$

which implies, in particular, that $m_+ \ge m_- > m_1$; remark that for $|\epsilon|$ large enough, $m_+$ approaches 0 and $m_-$ approaches $m_1$. By letting $m = m_{\pm}$ in (7), we get

$$p = -\epsilon + ab[m_{\pm}(\epsilon) - b\epsilon] - \alpha b^2 \epsilon m_{\pm}^2(\epsilon) + \frac{2}{3}\alpha b m_{\pm}^3(\epsilon) =: f_{\pm}(\epsilon) \qquad (11)$$

Since $f_{\pm}$ are not defined around $\epsilon = 0$, namely, for $|\epsilon| < 2/(b\sqrt{\alpha/a})$, we can hope to describe some phase transition driven by $p$; indeed we guess that for $p$ small enough the equation (11) will have no real solution so that $(m_1, \epsilon_1)$ will be the sole stationary point of the overall potential. To explain our guess we draw, in Figure 2, $f_1, f_{\pm}$ as function of $\epsilon$ for the parameters specified in the captions; note that $f_-$ approaches $f_1$ for $|\epsilon|$ large, this is consistent with the behavior of $m_1$ and $m_-$. We let $p_c$ be the minimum of the function $f_+$ and $p'_c := f_+(\epsilon'_c) = f_-(\epsilon'_c)$ with $\epsilon'_c := -2/(b\sqrt{\alpha/a})$, the largest value of $\epsilon$ for which the additional stationary point of $\Phi$ appears. For $p < p_c$ the
Letting $\Phi_m = 0$ we get the following three solutions: $m_1 = b\epsilon$ and

$$m_{\pm} = \frac{b}{2} \left[ \epsilon \pm \sqrt{\epsilon^2 - \frac{4a}{\alpha b^2}} \right] \qquad (8)$$

To get the corresponding values of $\epsilon$, the equation $\Phi_{\epsilon} = 0$ must be solved with $m = m_1, m_{\pm}$. We first note that by inserting $m = m_1$ in $\Phi_{\epsilon} = 0$, we get

$$p = -\epsilon - \frac{1}{3}\alpha b^4 \epsilon^3 =: f_1(\epsilon) \qquad (9)$$

Studying the function $f_1$ it is immediate to deduce that (9) has a single real solution for any $p > 0$; we let $\epsilon_1$ such a unique solution and remark that the stationary point $(m_1, \epsilon_1)$ of the overall potential does exist for any choice of the parameter of the model.

To solve (7) with $m = m_{\pm}$, we bound the discussion to the physically relevant region $\epsilon < 0$. Note that, if $|\epsilon| \ge 2/(b\sqrt{\alpha/a})$ then

$$m_{\pm} = \frac{1}{2} b \epsilon \left[ 1 \mp \sqrt{1 - \frac{4a}{\alpha b^2 \epsilon^2}} \right] < 0 \qquad (10)$$

which implies, in particular, that $m_+ \ge m_- > m_1$; remark that for $|\epsilon|$ large enough, $m_+$ approaches 0 and $m_-$ approaches $m_1$. By letting $m = m_{\pm}$ in (7), we get

$$p = -\epsilon + ab[m_{\pm}(\epsilon) - b\epsilon] - \alpha b^2 \epsilon m_{\pm}^2(\epsilon) + \frac{2}{3}\alpha b m_{\pm}^3(\epsilon) =: f_{\pm}(\epsilon) \qquad (11)$$

Since $f_{\pm}$ are not defined around $\epsilon = 0$, namely, for $|\epsilon| < 2/(b\sqrt{\alpha/a})$, we can hope to describe some phase transition driven by $p$; indeed we guess that for $p$ small enough the equation (11) will have no real solution so that $(m_1, \epsilon_1)$ will be the sole stationary point of the overall potential. To explain our guess we draw, in Figure 2, $f_1, f_{\pm}$ as function of $\epsilon$ for the parameters specified in the captions; note that $f_-$ approaches $f_1$ for $|\epsilon|$ large, this is consistent with the behavior of $m_1$ and $m_-$. We let $p_c$ be the minimum of the function $f_+$ and $p'_c := f_+(\epsilon'_c) = f_-(\epsilon'_c)$ with $\epsilon'_c := -2/(b\sqrt{\alpha/a})$, the largest value of $\epsilon$ for which the additional stationary point of $\Phi$ appears. For $p < p_c$ the
<figure><img src="image_3.png" /><figcaption>Figure 2: The graph of $f_{\pm}$ and $f_1$ is plotted. Note that the dashed line corresponds to the graph of $f_{-}$. The values of these functions at $p$ determine the states of equilibrium of the porous medium. The used values of the constitutive parameters are $a = 1/2$, $b = 1$, $\alpha = 100$.</figcaption></figure>

solutions $(m_{\pm}, \varepsilon_{\pm})$ are not real, hence the system has a single phase, the one essentially due to the Biot model. For $p_c < p < p'_c$, the solution $(m_+, \varepsilon_+)$ with the good $p$-behavior (the one with smallest $\varepsilon_+$) should be interpreted as the second phase, while the other should be a saddle point of the overall potential. For $p > p'_c$, the unique solution $(m_+, \varepsilon_+)$ should be interpreted as the second phase, while $(m_-, \varepsilon_-)$ should be a saddle point of the overall potential. To prove that this interpretation is correct one should study the second order derivatives of $\Phi$ (see §3.2).

To study the equation (11) we recall that the two functions $f_{\pm}$ are defined for $\varepsilon < \varepsilon'_c := -2/(b\sqrt{\alpha/a})$ and note that

$$p'_c := f_{\pm}(\varepsilon'_c) = \frac{1}{b\sqrt{\alpha/a}}\left(2 + \frac{7}{3}b^2a\right) > 0 \quad (12)$$

To study the $\varepsilon \to -\infty$ limit, we recall that $\sqrt{1-x} = 1 - x/2 - x^2/8 + O(x^3)$, and using (10) we get

$$m_{\pm}(\varepsilon) = \frac{1}{2}b\varepsilon\left[2\delta_{-,\pm} \pm \frac{2a}{\alpha b^2 \varepsilon^2} \pm \frac{2a^2}{\alpha^2 b^4 \varepsilon^4} + O(1/\varepsilon^6)\right] \quad (13)$$

where $\delta_{-,\pm}$ is the Kronecker delta, such that $\delta_{-,+} = 1$ and $\delta_{+,+} = 0$. Substituting those expansions
<figure><img src="image_3.png" /><figcaption>Figure 2: The graph of $f_{\pm}$ and $f_1$ is plotted. Note that the dashed line corresponds to the graph of $f_{-}$. The values of these functions at $p$ determine the states of equilibrium of the porous medium. The used values of the constitutive parameters are $a = 1/2$, $b = 1$, $\alpha = 100$.</figcaption></figure>

solutions $(m_{\pm}, \varepsilon_{\pm})$ are not real, hence the system has a single phase, the one essentially due to the Biot model. For $p_c < p < p'_c$, the solution $(m_+, \varepsilon_+)$ with the good $p$-behavior (the one with smallest $\varepsilon_+$) should be interpreted as the second phase, while the other should be a saddle point of the overall potential. For $p > p'_c$, the unique solution $(m_+, \varepsilon_+)$ should be interpreted as the second phase, while $(m_-, \varepsilon_-)$ should be a saddle point of the overall potential. To prove that this interpretation is correct one should study the second order derivatives of $\Phi$ (see §3.2).

To study the equation (11) we recall that the two functions $f_{\pm}$ are defined for $\varepsilon < \varepsilon'_c := -2/(b\sqrt{\alpha/a})$ and note that

$$p'_c := f_{\pm}(\varepsilon'_c) = \frac{1}{b\sqrt{\alpha/a}}\left(2 + \frac{7}{3}b^2a\right) > 0 \quad (12)$$

To study the $\varepsilon \to -\infty$ limit, we recall that $\sqrt{1-x} = 1 - x/2 - x^2/8 + O(x^3)$, and using (10) we get

$$m_{\pm}(\varepsilon) = \frac{1}{2}b\varepsilon\left[2\delta_{-,\pm} \pm \frac{2a}{\alpha b^2 \varepsilon^2} \pm \frac{2a^2}{\alpha^2 b^4 \varepsilon^4} + O(1/\varepsilon^6)\right] \quad (13)$$

where $\delta_{-,\pm}$ is the Kronecker delta, such that $\delta_{-,+} = 1$ and $\delta_{+,+} = 0$. Substituting those expansions
in the definition (11) of $f_{\pm}$ and keeping track of the not vanishing terms we get

$$f_{+}(\varepsilon) \sim -\varepsilon(1 + ab^2) \quad \text{and} \quad f_{-}(\varepsilon) \sim -\varepsilon - \frac{1}{3}\alpha b^4 \varepsilon^3 \qquad (14)$$

Note that for $\varepsilon \to -\infty$ the function $f_{-}$ behaves precisely as $f_1$.

We now compute the first order derivative. Recalling the definition (11) of $f_{\pm}$, see equation (6), and the fact that $m_{\pm}$ are solutions of the equation $\Phi_m = 0$, we get

$$ \begin{aligned} f'_{\pm}(\varepsilon) &= -1 - ab^2 - \alpha b^2 m_{\pm}^2(\varepsilon) + b m'_{\pm}(\varepsilon) [a - 2\alpha b \varepsilon m_{\pm}(\varepsilon) + 2\alpha m_{\pm}^2(\varepsilon)] \\ &= -1 - ab^2 - \alpha b^2 m_{\pm}^2(\varepsilon) - ab m'_{\pm}(\varepsilon) \end{aligned} \qquad (15) $$

Note the the first three terms are clearly negative, so that the sign of the derivatives depends essentially on the sign of the fourth term. We then recall (10), which is valid for $\varepsilon < 0$, and compute

$$m'_{\pm}(\varepsilon) = \frac{1}{\varepsilon}m_{\pm}(\varepsilon) \mp \frac{2a}{\alpha b \varepsilon^2} \left[1 - \frac{4a}{\alpha b^2 \varepsilon^2}\right]^{-1/2} \qquad (16)$$

By using (16) and (15) it follows immediately that $f'_{\pm}(\varepsilon) \to \pm\infty$ for $\varepsilon \to \varepsilon'_{c}$ from the left.

Notice that $m'_{-}(\varepsilon)$, for $\varepsilon < 0$, is positive. Hence, by (15) it follows that $f'_{-}(\varepsilon) < 0$ for $\varepsilon < \varepsilon'_{c}$. The function $f_{-}$ decreases monotonously from $+\infty$ to $p'_{c}$ as $\varepsilon$ goes from $-\infty$ to $\varepsilon'_{c}$. We then conclude that the stationary point $(m_{-}, \varepsilon_{-})$ of the overall potential exists and is unique for $p \ge p'_{c}$.

The study of $f_{+}$ is slightly more difficult. First of all we recall that $0 > m_{+}(\varepsilon) > b\varepsilon$, see (10), hence, using the definition (11), we get

$$f_{+}(\varepsilon) = -\varepsilon + ab(m_{+}(\varepsilon) - b\varepsilon) + \alpha b m_{+}^{2}(\varepsilon) \left[\frac{2}{3}m_{+}(\varepsilon) - b\varepsilon\right] > 0 \qquad (17)$$

Thus $f_{+}(\varepsilon)$ is a positive function tending to $+\infty$ for $\varepsilon \to -\infty$ and approaching $f_{+}(\varepsilon'_{c}) > 0$, see (12), with positive slope. Hence there must be at least a minimum of the function $f_{+}(\varepsilon)$ in the region $\varepsilon \in (-\infty, \varepsilon'_{c}]$. Moreover, it is possible easy to see that the equation $f'_{+}(\varepsilon) = 0$ is biquadratic in $\varepsilon$. Hence it can have either zero or one or two negative solutions, the only possible case, compatible
in the definition (11) of $f_{\pm}$ and keeping track of the not vanishing terms we get

$$f_{+}(\varepsilon) \sim -\varepsilon(1 + ab^2) \quad \text{and} \quad f_{-}(\varepsilon) \sim -\varepsilon - \frac{1}{3}\alpha b^4 \varepsilon^3 \qquad (14)$$

Note that for $\varepsilon \to -\infty$ the function $f_{-}$ behaves precisely as $f_1$.

We now compute the first order derivative. Recalling the definition (11) of $f_{\pm}$, see equation (6), and the fact that $m_{\pm}$ are solutions of the equation $\Phi_m = 0$, we get

$$\begin{aligned} f'_{\pm}(\varepsilon) &= -1 - ab^2 - \alpha b^2 m_{\pm}^2(\varepsilon) + b m'_{\pm}(\varepsilon) [a - 2\alpha b \varepsilon m_{\pm}(\varepsilon) + 2\alpha m_{\pm}^2(\varepsilon)] \\ &= -1 - ab^2 - \alpha b^2 m_{\pm}^2(\varepsilon) - ab m'_{\pm}(\varepsilon) \end{aligned} \qquad (15)$$

Note the the first three terms are clearly negative, so that the sign of the derivatives depends essentially on the sign of the fourth term. We then recall (10), which is valid for $\varepsilon < 0$, and compute

$$m'_{\pm}(\varepsilon) = \frac{1}{\varepsilon}m_{\pm}(\varepsilon) \mp \frac{2a}{\alpha b \varepsilon^2} \left[1 - \frac{4a}{\alpha b^2 \varepsilon^2}\right]^{-1/2} \qquad (16)$$

By using (16) and (15) it follows immediately that $f'_{\pm}(\varepsilon) \to \pm\infty$ for $\varepsilon \to \varepsilon'_{c}$ from the left.

Notice that $m'_{-}(\varepsilon)$, for $\varepsilon < 0$, is positive. Hence, by (15) it follows that $f'_{-}(\varepsilon) < 0$ for $\varepsilon < \varepsilon'_{c}$. The function $f_{-}$ decreases monotonously from $+\infty$ to $p'_{c}$ as $\varepsilon$ goes from $-\infty$ to $\varepsilon'_{c}$. We then conclude that the stationary point $(m_{-}, \varepsilon_{-})$ of the overall potential exists and is unique for $p \ge p'_{c}$.

The study of $f_{+}$ is slightly more difficult. First of all we recall that $0 > m_{+}(\varepsilon) > b\varepsilon$, see (10), hence, using the definition (11), we get

$$f_{+}(\varepsilon) = -\varepsilon + ab(m_{+}(\varepsilon) - b\varepsilon) + \alpha b m_{+}^{2}(\varepsilon) \left[\frac{2}{3}m_{+}(\varepsilon) - b\varepsilon\right] > 0 \qquad (17)$$

Thus $f_{+}(\varepsilon)$ is a positive function tending to $+\infty$ for $\varepsilon \to -\infty$ and approaching $f_{+}(\varepsilon'_{c}) > 0$, see (12), with positive slope. Hence there must be at least a minimum of the function $f_{+}(\varepsilon)$ in the region $\varepsilon \in (-\infty, \varepsilon'_{c}]$. Moreover, it is possible easy to see that the equation $f'_{+}(\varepsilon) = 0$ is biquadratic in $\varepsilon$. Hence it can have either zero or one or two negative solutions, the only possible case, compatible
with the above mentioned properties of $f_+$, is that the negative solution is unique.

Finally, we set $p_c := f_+(\epsilon_c) > 0$ and remark that for $p < p_c$ the equations $f_\pm(\epsilon) = p$ have no real solution. For $p_c < p < p'_c$, the equation $f_-(\epsilon) = p$ has no real solution, while $f_+(\epsilon) = p$ has two real solutions $(m_+^1, \epsilon_+^1)$ and $(m_+^2, \epsilon_+^2)$ with $\epsilon_+^1 < \epsilon_+^2$. For $p > p'_c$, both $f_-(\epsilon) = p$ and $f_+(\epsilon) = p$ have a unique real solution, respectively denoted by $(m_-, \epsilon_-)$ and $(m_+^1, \epsilon_+^1)$. According to the above depicted scenario, see the discussion below (11), we expect that $(m_+^1, \epsilon_+^1)$ is a minimum of the overall potential, while $(m_+^2, \epsilon_+^2)$ and $(m_-, \epsilon_-)$ are saddle points. The validity of this guess will be proven in the next section.

## 3.2 Character of the stationary points

We study, now, the character of the stationary points of the overall potential. Computing the second order derivatives of $\Phi$ with respect to $m$ and $\epsilon$ we get

$$ 
\begin{aligned}
\Phi_{mm} &= (\alpha m^2 - \alpha b \epsilon m + a) + (m - b \epsilon)(2\alpha m - \alpha b \epsilon) \\
\Phi_{m\epsilon} &= -ab + 2\alpha b m(b\epsilon - m) \\
\Phi_{\epsilon\epsilon} &= 1 + ab^2 + \alpha b^2 m^2 > 0
\end{aligned}
\qquad (18) 
$$

Using that $m_1 = b\epsilon_1$, we are able to compute the Hessian $H(m_1, \epsilon_1) = a(1 + \alpha m_1^2 b^2) > 0$ and, thus, conclude that $(m_1, \epsilon_1)$ is a local minimum of (5). Hence it represents a stationary state of the model.

We have to study, now, the properties of the stationary points $(m_\pm, \epsilon_\pm)$. In view of this we give a nice expression of the Hessian computed in $(m_\pm(\epsilon), \epsilon)$. By using (18) and recalling that $m_\pm(\epsilon)$ is obtained by solving the equation $\Phi_m = 0$, see (6), we get

$$ 
\begin{aligned}
\Phi_{mm}(m_\pm(\epsilon), \epsilon) &= \alpha(m_\pm(\epsilon) - b\epsilon)(2m_\pm(\epsilon) - b\epsilon) \\
\Phi_{m\epsilon}(m_\pm(\epsilon), \epsilon) &= ab \\
\Phi_{\epsilon\epsilon}(m_\pm(\epsilon), \epsilon) &= 1 + b^2(a + \alpha(m_\pm(\epsilon))^2) > 0
\end{aligned}
\qquad (19) 
$$
with the above mentioned properties of $f_+$, is that the negative solution is unique.

Finally, we set $p_c := f_+(\epsilon_c) > 0$ and remark that for $p < p_c$ the equations $f_\pm(\epsilon) = p$ have no real solution. For $p_c < p < p'_c$, the equation $f_-(\epsilon) = p$ has no real solution, while $f_+(\epsilon) = p$ has two real solutions $(m_+^1, \epsilon_+^1)$ and $(m_+^2, \epsilon_+^2)$ with $\epsilon_+^1 < \epsilon_+^2$. For $p > p'_c$, both $f_-(\epsilon) = p$ and $f_+(\epsilon) = p$ have a unique real solution, respectively denoted by $(m_-, \epsilon_-)$ and $(m_+^1, \epsilon_+^1)$. According to the above depicted scenario, see the discussion below (11), we expect that $(m_+^1, \epsilon_+^1)$ is a minimum of the overall potential, while $(m_+^2, \epsilon_+^2)$ and $(m_-, \epsilon_-)$ are saddle points. The validity of this guess will be proven in the next section.

## 3.2 Character of the stationary points

We study, now, the character of the stationary points of the overall potential. Computing the second order derivatives of $\Phi$ with respect to $m$ and $\epsilon$ we get

$$ 
\begin{aligned}
\Phi_{mm} &= (\alpha m^2 - \alpha b \epsilon m + a) + (m - b \epsilon)(2\alpha m - \alpha b \epsilon) \\
\Phi_{m\epsilon} &= -ab + 2\alpha b m (b\epsilon - m) \\
\Phi_{\epsilon\epsilon} &= 1 + ab^2 + \alpha b^2 m^2 > 0
\end{aligned}
\qquad (18) 
$$

Using that $m_1 = b\epsilon_1$, we are able to compute the Hessian $H(m_1, \epsilon_1) = a(1 + \alpha m_1^2 b^2) > 0$ and, thus, conclude that $(m_1, \epsilon_1)$ is a local minimum of (5). Hence it represents a stationary state of the model.

We have to study, now, the properties of the stationary points $(m_\pm, \epsilon_\pm)$. In view of this we give a nice expression of the Hessian computed in $(m_\pm(\epsilon), \epsilon)$. By using (18) and recalling that $m_\pm(\epsilon)$ is obtained by solving the equation $\Phi_m = 0$, see (6), we get

$$ 
\begin{aligned}
\Phi_{mm}(m_\pm(\epsilon), \epsilon) &= \alpha(m_\pm(\epsilon) - b\epsilon)(2m_\pm(\epsilon) - b\epsilon) \\
\Phi_{m\epsilon}(m_\pm(\epsilon), \epsilon) &= ab \\
\Phi_{\epsilon\epsilon}(m_\pm(\epsilon), \epsilon) &= 1 + b^2(a + \alpha(m_\pm(\epsilon))^2) > 0
\end{aligned}
\qquad (19) 
$$
so that

$$
\begin{align}
\mathcal{H}_{\pm}(\varepsilon) := \mathcal{H}(m_{\pm}(\varepsilon), \varepsilon) &= [1 + b^2(a + \alpha(m_{\pm}(\varepsilon))^2)][\alpha(m_{\pm}(\varepsilon) - b\varepsilon)(2m_{\pm}(\varepsilon) - b\varepsilon)] - a^2b^2 \nonumber \\
&= [1 + b^2\alpha b\varepsilon m_{\pm}(\varepsilon)][\alpha(m_{\pm}(\varepsilon) - b\varepsilon)(2m_{\pm}(\varepsilon) - b\varepsilon)] - a^2b^2 \tag{20}
\end{align}
$$

By using (8) we get

$$
\mathcal{H}_{\pm}(\varepsilon) = -2a - a^2b^2 + \alpha b^2 \varepsilon^2 / 2 \mp \alpha b^2 (1/2 + ab^2) \varepsilon \sqrt{\varepsilon^2 - 4a / (\alpha b^2)} \quad (21)
$$

Since we have focused our discussion on the case $\epsilon < 0$, we have that

$$
\mathcal{H}_{\pm}(\varepsilon) = -2a - a^2b^2 + \alpha b^2 \varepsilon^2 / 2 \pm \alpha b^2 (1/2 + ab^2) \varepsilon^2 \sqrt{1 - 4a / (\alpha b^2 \varepsilon^2)} \quad (22)
$$

We use the above expression (22) to prove that the stationary point ($m_-, \varepsilon_-$) is a saddle point.
First we note that if such a stationary point exists then it must be necessarily $\varepsilon_- \le \varepsilon'_c$, so that
$4a/(\alpha b^2 \varepsilon_-^2) \le 1$. Remarking that, for $0 \le x \le 1$, one has $\sqrt{1-x} > 1-x$, by (22) we get the
bound

$$
\begin{align}
\mathcal{H}_{-}(\varepsilon_{-}) &\leq -2a - a^2b^2 + \alpha b^2 \varepsilon_{-}^2 / 2 - \alpha b^2 (1/2 + ab^2) \varepsilon_{-}^2 [1 - 4a / (\alpha b^2 \varepsilon_{-}^2)] \tag{23} \\
&= ab^2 [3a - \alpha b^2 \varepsilon_{-}^2] \leq ab^2 (3a - 4a) = -a^2b^2 < 0 \nonumber
\end{align}
$$

Where the last bound follows from the fact that, since $\varepsilon_- \le \varepsilon'_c$, we have $\alpha b^2 \varepsilon_-^2 \ge 4a$.

In order to study the character of the stationary points ($m_+^i, \varepsilon_+^i$), with $i = 1, 2$, we study the function $\mathcal{H}_+(\varepsilon)$ for $\varepsilon < 0$. Recall (22) and note that $\mathcal{H}_+(\varepsilon'_c) = -a^2b^2 < 0$ and $\mathcal{H}_+(\varepsilon) \to +\infty$ for $\varepsilon \to -\infty$. By computing the first derivative of $\mathcal{H}_+(\varepsilon)$ with respect to $\varepsilon$ we get

$$
\mathcal{H}'_+(\varepsilon) = \alpha b^2 \varepsilon \left[ 1 + \left(\frac{1}{2} + ab^2\right) \left(1 - \frac{4a}{\alpha b^2 \varepsilon^2}\right)^{-1/2} \left(2 - \frac{4a}{\alpha b^2 \varepsilon^2}\right) \right] \quad (24)
$$

which is clearly negative for $\varepsilon < \varepsilon'_c$. Thus, we have that $\mathcal{H}_+(\varepsilon)$ decreases from $+\infty$ to $-a^2b^2$ as $\varepsilon$
goes from $-\infty$ to $\varepsilon'_c$. Denoted by $\bar{\varepsilon}$ the unique negative zero of $\mathcal{H}_+(\varepsilon)$, we prove that $\bar{\varepsilon} = \varepsilon_c$, that
so that

$$ \begin{aligned} \mathcal{H}_{\pm}(\varepsilon) := \mathcal{H}(m_{\pm}(\varepsilon), \varepsilon) &= [1 + b^2(a + \alpha(m_{\pm}(\varepsilon))^2)][\alpha(m_{\pm}(\varepsilon) - b\varepsilon)(2m_{\pm}(\varepsilon) - b\varepsilon)] - a^2b^2 \\ &= [1 + b^2\alpha b\varepsilon m_{\pm}(\varepsilon)][\alpha(m_{\pm}(\varepsilon) - b\varepsilon)(2m_{\pm}(\varepsilon) - b\varepsilon)] - a^2b^2 \end{aligned} \quad (20) $$

By using (8) we get

$$ \mathcal{H}_{\pm}(\varepsilon) = -2a - a^2b^2 + \alpha b^2 \varepsilon^2 / 2 \mp \alpha b^2 (1/2 + ab^2) \varepsilon \sqrt{\varepsilon^2 - 4a/(\alpha b^2)} \quad (21) $$

Since we have focused our discussion on the case $\varepsilon < 0$, we have that

$$ \mathcal{H}_{\pm}(\varepsilon) = -2a - a^2b^2 + \alpha b^2 \varepsilon^2 / 2 \pm \alpha b^2 (1/2 + ab^2) \varepsilon^2 \sqrt{1 - 4a/(\alpha b^2 \varepsilon^2)} \quad (22) $$

We use the above expression (22) to prove that the stationary point $(m_{-}, \varepsilon_{-})$ is a saddle point. First we note that if such a stationary point exists then it must be necessarily $\varepsilon_{-} \le \varepsilon'_{c}$, so that $4a/(\alpha b^2 \varepsilon_{-}^2) \le 1$. Remarking that, for $0 \le x \le 1$, one has $\sqrt{1-x} > 1-x$, by (22) we get the bound

$$ \begin{aligned} \mathcal{H}_{-}(\varepsilon_{-}) &\le -2a - a^2b^2 + \alpha b^2 \varepsilon_{-}^2 / 2 - \alpha b^2 (1/2 + ab^2) \varepsilon_{-}^2 [1 - 4a/(\alpha b^2 \varepsilon_{-}^2)] \\ &= ab^2 [3a - \alpha b^2 \varepsilon_{-}^2] \le ab^2 (3a - 4a) = -a^2b^2 < 0 \end{aligned} \quad (23) $$

Where the last bound follows from the fact that, since $\varepsilon_{-} \le \varepsilon'_{c}$, we have $\alpha b^2 \varepsilon_{-}^2 \ge 4a$.

In order to study the character of the stationary points $(m_{+}^{i}, \varepsilon_{+}^{i})$, with $i = 1, 2$, we study the function $\mathcal{H}_{+}(\varepsilon)$ for $\varepsilon < 0$. Recall (22) and note that $\mathcal{H}_{+}(\varepsilon'_{c}) = -a^2b^2 < 0$ and $\mathcal{H}_{+}(\varepsilon) \to +\infty$ for $\varepsilon \to -\infty$. By computing the first derivative of $\mathcal{H}_{+}(\varepsilon)$ with respect to $\varepsilon$ we get

$$ \mathcal{H}'_{+}(\varepsilon) = \alpha b^2 \varepsilon \left[ 1 + \left(\frac{1}{2} + ab^2\right) \left(1 - \frac{4a}{\alpha b^2 \varepsilon^2}\right)^{-1/2} \left(2 - \frac{4a}{\alpha b^2 \varepsilon^2}\right) \right] \quad (24) $$

which is clearly negative for $\varepsilon < \varepsilon'_{c}$. Thus, we have that $\mathcal{H}_{+}(\varepsilon)$ decreases from $+\infty$ to $-a^2b^2$ as $\varepsilon$ goes from $-\infty$ to $\varepsilon'_{c}$. Denoted by $\bar{\varepsilon}$ the unique negative zero of $\mathcal{H}_{+}(\varepsilon)$, we prove that $\bar{\varepsilon} = \varepsilon_{c}$, that
is $\bar{\epsilon}$ is the minimum of the function $f_+(\epsilon)$. Note, indeed, that by definition

$$ \mathcal{H}_+(\epsilon) = \Phi_{mm}(m_+(\epsilon), \epsilon)\Phi_{\epsilon\epsilon}(m_+(\epsilon), \epsilon) - \Phi_{m\epsilon}^2(m_+(\epsilon), \epsilon) \quad (25) $$

Recall that the function $m_+(\epsilon)$ is implicitly defined by the equation $\Phi_m(m_+(\epsilon), \epsilon) = 0$ and note that (11) is equivalent to $f_+(\epsilon) = -\Phi_\epsilon(m_+(\epsilon), \epsilon) + p$. By using the chain rule and the implicit function theorem we get that

$$ \begin{aligned} f'_+(\epsilon) &= -\Phi_{\epsilon\epsilon}(m_+(\epsilon), \epsilon) - \Phi_{\epsilon m}(m_+(\epsilon), \epsilon) m'_+(\epsilon) = -\Phi_{\epsilon\epsilon}(m_+(\epsilon), \epsilon) + \frac{\Phi_{\epsilon m}^2(m_+(\epsilon), \epsilon)}{\Phi_{mm}(m_+(\epsilon), \epsilon)} \\ &= -\frac{1}{\Phi_{mm}(m_+(\epsilon), \epsilon)} \mathcal{H}_+(\epsilon) \end{aligned} \quad (26) $$

where in the last step we have used (25).

The above equality ensures that $\bar{\epsilon} = \epsilon_c$. Finally, since, $\epsilon_+^1 < \epsilon_c < \epsilon_+^2$, we have that $(m_+^1, \epsilon_+^1)$ is a minimum of (5) while $(m_+^2, \epsilon_+^2)$ is a saddle.

# 4 Results

In this section a discussion on the behavior of the minima of the overall potential (5), when varying the external pressure, will be presented. A parametric analysis in dependence of the coefficients introduced in the Biot model, a and b, as well as of the additional coefficient $\alpha(a, b)$ multiplying the fourth order terms of (5) will be also developed.

According to the experimental results on fluidization of soils, cited in the introduction (see Nichols et al., 1994; Holcomb and Olsson, 2003; Holcomb et al., 2007; Vardoulakis, 2004a,b), the proposed model is capable to exhibit, in the presence of an applied consolidating external pressure, an additional stationary state in the $(\epsilon, m)$ plane beyond the one corresponding to classical consolidation. Increasing of the external pressure, acting on the porous medium, induces, on one hand, part of the fluid to flow out of the skeleton, on the other, part of the fluid to be segregated into not connected cavities of the solid matrix, when the external pressure overwhelms a critical value.
is $\bar{\epsilon}$ is the minimum of the function $f_+(\epsilon)$. Note, indeed, that by definition

$$ \mathcal{H}_+(\epsilon) = \Phi_{mm}(m_+(\epsilon), \epsilon)\Phi_{\epsilon\epsilon}(m_+(\epsilon), \epsilon) - \Phi_{m\epsilon}^2(m_+(\epsilon), \epsilon) \quad (25) $$

Recall that the function $m_+(\epsilon)$ is implicitly defined by the equation $\Phi_m(m_+(\epsilon), \epsilon) = 0$ and note that (11) is equivalent to $f_+(\epsilon) = -\Phi_\epsilon(m_+(\epsilon), \epsilon) + p$. By using the chain rule and the implicit function theorem we get that

$$ \begin{aligned} f'_+(\epsilon) &= -\Phi_{\epsilon\epsilon}(m_+(\epsilon), \epsilon) - \Phi_{\epsilon m}(m_+(\epsilon), \epsilon) m'_+(\epsilon) = -\Phi_{\epsilon\epsilon}(m_+(\epsilon), \epsilon) + \frac{\Phi_{\epsilon m}^2(m_+(\epsilon), \epsilon)}{\Phi_{mm}(m_+(\epsilon), \epsilon)} \quad (26) \\ &= -\frac{1}{\Phi_{mm}(m_+(\epsilon), \epsilon)} \mathcal{H}_+(\epsilon) \end{aligned} $$

where in the last step we have used (25).

The above equality ensures that $\bar{\epsilon} = \epsilon_c$. Finally, since, $\epsilon_+^1 < \epsilon_c < \epsilon_+^2$, we have that $(m_+^1, \epsilon_+^1)$ is a minimum of (5) while $(m_+^2, \epsilon_+^2)$ is a saddle.

# 4 Results

In this section a discussion on the behavior of the minima of the overall potential (5), when varying the external pressure, will be presented. A parametric analysis in dependence of the coefficients introduced in the Biot model, a and b, as well as of the additional coefficient $\alpha(a, b)$ multiplying the fourth order terms of (5) will be also developed.

According to the experimental results on fluidization of soils, cited in the introduction (see Nichols et al., 1994; Holcomb and Olsson, 2003; Holcomb et al., 2007; Vardoulakis, 2004a,b), the proposed model is capable to exhibit, in the presence of an applied consolidating external pressure, an additional stationary state in the $(\epsilon, m)$ plane beyond the one corresponding to classical consolidation. Increasing of the external pressure, acting on the porous medium, induces, on one hand, part of the fluid to flow out of the skeleton, on the other, part of the fluid to be segregated into not connected cavities of the solid matrix, when the external pressure overwhelms a critical value.
The arising of the additional state of equilibrium, associated with pore-fluid segregation, is explicitly depicted in Fig. (3), where the opposite of the total stress $\sigma = \partial\Psi/\partial\varepsilon$ and the fluid chemical potential $\mu = \partial\Psi/\partial m$ are plotted against the kinematical parameters $\varepsilon$ and $m$. The projection over the $(\varepsilon, m)$ plane of the curves obtained cutting the $-\sigma$ surface with the horizontal plane at $p$ and the $\mu$ surface with the $(\varepsilon, m)$ plane itself have more than one mutual intersection, each of them identifies one of the stationary points of the overall potential (5). In particular the dotted (dashed) lines in Fig. (3) (Fig. (4)) correspond to the solutions of $\Phi_m = 0$, on the other hand the solid ones correspond to the solution of $\Phi_\varepsilon = 0$, for different values of $p$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Cutting the $-\sigma$ surface with the horizontal plane at $p$ identifies the root locus of $\Phi_\varepsilon = 0$. In particular the two solid lines identify the root loci associated with the critical value of the applied pressure and that associated to a larger value of $p$. The dotted lines describe the root locus of $\Phi_m = 0$, obtained cutting the $\mu$ surface with the $(\varepsilon, m)$ plane.</figcaption></figure>

As already noticed, the critical pressure $p_c$ is the smallest value of the consolidating loading for which the pore-fluid segregation solution shows up.

With another point of view, the stationarity conditions stated by eq. (6) and (7) describe two
The arising of the additional state of equilibrium, associated with pore-fluid segregation, is explicitly depicted in Fig. (3), where the opposite of the total stress $\sigma = \partial\Psi/\partial\varepsilon$ and the fluid chemical potential $\mu = \partial\Psi/\partial m$ are plotted against the kinematical parameters $\varepsilon$ and $m$. The projection over the $(\varepsilon, m)$ plane of the curves obtained cutting the $-\sigma$ surface with the horizontal plane at $p$ and the $\mu$ surface with the $(\varepsilon, m)$ plane itself have more than one mutual intersection, each of them identifies one of the stationary points of the overall potential (5). In particular the dotted (dashed) lines in Fig. (3) (Fig. (4)) correspond to the solutions of $\Phi_m = 0$, on the other hand the solid ones correspond to the solution of $\Phi_\varepsilon = 0$, for different values of $p$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Cutting the $-\sigma$ surface with the horizontal plane at $p$ identifies the root locus of $\Phi_\varepsilon = 0$. In particular the two solid lines identify the root loci associated with the critical value of the applied pressure and that associated to a larger value of $p$. The dotted lines describe the root locus of $\Phi_m = 0$, obtained cutting the $\mu$ surface with the $(\varepsilon, m)$ plane.</figcaption></figure>

As already noticed, the critical pressure $p_c$ is the smallest value of the consolidating loading for which the pore-fluid segregation solution shows up.

With another point of view, the stationarity conditions stated by eq. (6) and (7) describe two
<figure><img src="image_5.png" /><figcaption>Figure 4: This is a detail of the projection on the $(\epsilon, m)$ plane of the intersection of $-\sigma$ and $\mu$ with the horizontal plane at $p$ and zero, respectively. In particular $(\epsilon_i^c, m_i^c)$, $i = 1, 2$, identify the two stationary points of (5) when $p = p_c$, while $(\epsilon_i^{sc}, m_i^{sc})$, $i = 1, 2, 3$, identify the three stationary points of (5) when $p > p_c$ (super-critical conditions).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: This is a detail of the projection on the $(\epsilon, m)$ plane of the intersection of $-\sigma$ and $\mu$ with the horizontal plane at $p$ and zero, respectively. In particular $(\epsilon_i^c, m_i^c)$, $i = 1, 2$, identify the two stationary points of (5) when $p = p_c$, while $(\epsilon_i^{sc}, m_i^{sc})$, $i = 1, 2, 3$, identify the three stationary points of (5) when $p > p_c$ (super-critical conditions).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: The dotted lines correspond to the intersection between $-\sigma$ and the root locus of $\Phi_m = 0$, the solid lines are the projections of these curves on the stress–strain plane, i.e. the $f_1$ and $f_\pm$ curves.</figcaption></figure>

different surfaces in the $(\varepsilon, m)$ plane, the intersection of which gives rise to a curve in the three–dimensional space, the dotted curve in Fig. (5), playing the role of the curve of the fluid pressure in terms of the specific volume, in the case of Van der Waals’ model of liquid–vapor coexistence. It is worth to notice that the projection of these curves on the $(-\sigma, \varepsilon)$ plane provides the plot of the $f_1$ and $f_\pm$ functions, depicted in Fig. (2). As in the case of Van der Waals’ model sectioning this curve in the three dimensional space with an horizzontal plane at $p$ identifies one stationary (equilibrium) point, if $p < p_c$, or three stationary points, if $p \ge p_c$, two of which corresponds to classical consolidation and pore–fluid segregation solutions, respectively.

Fig. (6) explicitly shows the parametrization with pressure of the states of equilibrium: when the applied external pressure reaches the critical value $p_c$ a new state of equilibrium appears in which the variation of the fluid mass and the negative strain increase. This is what is generally called fluid–segregation: under the applied external pressure the porous material undergoes to consolidation, however part of the fluid remains trapped into non–connected cavities of the skeleton, the formation of which is possibly due to the thinnering of duct connections of the pore–space.

Note that all the pictures have been drawn for the parameters specified in the caption of Fig.(2). However it is important to stress that, according to the general discussion of §3, no variation in
<figure><img src="image_6.png" /><figcaption>Figure 5: The dotted lines correspond to the intersection between $-\sigma$ and the root locus of $\Phi_m = 0$, the solid lines are the projections of these curves on the stress–strain plane, i.e. the $f_1$ and $f_\pm$ curves.</figcaption></figure>

different surfaces in the $(\varepsilon, m)$ plane, the intersection of which gives rise to a curve in the three–dimensional space, the dotted curve in Fig. (5), playing the role of the curve of the fluid pressure in terms of the specific volume, in the case of Van der Waals’ model of liquid–vapor coexistence. It is worth to notice that the projection of these curves on the $(-\sigma, \varepsilon)$ plane provides the plot of the $f_1$ and $f_\pm$ functions, depicted in Fig. (2). As in the case of Van der Waals’ model sectioning this curve in the three dimensional space with an horizzontal plane at $p$ identifies one stationary (equilibrium) point, if $p < p_c$, or three stationary points, if $p \ge p_c$, two of which corresponds to classical consolidation and pore–fluid segregation solutions, respectively.

Fig. (6) explicitly shows the parametrization with pressure of the states of equilibrium: when the applied external pressure reaches the critical value $p_c$ a new state of equilibrium appears in which the variation of the fluid mass and the negative strain increase. This is what is generally called fluid–segregation: under the applied external pressure the porous material undergoes to consolidation, however part of the fluid remains trapped into non–connected cavities of the skeleton, the formation of which is possibly due to the thinnering of duct connections of the pore–space.

Note that all the pictures have been drawn for the parameters specified in the caption of Fig.(2). However it is important to stress that, according to the general discussion of §3, no variation in
<figure><img src="image_7.png" /><figcaption>Figure 6: The graph of $f_1$ corresponds to the classical behavior of consolidating soils, the negative strain progressively increases and the mass of the fluid progressively decreases when the external pressure raises up; the graph of $f_+$, on the other hand, describes the phenomenon of fluid-segregation, the dashed line (graph of $f_-$) corresponds to the saddle point of the overall potential.</figcaption></figure>

the main features of the considered phase transition, from the classical consolidating soil phase towards the pore-fluid segregated one, can arise when tuning the constitutive parameters $a$, $b$ and $\alpha$. The only constraint on the showing up of the new equilibrium being the impenetrability of the skeleton ($\varepsilon > -1/2$). It is worth to notice that the absolute values of both the free energy of the pure fluid and that associated with solid-fluid coupling in the overall free energy (5)

$$ \Psi_{\text{ff}}(\varepsilon, m) := 1/4\alpha m^4 + 1/2am^2 \quad \text{and} \quad \Psi_{\text{sf}}(\varepsilon, m) := -2/3\alpha bm^3\varepsilon + 1/2\alpha b^2m^2\varepsilon^2 - abm\varepsilon \quad (27) $$

progressively decrease for increasing values of the applied external pressure, along the pore-fluid segregating equilibrium path, see Fig. (7). Apparently this is due to the increasing incapability of the pore fluid to flow out of the skeleton. Conversely, the free energy of the skeleton, $\Psi_{\text{ss}} := 1/2(1+ab^2)\varepsilon^2$, along the pore-fluid segregating equilibrium is larger than the one associated to classical consolidation, the trapped fluid inducing a virtual stiffening of the solid matrix.

Phase transition is therefore completely characterized in terms of the dependence of the critical pressure on the ratio between the Biot modulus and the bulk modulus of the skeleton, $a$, and the Biot coefficient $b$. In the frame of physically meaningful assumptions on these parameters
<figure><img src="image_7.png" /><figcaption>Figure 6: The graph of $f_1$ corresponds to the classical behavior of consolidating soils, the negative strain progressively increases and the mass of the fluid progressively decreases when the external pressure raises up; the graph of $f_+$, on the other hand, describes the phenomenon of fluid-segregation, the dashed line (graph of $f_-$) corresponds to the saddle point of the overall potential.</figcaption></figure>

the main features of the considered phase transition, from the classical consolidating soil phase towards the pore-fluid segregated one, can arise when tuning the constitutive parameters $a$, $b$ and $\alpha$. The only constraint on the showing up of the new equilibrium being the impenetrability of the skeleton ($\epsilon > -1/2$). It is worth to notice that the absolute values of both the free energy of the pure fluid and that associated with solid-fluid coupling in the overall free energy (5)

$$ \Psi_{\text{ff}}(\epsilon, m) := 1/4\alpha m^4 + 1/2am^2 \quad \text{and} \quad \Psi_{\text{sf}}(\epsilon, m) := -2/3\alpha bm^3\epsilon + 1/2\alpha b^2m^2\epsilon^2 - abm\epsilon \quad (27) $$

progressively decrease for increasing values of the applied external pressure, along the pore-fluid segregating equilibrium path, see Fig. (7). Apparently this is due to the increasing incapability of the pore fluid to flow out of the skeleton. Conversely, the free energy of the skeleton, $\Psi_{\text{ss}} := 1/2(1 + ab^2)\epsilon^2$, along the pore-fluid segregating equilibrium is larger than the one associated to classical consolidation, the trapped fluid inducing a virtual stiffening of the solid matrix.

Phase transition is therefore completely characterized in terms of the dependence of the critical pressure on the ratio between the Biot modulus and the bulk modulus of the skeleton, $a$, and the Biot coefficient $b$. In the frame of physically meaningful assumptions on these parameters
<figure><img src="image_8.png" /><figcaption>Figure 7: The solid lines indicate the *pure solid*, the *pure fluid* and the *coupling* energies associated with the soil consolidating equilibrium path; the dashed ones indicate the same energy contributions associated with the pore-fluid segregation equilibrium.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Picture of the variation of the critical pressure with the Biot coefficient *b*, parametrized by *a*; both *a* and *b* range in the open interval (0, 1).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: The solid lines indicate the *pure solid*, the *pure fluid* and the *coupling* energies associated with the soil consolidating equilibrium path; the dashed ones indicate the same energy contributions associated with the pore-fluid segregation equilibrium.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Picture of the variation of the critical pressure with the Biot coefficient *b*, parametrized by *a*; both *a* and *b* range in the open interval (0, 1).</figcaption></figure>
(a ≤ 1 and b ≤ 1) the critical pressure $p_c$ exhibits the behavior depicted in Fig. (8). Increasing a, keeping b fixed, implies the critical pressure to increase and therefore the phase transition to occur for higher values of the consolidating pressure. Indeed increasing a means increasing of the bulk modulus of the fluid; the inverse of Biot's modulus is an affine function of the inverse of the fluid bulk modulus, with respect to that of the solid. The fluid becomes more stiff and therefore does not allow the applied external pressure to force shrinkage of the ducts inside the porous medium; as a consequence flowing of the fluid out of the solid matrix is enhanced.

Increasing b, keeping a fixed, on the other hand, implies the critical pressure to decrease and therefore the phase transition to occur for lower levels of the consolidating loading. This is due to the fact that increasing b enhances solid-fluid coupling.

Apparently no physical reasoning can be developed in order to discuss the behavior of the system when tuning α, however it is clear that only for sufficiently high values of this parameter phase transition can be appreciated.

## References

M. Biot. General theory of three dimensional consolidation. *Journal of Applied Physics*, 12 (2): 155-164, 1941.

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

R. Bürger, W.L. Wendland, and F. Concha. A mathematical model for sedimentation-consolidation processes. *Zeitschrift für Angewandte Mathematik und Mechanik*, 80:79-92, 2000.

O. Coussy. *Poromechanics*. John Wiley & Sons, Ltd, 2004.

C.W. Cryer. A comparison of the three-dimensional consolidation theories of Biot and Terzaghi. *Quart. J. Mech. and Applied Math.*, 16 (4):401-412, 1963.
(a ≤ 1 and b ≤ 1) the critical pressure $p_c$ exhibits the behavior depicted in Fig. (8). Increasing a, keeping b fixed, implies the critical pressure to increase and therefore the phase transition to occur for higher values of the consolidating pressure. Indeed increasing a means increasing of the bulk modulus of the fluid; the inverse of Biot's modulus is an affine function of the inverse of the fluid bulk modulus, with respect to that of the solid. The fluid becomes more stiff and therefore does not allow the applied external pressure to force shrinkage of the ducts inside the porous medium; as a consequence flowing of the fluid out of the solid matrix is enhanced.

Increasing b, keeping a fixed, on the other hand, implies the critical pressure to decrease and therefore the phase transition to occur for lower levels of the consolidating loading. This is due to the fact that increasing b enhances solid-fluid coupling.

Apparently no physical reasoning can be developed in order to discuss the behavior of the system when tuning α, however it is clear that only for sufficiently high values of this parameter phase transition can be appreciated.

## References

M. Biot. General theory of three dimensional consolidation. *Journal of Applied Physics*, 12 (2): 155-164, 1941.

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

R. Bürger, W.L. Wendland, and F. Concha. A mathematical model for sedimentation-consolidation processes. *Zeitschrift für Angewandte Mathematik und Mechanik*, 80:79-92, 2000.

O. Coussy. *Poromechanics*. John Wiley & Sons, Ltd, 2004.

C.W. Cryer. A comparison of the three-dimensional consolidation theories of Biot and Terzaghi. *Quart. J. Mech. and Applied Math.*, 16 (4):401-412, 1963.
D. Holcomb and W.A. Olsson. Compaction localization and fluid flow. *Journal of Geophysical Research*, 108 (B6):2290, 2003.

D. Holcomb, J. W. Rudnicki, K. A. Issen, and K. Sternlof. Compaction localization in the earth and the laboratory: State of the research and research directions. *Acta Geotechnica*, 2 (1):1-15, 2007.

D. Kolymbas. Behaviour of liquified sand. *Philosophical Transactions of the Royal Society of London A*, 356:2609-2622, 1998.

J. Mandel. Consolidation des sols (étude mathématique). *Géotechnique*, 3:287-299, 1953.

P. N. Mollema and M. A. Antonellini. Compaction bands: a structural analog for anti-mode I cracks in aeolian sandstone. *Tectonophysics*, 267:2092-228, 1996.

R.J. Nichols, R.S.J. Sparks, and C.J.N. Wilson. Experimental studies of the fluidization of layered sediments and the formation of fluid escape structures. *Sedimentology*, 41:233-253, 1994.

G. Sciarra, F. dell'Isola, N. Ianiro, and A. Madeo. A variational deduction of second gradient poroelasticity: Part I: General theory. *Journal of Mechanics of Materials and Structures*, 3 (3): 507-526, 2008.

K. Terzaghi. *Theoretical Soil Mechanics*. John Wiley & Sons, Ltd, 1946.

I. Vardoulakis. Fluidization in artesian flow conditions: hydromechanically stable granular media. *Géotechnique*, 54:117-130, 2004a.

I. Vardoulakis. Fluidization in artesian flow conditions: hydromechanically unstable granular media. *Géotechnique*, 54:165-177, 2004b.
D. Holcomb and W.A. Olsson. Compaction localization and fluid flow. *Journal of Geophysical Research*, 108 (B6):2290, 2003.

D. Holcomb, J. W. Rudnicki, K. A. Issen, and K. Sternlof. Compaction localization in the earth and the laboratory: State of the research and research directions. *Acta Geotechnica*, 2 (1):1-15, 2007.

D. Kolymbas. Behaviour of liquified sand. *Philosophical Transactions of the Royal Society of London A*, 356:2609-2622, 1998.

J. Mandel. Consolidation des sols (étude mathématique). *Géotechnique*, 3:287-299, 1953.

P. N. Mollema and M. A. Antonellini. Compaction bands: a structural analog for anti-mode I cracks in aeolian sandstone. *Tectonophysics*, 267:2092-228, 1996.

R.J. Nichols, R.S.J. Sparks, and C.J.N. Wilson. Experimental studies of the fluidization of layered sediments and the formation of fluid escape structures. *Sedimentology*, 41:233-253, 1994.

G. Sciarra, F. dell'Isola, N. Ianiro, and A. Madeo. A variational deduction of second gradient poroelasticity: Part I: General theory. *Journal of Mechanics of Materials and Structures*, 3 (3): 507-526, 2008.

K. Terzaghi. *Theoretical Soil Mechanics*. John Wiley & Sons, Ltd, 1946.

I. Vardoulakis. Fluidization in artesian flow conditions: hydromechanically stable granular media. *Géotechnique*, 54:117-130, 2004a.

I. Vardoulakis. Fluidization in artesian flow conditions: hydromechanically unstable granular media. *Géotechnique*, 54:165-177, 2004b.