# Allen–Cahn and Cahn–Hilliard-like equations for dissipative dynamics of saturated porous media

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

► To cite this version:

Emilio N.M. Cirillo, Nicoletta Ianiro, Giulio Sciarra. Allen–Cahn and Cahn–Hilliard-like equations for dissipative dynamics of saturated porous media. Journal of the Mechanics and Physics of Solids, 2013, 61 (2), 10.1016/j.jmps.2012.08.014. hal-01362182

HAL Id: hal-01362182

https://hal.science/hal-01362182

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
# Allen–Cahn and Cahn–Hilliard–like equations for dissipative dynamics of saturated porous media

Emilio N.M. Cirillo<sup>a,*</sup>, Nicoletta Ianiro<sup>b</sup>, Giulio Sciarra<sup>c</sup>

<sup>a</sup>Dipartimento di Scienze di Base e Applicate per l'Ingegneria, Sapienza Università di Roma, via A. Scarpa 16, I-00161, Roma, Italy.

<sup>b</sup>Dipartimento di Scienze di Base e Applicate per l'Ingegneria, Sapienza Università di Roma, via A. Scarpa 16, I-00161, Roma, Italy.

<sup>c</sup>Dipartimento di Ingegneria Chimica Materiali Ambiente, Sapienza Università di Roma, via Eudossiana 18, 00184 Roma, Italy

## Abstract

We consider a saturated porous medium in the solid-fluid segregation regime under the effect of an external pressure applied on the solid constituent. We prove that, depending on the dissipation mechanism, the dynamics is described either by a Cahn-Hilliard or by an Allen-Cahn-like equation. More precisely, when the dissipation is modeled via the Darcy law we find that, provided the solid deformation and the fluid density variations are small, the evolution equation is very similar to the Cahn-Hilliard one. On the other hand, when only the Stokes dissipation term is considered, we find that the evolution is governed by an Allen-Cahn-like equation. We use this theory to describe the formation of interfaces inside porous media. We consider a recently developed model proposed to study the solid-liquid segregation in consolidation and we give a complete description of the formation of an interface between the fluid-rich and the fluid-poor phase.

**Keywords:** phase transformation (A), porous material (B), finite differences (C), energy methods (C)

## 1. Introduction

Many interesting swelling and shrinking phenomena are observed when a porous medium deforms as a consequence of the variation of the fluid mass content. The behavior of the system can be controlled via different external parameters, e.g. the fluid chemical potential [Gelb et al. 1999, Convertino et al. 2003, Convertino et al. 2002] or a mechanical pressure exerted on the solid component [Cirillo et al. 2011, Cirillo et al. 2010].

By focusing the attention on this last class of phenomena we are led to consider the so-called consolidation problem. Despite the existence of classical and well known theories [von Terzaghi 1946, Biot 1941, Cryer 1963, Mandel 1953],

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: emilio.cirillo@uniroma1.it (Emilio N.M. Cirillo), nicoletta.ianiro@uniroma1.it (Nicoletta Ianiro), giulio.sciarra@uniroma1.it (Giulio Sciarra)</footnote>
Allen–Cahn and Cahn–Hilliard–like equations for
dissipative dynamics of saturated porous media

Emilio N.M. Cirillo<sup>a,*</sup>, Nicoletta Ianiro<sup>b</sup>, Giulio Sciarra<sup>c</sup>

<sup>a</sup>Dipartimento di Scienze di Base e Applicate per l'Ingegneria, Sapienza Università di Roma, via A. Scarpa 16, I-00161, Roma, Italy.

<sup>b</sup>Dipartimento di Scienze di Base e Applicate per l'Ingegneria, Sapienza Università di Roma, via A. Scarpa 16, I-00161, Roma, Italy.

<sup>c</sup>Dipartimento di Ingegneria Chimica Materiali Ambiente, Sapienza Università di Roma, via Eudossiana 18, 00184 Roma, Italy

**Abstract**

We consider a saturated porous medium in the solid-fluid segregation regime under the effect of an external pressure applied on the solid constituent. We prove that, depending on the dissipation mechanism, the dynamics is described either by a Cahn-Hilliard or by an Allen-Cahn-like equation. More precisely, when the dissipation is modeled via the Darcy law we find that, provided the solid deformation and the fluid density variations are small, the evolution equation is very similar to the Cahn-Hilliard one. On the other hand, when only the Stokes dissipation term is considered, we find that the evolution is governed by an Allen-Cahn-like equation. We use this theory to describe the formation of interfaces inside porous media. We consider a recently developed model proposed to study the solid-liquid segregation in consolidation and we give a complete description of the formation of an interface between the fluid-rich and the fluid-poor phase.

**Keywords:** phase transformation (A), porous material (B), finite differences (C), energy methods (C)

**1. Introduction**

Many interesting swelling and shrinking phenomena are observed when a porous medium deforms as a consequence of the variation of the fluid mass content. The behavior of the system can be controlled via different external parameters, e.g. the fluid chemical potential [Gelb et al. 1999, Convertino et al. 2003, Convertino et al. 2002] or a mechanical pressure exerted on the solid component [Cirillo et al. 2011, Cirillo et al. 2010].

By focusing the attention on this last class of phenomena we are led to consider the so-called consolidation problem. Despite the existence of classical and well known theories [von Terzaghi 1946, Biot 1941, Cryer 1963, Mandel 1953],

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: emilio.cirillo@uniroma1.it (Emilio N.M. Cirillo),
nicoletta.ianiro@uniroma1.it (Nicoletta Ianiro), giulio.sciarra@uniroma1.it (Giulio
Sciarra)</footnote>
the dilatant/contractant behavior of porous solids needs, however, deeper investigation. We refer, in particular, to undrained conditions, i.e., to the case in which the loading rate does not allow pore-water pressure dissipation. In these conditions, the coupling between shearing and dilatancy typically induces the formation of strain localization bands. In particular shear bands with transverse size depending on the texture and the constitutive properties of the material, see e.g. [Besuelle 2001], are seen. For instance, in simple shear test, densely or loosely packed granular media exhibit dilatancy or compaction behavior. Moreover due to porosity change, the fluid can migrate through the pores and eventually remain segregated, possibly enhancing localised overpressurization and fluidization of the soil, see e.g. [Kolymbas 1994, Nichols et al. 1994].

Local dilatancy is typically modeled in the framework of (poro-)plasticity by assuming that irreversible deformations take place inside the porous medium, see e.g. the general treatises [Coussy 2004, Coussy 2010]. Special important applications to multi-phase porous media and geomaterials have been considered in the recent literature [Zhang et al. 1999, Collin et al. 2006, Chambon et al. 2001, Andrade et al. 2011] by approaching plasticity with strain gradient [Fleck et al. 1997, Fleck et al. 2001] or multi-scale analysis [Christoffersen et al. 1981, Nemat-Nasser et al. 1993].

In the papers [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] we have attacked this problem from the point of view of bifurcation theory and we have shown that it is possible to describe interesting phenomena (still in the range of non-linear elasticity) taking place when the pressure exerted on the solid exceeds a suitable limiting value. More precisely, we have proven that enhanced models based on gradient elasticity allow for describing the transition between phases of the porous medium associated with different fluid content, say a fluid-rich and a fluid-poor phase. In other words the theory proposed in those papers is able to describe the occurrence of segregation and local overpressurization of the pore-water, see the above quoted references for more details.

Experiments demonstrate that the presence of a pressurized fluid-rich phase at the base of a layered granular material can give rise to an unstable (with respect to possible fluidization) configuration of the system [Kolymbas 1998]. This result is achieved in laboratory, see [Nichols et al. 1994], considering a fluidized column test setup where a fluid is forced to flow through a saturated sample from the bottom. By tuning the velocity of the fluid, the drag force acting on the solid grains, and hence the possible unbalance of gravity, is controlled [Vardoulakis 2004a, Vardoulakis 2004b]. The model in [Cirillo et al. 2009] is able to explain the formation of fluid-rich regions inside a porous material via a sort of segregation process, due to different initial data and parametrized by the pressure applied to the solid; note that the control parameter in the experiment quoted above is, instead, the velocity of the injected fluid. Different dissipation mechanisms, driven by the microstructural properties of the material, give rise to different deformation paths and fluid mass content evolutions.

In this paper we show that, depending on the dissipative mechanism which is taken into account, the evolution is described either by an Allen-Cahn or by a Cahn-Hilliard-like system of partial differential equations, which indeed play an important role in the theory of phase ordering dynamics [Bray 1994, Langer 1992, Eyre 1993]. When a system is quenched from the homogeneous phase into a broken-symmetry one, think to a ferromagnet or to a gas abruptly cooled below their critical temperature, the two phases have to separate and order has to increase throughout the system via domain coarsening. This phe-
the dilatant/contractant behavior of porous solids needs, however, deeper in-
vestigation. We refer, in particular, to undrained conditions, i.e., to the case
in which the loading rate does not allow pore-water pressure dissipation. In
these conditions, the coupling between shearing and dilatancy typically induces
the formation of strain localization bands. In particular shear bands with tran-
severe size depending on the texture and the constitutive properties of the
material, see e.g. [Besuelle 2001], are seen. For instance, in simple shear test,
densely or loosely packed granular media exhibit dilatancy or compaction behav-
ior. Moreover due to porosity change, the fluid can migrate through the pores
and eventually remain segregated, possibly enhancing localised overpressuriza-
tion and fluidization of the soil, see e.g. [Kolymbas 1994, Nichols et al. 1994].

Local dilatancy is typically modeled in the framework of (poro-)plasticity by
assuming that irreversible deformations take place inside the porous medium, see
e.g. the general treatises [Coussy 2004, Coussy 2010]. Special important appli-
cations to multi-phase porous media and geomaterials have been considered in
the recent literature [Zhang et al. 1999, Collin et al. 2006, Chambon et al. 2001,
Andrade et al. 2011] by approaching plasticity with strain gradient [Fleck et al. 1997,
Fleck et al. 2001] or multi-scale analysis [Christoffersen et al. 1981, Nemat-Nasser et al. 1993].

In the papers [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] we have attacked this problem from the point of view of bifurcation theory and we have shown that it is possible to describe interesting phenomena (still in the range of non-linear elasticity) taking place when the pressure exerted on the solid exceeds a suitable limiting value. More precisely, we have proven that enhanced models based on gradient elasticity allow for describing the transition between phases of the porous medium associated with different fluid content, say a fluid-rich and a fluid-poor phase. In other words the theory proposed in those papers is able to describe the occurrence of segregation and local overpres- surization of the pore-water, see the above quoted references for more details.

Experiments demonstrate that the presence of a pressurized fluid-rich phase at the base of a layered granular material can give rise to an unstable (with respect to possible fluidization) configuration of the system [Kolymbas 1998]. This result is achieved in laboratory, see [Nichols et al. 1994], considering a fluidized column test setup where a fluid is forced to flow through a saturated sample from the bottom. By tuning the velocity of the fluid, the drag force acting on the solid grains, and hence the possible unbalance of gravity, is controlled [Vardoulakis 2004a, Vardoulakis 2004b]. The model in [Cirillo et al. 2009] is able to explain the formation of fluid-rich regions inside a porous material via a sort of segregation process, due to different initial data and parametrized by the pressure applied to the solid; note that the control parameter in the experiment quoted above is, instead, the velocity of the injected fluid. Different dissipation mechanisms, driven by the microstructural properties of the material, give rise to different deformation paths and fluid mass content evolutions.

In this paper we show that, depending on the dissipative mechanism which
is taken into account, the evolution is described either by an Allen–Cahn or
by a Cahn–Hilliard–like system of partial differential equations, which indeed
play an important role in the theory of phase ordering dynamics [Bray 1994,
Langer 1992, Eyre 1993]. When a system is quenched from the homogeneous
phase into a broken–symmetry one, think to a ferromagnet or to a gas abruptly
cooled below their critical temperature, the two phases have to separate and
order has to increase throughout the system via domain coarsening. This phe-
nomenon can be described via a field $u$ on the physical space $\Omega \subset \mathbb{R}^d$, with $d$ the physical dimension, having the proper physical interpretation, for instance local magnetization in ferromagnets, density in liquid-vapor systems, and concentration in alloys. A possible model for the evolution of this field is the Allen-Cahn equation

$$ \frac{\partial u}{\partial t} = \varepsilon^2 \Delta u - W'(u) \quad (1) $$

where $\varepsilon$ is a positive constant and $W(u)$ a double well regular function. The Allen-Cahn equation, also called the time-dependent Ginzburg-Landau equation, was introduced in [Allen et al. 1979] to describe the motion of anti-phase boundaries in crystalline solids. In this context $u$ represents the concentration of one of the two components of the alloy and the parameter $\varepsilon^2$ the interface width.

The Allen-Cahn equation is not suited to describe phase separation when the order parameter is conserved, namely, when the integral of the field $u$ on the whole space $\Omega$ is constant. On the other hand this is possible in the framework of the Cahn-Hilliard equation

$$ \frac{\partial u}{\partial t} = -\Delta(\varepsilon^2 \Delta u - W'(u)) \quad (2) $$

by assuming $\partial_{\nu}u = \partial_{\nu}\Delta u = 0$ on $\partial\Omega$, that is by requiring the derivative of $u$ and that of the Laplacian of $u$ orthogonal to the boundary to vanish on the boundary $\partial\Omega$ itself. The Cahn-Hilliard equation was firstly introduced to describe spinodal decomposition in alloys [Cahn 1961].

A straightforward way to derive the Allen-Cahn and the Cahn-Hilliard equations is that of assuming the evolution of the field $u$ to be governed by a gradient equation [Fife 2002, Alikakos et al. 1991] $\partial u/\partial t = -\delta F/\delta u$ associated with the Landau energy functional

$$ F(u) := \int_{\Omega} \left[ \frac{1}{2} \varepsilon^2 \|\nabla u\|^2 + W(u) \right] dx \quad (3) $$

The functional above has a clear physical interpretation: the term $W$ has two minima corresponding to the two phases, while the gradient-squared term associates an energy cost with an interface between the two phases. If no constraint to the total value of the order parameter is imposed, it is possible to compute the gradient of the Landau functional in the Hilbert space $L^2(\Omega)$ to get the Allen-Cahn equation. When the integral of the field $u$ is assumed to be constant throughout the evolution, computing the gradient in the $L^2(\Omega)$ results into a non-local evolution equation, while by using the Hilbert space $H^{-1}(\Omega)$ the Cahn-Hilliard equation is found.

The above described approach is general and applies to any physical interpretation of the two equations. When the discussion is bounded to particular physical systems, phenomenological approaches can also be used requiring the a priori specification of the constitutive equations. The related physical literature is huge and, among others, we refer to [Gurtin 1996]. In the context of binary fluids, for example, $u$ is the fluid concentration. Assume the continuity equation in the form $u_t + \nabla \cdot J = 0$, where $J$ is the fluid current prescribed in terms of the chemical potential $\mu$ by Fick's law $J = -k\nabla\mu$, with $k$ a positive constant.
nomenon can be described via a field $u$ on the physical space $\Omega \subset \mathbb{R}^d$, with $d$ the physical dimension, having the proper physical interpretation, for instance local magnetization in ferromagnets, density in liquid-vapor systems, and concentration in alloys. A possible model for the evolution of this field is the Allen-Cahn equation

$$ \frac{\partial u}{\partial t} = \varepsilon^2 \Delta u - W'(u) \quad (1) $$

where $\varepsilon$ is a positive constant and $W(u)$ a double well regular function. The Allen-Cahn equation, also called the time-dependent Ginzburg-Landau equation, was introduced in [Allen et al. 1979] to describe the motion of anti-phase boundaries in crystalline solids. In this context $u$ represents the concentration of one of the two components of the alloy and the parameter $\varepsilon^2$ the interface width.

The Allen-Cahn equation is not suited to describe phase separation when the order parameter is conserved, namely, when the integral of the field $u$ on the whole space $\Omega$ is constant. On the other hand this is possible in the framework of the Cahn-Hilliard equation

$$ \frac{\partial u}{\partial t} = -\Delta(\varepsilon^2 \Delta u - W'(u)) \quad (2) $$

by assuming $\partial_\nu u = \partial_\nu \Delta u = 0$ on $\partial\Omega$, that is by requiring the derivative of $u$ and that of the Laplacian of $u$ orthogonal to the boundary to vanish on the boundary $\partial\Omega$ itself. The Cahn-Hilliard equation was firstly introduced to describe spinodal decomposition in alloys [Cahn 1961].

A straightforward way to derive the Allen-Cahn and the Cahn-Hilliard equations is that of assuming the evolution of the field $u$ to be governed by a gradient equation [Fife 2002, Alikakos et al. 1991] $\partial u/\partial t = -\delta F/\delta u$ associated with the Landau energy functional

$$ F(u) := \int_{\Omega} \left[ \frac{1}{2} \varepsilon^2 \|\nabla u\|^2 + W(u) \right] dx \quad (3) $$

The functional above has a clear physical interpretation: the term $W$ has two minima corresponding to the two phases, while the gradient-squared term associates an energy cost with an interface between the two phases. If no constraint to the total value of the order parameter is imposed, it is possible to compute the gradient of the Landau functional in the Hilbert space $L^2(\Omega)$ to get the Allen-Cahn equation. When the integral of the field $u$ is assumed to be constant throughout the evolution, computing the gradient in the $L^2(\Omega)$ results into a non-local evolution equation, while by using the Hilbert space $H^{-1}(\Omega)$ the Cahn-Hilliard equation is found.

The above described approach is general and applies to any physical interpretation of the two equations. When the discussion is bounded to particular physical systems, phenomenological approaches can also be used requiring the a priori specification of the constitutive equations. The related physical literature is huge and, among others, we refer to [Gurtin 1996]. In the context of binary fluids, for example, $u$ is the fluid concentration. Assume the continuity equation in the form $u_t + \nabla \cdot J = 0$, where $J$ is the fluid current prescribed in terms of the chemical potential $\mu$ by Fick's law $J = -k\nabla\mu$, with $k$ a positive constant.
By letting $\mu$ be the functional derivative of $F$, one gets

$$ \frac{\partial u}{\partial t} = -\nabla \cdot J = -\nabla \cdot (-k\nabla\mu) = k\Delta\mu = k\Delta\frac{\delta F}{\delta u} = k\Delta(-\epsilon^2\Delta u + W'(u)) $$

which is the Cahn-Hilliard equation.

In our problem the Allen-Cahn and the Cahn-Hilliard-like equations appear in a completely natural way by using the variational description of continuum mechanics. By following [Coussy 2004, Sciarra et al. 2008] we describe the porous material via two fields, the strain $\epsilon$ and the fluid density (measured with respect to the solid reference volume) $m$. We then assume that the possible motions of the system are those such that the variation of the action equals the opposite of the time integral of the virtual work of the dissipative forces corresponding to the variation of the fields. Solving the variational problem we get the equations of motion.

The explicit form of the equations of motion depends not only on the form of the action but also on the way in which the dissipation phenomenon is modeled [Nield 2007, Lopatnikov et al. 2004]. Assuming the potential energy to be quadratic in the first derivatives of the strain and of the fluid density, the evolution is described by an Allen-Cahn-like equation for the field $m$ provided that dissipative forces are assumed to be proportional to the second derivative of the velocity of the fluid. We refer to this assumption as to the Stokes dissipation. On the other hand, we find that the evolution is described by a Cahn-Hilliard-like equation for the field $m$ provided dissipation is modeled via the Darcy law, namely, the dissipative forces are proportional to the velocity of the fluid.

It is possible to give a nice physical interpretation to the fact that the evolution is described by the Allen-Cahn or the Cahn-Hilliard-like equations depending on the way in which dissipation is put into the game. With homogenization arguments it has been proven that when the characteristic length at the local ($l$) and at the macroscopic ($L$) scales separate, say $l/L \ll 1$, then the resistance experienced by the fluid when flowing through the porous material can be described by the Darcy law [Ene et al. 1975]. On the other hand when the separation of scales is poor, the Darcy limit unfairly approximates the flow, since the macroscopic characteristic length is of the same order of magnitude as the pore characteristic size [Auriault et al. 2005]. In other words in this case the solid grains have to be considered as microporous themselves.

Moreover, we note that the steady Stokes equations, which describe the behavior of the fluid at the microscopic level, can be upscaled to a set of equations describing a non-Stokesian flow only if the size of the solid obstacles to the fluid motion have a critical size with respect to the inter-obstacle distance, see for more details [Allaire 1997, Theorem 3.1]. This condition provide the so-called Brinkman limit. Conversely when the ostables are too small the homogenized equations coincide with the Stokes equations. Thus the Allen-Cahn-like equations, which correspond to a purely Stokesian flow, are controlled by the viscosity of the fluid; on the other hand the Cahn-Hilliard-like equations, which correspond to a Darcy flow, are controlled by the permeability of the solid skeleton. In particular the Cahn-Hilliard-like model reads, within the framework of soil consolidation theory, as a generalization to strain gradient materials of the one-dimensional Terzaghi problem, see [Madeo et al. 2008].

The paper is organized as follows. We first introduce the model in Section 2. In Section 3 we assume small variations of the fields with respect to some ref-
By letting $\mu$ be the functional derivative of $F$, one gets

$$ \frac{\partial u}{\partial t} = -\nabla \cdot J = -\nabla \cdot (-k\nabla\mu) = k\Delta\mu = k\Delta\frac{\delta F}{\delta u} = k\Delta(-\epsilon^2\Delta u + W'(u)) $$

which is the Cahn-Hilliard equation.

In our problem the Allen-Cahn and the Cahn-Hilliard-like equations appear in a completely natural way by using the variational description of continuum mechanics. By following [Coussy 2004, Sciarra et al. 2008] we describe the porous material via two fields, the strain $\epsilon$ and the fluid density (measured with respect to the solid reference volume) $m$. We then assume that the possible motions of the system are those such that the variation of the action equals the opposite of the time integral of the virtual work of the dissipative forces corresponding to the variation of the fields. Solving the variational problem we get the equations of motion.

The explicit form of the equations of motion depends not only on the form of the action but also on the way in which the dissipation phenomenon is modeled [Nield 2007, Lopatnikov et al. 2004]. Assuming the potential energy to be quadratic in the first derivatives of the strain and of the fluid density, the evolution is described by an Allen-Cahn-like equation for the field $m$ provided that dissipative forces are assumed to be proportional to the second derivative of the velocity of the fluid. We refer to this assumption as to the Stokes dissipation. On the other hand, we find that the evolution is described by a Cahn-Hilliard-like equation for the field $m$ provided dissipation is modeled via the Darcy law, namely, the dissipative forces are proportional to the velocity of the fluid.

It is possible to give a nice physical interpretation to the fact that the evolution is described by the Allen-Cahn or the Cahn-Hilliard-like equations depending on the way in which dissipation is put into the game. With homogenization arguments it has been proven that when the characteristic length at the local ($\ell$) and at the macroscopic ($L$) scales separate, say $\ell/L \ll 1$, then the resistance experienced by the fluid when flowing through the porous material can be described by the Darcy law [Ene et al. 1975]. On the other hand when the separation of scales is poor, the Darcy limit unfairly approximates the flow, since the macroscopic characteristic length is of the same order of magnitude as the pore characteristic size [Auriault et al. 2005]. In other words in this case the solid grains have to be considered as microporous themselves.

Moreover, we note that the steady Stokes equations, which describe the behavior of the fluid at the microscopic level, can be upscaled to a set of equations describing a non-Stokesian flow only if the size of the solid obstacles to the fluid motion have a critical size with respect to the inter-obstacle distance, see for more details [Allaire 1997, Theorem 3.1]. This condition provide the so-called Brinkman limit. Conversely when the ostables are too small the homogenized equations coincide with the Stokes equations. Thus the Allen-Cahn-like equations, which correspond to a purely Stokesian flow, are controlled by the viscosity of the fluid; on the other hand the Cahn-Hilliard-like equations, which correspond to a Darcy flow, are controlled by the permeability of the solid skeleton. In particular the Cahn-Hilliard-like model reads, within the framework of soil consolidation theory, as a generalization to strain gradient materials of the one-dimensional Terzaghi problem, see [Madeo et al. 2008].

The paper is organized as follows. We first introduce the model in Section 2. In Section 3 we assume small variations of the fields with respect to some ref-
ference values, and we find the Allen-Cahn and Cahn-Hilliard-like equations. In Section 4 we consider a special model allowing to describe solid-fluid segregation in consolidation [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] and we give a numerical study of the equations deduced in the previous sections.

## 2. The model

In this section we introduce the one dimensional poromechanical model whose geometrically linearized version will be studied in the following sections. Kinematics will be briefly resumed starting from the general statement of the model [Coussy 2004] together with some particular issue introduced in [Sciarra et al. 2008]. The equations governing the behavior of the porous system will then be deduced prescribing the conservative part of the constitutive law through a suitable potential energy density $\Phi$ and the dissipative contributions through purely Stokes or Darcy terms. Special emphasis will be given to the boundary conditions and the extended definition of the essential and the natural ones will be discussed.

### 2.1. Poromechanics setup

Let $B_s := [\ell_1, \ell_2] \subset \mathbb{R}$, with $\ell_1, \ell_2 \in \mathbb{R}$, and $B_f := \mathbb{R}$ be the reference configurations for the solid and fluid components [Coussy 2004]. The solid placement $\chi_s : B_s \times \mathbb{R} \to \mathbb{R}$ is a $C^2$ function such that the map $\chi_s(\cdot, t)$, associating to each $X_s \in B_s$ the position occupied at time $t$ by the particle labeled by $X_s$ in the reference configuration $B_s$, is a $C^2$-diffeomorphism. The fluid placement map $\chi_f : B_f \times \mathbb{R} \to \mathbb{R}$ is defined analogously. The current configuration $B_t := \chi_s(B_s, t)$ at time $t$ is the set of positions of the superposed solid and fluid particles.

Consider the $C^2$ function $\phi : B_s \times \mathbb{R} \to B_f$ such that $\phi(X_s, t)$ is the fluid particle that at time $t$ occupies the same position of the solid particle $X_s$; assume, also, that $\phi(\cdot, t)$ is a $C^2$-diffeomorphism mapping univocally a solid particle into a fluid one. The three fields $\chi_s, \chi_f$, and $\phi$ are not at all independent; indeed, by definition, we immediately have that $\chi_f(\phi(X_s, t), t) = \chi_s(X_s, t)$ for any $X_s \in B_s$ and $t \in \mathbb{R}$.

In the sequel we shall often use the inverse functions of the space sections of the field $\chi_f, \chi_s$, and $\phi$. We shall misuse the notation and let $\phi^{-1}(\cdot, t)$ be the inverse of the map $X_s \to \phi(X_s, t)$ at a given time $t$. Similarly we shall also consider $\chi_s^{-1}(\cdot, t)$ and $\chi_f^{-1}(\cdot, t)$.

The Lagrangian velocities are two maps associating with each time and each point in the solid and fluid reference space the velocities of the corresponding solid and fluid particles at the specified time. More precisely, the Lagrangian velocities are the two maps $u_\alpha : B_\alpha \times \mathbb{R} \to \mathbb{R}$ defined by setting

$$u_{\alpha}(X_{\alpha}, t) := \frac{\partial \chi_{\alpha}}{\partial t}(X_{\alpha}, t) \qquad (4)$$

for any $X_\alpha \in B_\alpha$, where $\alpha = \text{s,f}$. We also consider the Eulerian velocities $v_\alpha : B_t \times \mathbb{R} \to \mathbb{R}$ associating with each point $x \in B_t$ and for each time $t \in \mathbb{R}$ the velocities of the solid and fluid particle occupying the place $x$ at time $t$; more precisely we set $v_\alpha(x, t) := u_\alpha(\chi_\alpha^{-1}(x, t), t)$.

In studying the dynamics of the porous system one can arbitrarily choose two among the three fields $\chi_s, \chi_f$, and $\phi$. Since the reference configuration $B_s$
ference values, and we find the Allen-Cahn and Cahn-Hilliard-like equations. In Section 4 we consider a special model allowing to describe solid-fluid segregation in consolidation [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] and we give a numerical study of the equations deduced in the previous sections.

## 2. The model

In this section we introduce the one dimensional poromechanical model whose geometrically linearized version will be studied in the following sections. Kinematics will be briefly resumed starting from the general statement of the model [Coussy 2004] together with some particular issue introduced in [Sciarra et al. 2008]. The equations governing the behavior of the porous system will then be deduced prescribing the conservative part of the constitutive law through a suitable potential energy density $\Phi$ and the dissipative contributions through purely Stokes or Darcy terms. Special emphasis will be given to the boundary conditions and the extended definition of the essential and the natural ones will be discussed.

### 2.1. Poromechanics setup

Let $B_s := [\ell_1, \ell_2] \subset \mathbb{R}$, with $\ell_1, \ell_2 \in \mathbb{R}$, and $B_f := \mathbb{R}$ be the reference configurations for the solid and fluid components [Coussy 2004]. The solid placement $\chi_s : B_s \times \mathbb{R} \to \mathbb{R}$ is a $C^2$ function such that the map $\chi_s(\cdot, t)$, associating to each $X_s \in B_s$ the position occupied at time $t$ by the particle labeled by $X_s$ in the reference configuration $B_s$, is a $C^2$-diffeomorphism. The fluid placement map $\chi_f : B_f \times \mathbb{R} \to \mathbb{R}$ is defined analogously. The current configuration $B_t := \chi_s(B_s, t)$ at time $t$ is the set of positions of the superposed solid and fluid particles.

Consider the $C^2$ function $\phi : B_s \times \mathbb{R} \to B_f$ such that $\phi(X_s, t)$ is the fluid particle that at time $t$ occupies the same position of the solid particle $X_s$; assume, also, that $\phi(\cdot, t)$ is a $C^2$-diffeomorphism mapping univocally a solid particle into a fluid one. The three fields $\chi_s, \chi_f$, and $\phi$ are not at all independent; indeed, by definition, we immediately have that $\chi_f(\phi(X_s, t), t) = \chi_s(X_s, t)$ for any $X_s \in B_s$ and $t \in \mathbb{R}$.

In the sequel we shall often use the inverse functions of the space sections of the field $\chi_f, \chi_s$, and $\phi$. We shall misuse the notation and let $\phi^{-1}(\cdot, t)$ be the inverse of the map $X_s \to \phi(X_s, t)$ at a given time $t$. Similarly we shall also consider $\chi_s^{-1}(\cdot, t)$ and $\chi_f^{-1}(\cdot, t)$.

The Lagrangian velocities are two maps associating with each time and each point in the solid and fluid reference space the velocities of the corresponding solid and fluid particles at the specified time. More precisely, the Lagrangian velocities are the two maps $u_\alpha : B_\alpha \times \mathbb{R} \to \mathbb{R}$ defined by setting

$$u_{\alpha}(X_{\alpha}, t) := \frac{\partial \chi_{\alpha}}{\partial t}(X_{\alpha}, t) \qquad (4)$$

for any $X_\alpha \in B_\alpha$, where $\alpha = \text{s,f}$. We also consider the Eulerian velocities $v_\alpha : B_t \times \mathbb{R} \to \mathbb{R}$ associating with each point $x \in B_t$ and for each time $t \in \mathbb{R}$ the velocities of the solid and fluid particle occupying the place $x$ at time $t$; more precisely we set $v_\alpha(x, t) := u_\alpha(\chi_\alpha^{-1}(x, t), t)$.

In studying the dynamics of the porous system one can arbitrarily choose two among the three fields $\chi_s, \chi_f$, and $\phi$. Since the reference configuration $B_s$
of the solid component is known a priori, a good choice appears to be that of expressing all the dynamical observables in terms of the fields $\chi_s$ and $\phi$ which are defined on $B_s$. Consider, for instance, the Lagrangian velocity $u_f$ of the fluid component which is defined on $B_f \times \mathbb{R}$; we prove that for any $X_s \in B_s$ and $t \in \mathbb{R}$

$$u_f(\phi(X_s, t), t) = \dot{\chi}_s(X_s, t) - \frac{\dot{\phi}(X_s, t)}{\phi'(X_s, t)} \chi'_s(X_s, t) \quad (5)$$

where the dot and the prime denote, here and in the sequel, the partial derivative with respect to time and to the space variable $X_s$ respectively. The above equation gives the expression of the fluid Lagrangian velocity in terms of $\chi_s$ and $\phi$.

In order to prove (5), note that $\chi_f(X_f, t) = \chi_s(\phi^{-1}(X_f, t), t)$ for any $X_f$ and $t$, where we have used the definition of fluid placement map. By the definition (4) of the Lagrangian velocity of the fluid, we then get

$$u_f(X_f, t) = \chi'_s(\phi^{-1}(X_f, t), t) \dot{\phi}^{-1}(X_f, t) + \dot{\chi}_s(\phi^{-1}(X_f, t), t)$$

Since $\phi^{-1}(\cdot, t)$ is the inverse function of $\phi(\cdot, t)$, we have that $\phi(\phi^{-1}(X_f, t), t) = X_f$ for any $X_f \in B_f$. By deriving this equality with respect to time we get

$$\phi'(\phi^{-1}(X_f, t), t) \dot{\phi}^{-1}(X_f, t) + \dot{\phi}(\phi^{-1}(X_f, t), t) = 0$$

The last two equations yield (5).

## 2.2. Variational principle

In order to write the equations of motion of the system we use a variational approach much similar to that developed in [Sciarra et al. 2008]; the differences between the two computations will be pointed out. For the sake of self-sufficiency we shall report the computation in detail for the one-dimensional case.

It is natural to assume that, if the system is acted upon only by conservative forces, its dynamics is described by a *Lagrangian density* $\mathcal{L}$, relative to the solid reference configuration space volume, depending on the space variable $X_s$ and on time through (in principle) $\dot{\chi}_s$, $\dot{\phi}$, $\chi_s''$, $\phi''$, $\chi'_s$, $\phi'$, $\chi_s$, and $\phi$. The Lagrangian density is equal to the *kinetic energy density* minus the *overall potential energy density* accounting for both the internal and the external conservative forces. The equations of motion for the two fields $\chi_s$ and $\phi$ can be derived assuming that the possible motions of the system in an interval of time $(t_1, t_2) \subset \mathbb{R}$ are those such that the fields $\chi_s$ and $\phi$ are extremals for the *action functional*

$$A(\chi_s, \phi) := \int_{t_1}^{t_2} dt \int_{B_s} dX_s \mathcal{L}(\dot{\chi}_s(X_s, t), \dots, \phi(X_s, t)) \quad (6)$$

in correspondence of the independent variations of the two fields $\chi_s$ and $\phi$ on $B_s \times (t_1, t_2)$. In other words any possible motion of the system in the considered interval is a solution of the Euler-Lagrange equations associated to the variational principle $\delta A = 0$.

A different variational principle is needed if the fluid component of the system is acted upon by dissipative forces; the virtual work made by these forces must be taken into account. Consider the independent variations $\delta\chi_s$ and $\delta\phi$ of the two fields $\chi_s$ and $\phi$ and denote by $\delta W$ the corresponding elementary *virtual work*
of the solid component is known a priori, a good choice appears to be that of expressing all the dynamical observables in terms of the fields $\chi_s$ and $\phi$ which are defined on $B_s$. Consider, for instance, the Lagrangian velocity $u_f$ of the fluid component which is defined on $B_f \times \mathbb{R}$; we prove that for any $X_s \in B_s$ and $t \in \mathbb{R}$

$$u_f(\phi(X_s, t), t) = \dot{\chi}_s(X_s, t) - \frac{\dot{\phi}(X_s, t)}{\phi'(X_s, t)} \chi'_s(X_s, t) \quad (5)$$

where the dot and the prime denote, here and in the sequel, the partial derivative with respect to time and to the space variable $X_s$ respectively. The above equation gives the expression of the fluid Lagrangian velocity in terms of $\chi_s$ and $\phi$.

In order to prove (5), note that $\chi_f(X_f, t) = \chi_s(\phi^{-1}(X_f, t), t)$ for any $X_f$ and $t$, where we have used the definition of fluid placement map. By the definition (4) of the Lagrangian velocity of the fluid, we then get

$$u_f(X_f, t) = \chi'_s(\phi^{-1}(X_f, t), t) \dot{\phi}^{-1}(X_f, t) + \dot{\chi}_s(\phi^{-1}(X_f, t), t)$$

Since $\phi^{-1}(\cdot, t)$ is the inverse function of $\phi(\cdot, t)$, we have that $\phi(\phi^{-1}(X_f, t), t) = X_f$ for any $X_f \in B_f$. By deriving this equality with respect to time we get

$$\phi'(\phi^{-1}(X_f, t), t) \dot{\phi}^{-1}(X_f, t) + \dot{\phi}(\phi^{-1}(X_f, t), t) = 0$$

The last two equations yield (5).

## 2.2. Variational principle

In order to write the equations of motion of the system we use a variational approach much similar to that developed in [Sciarra et al. 2008]; the differences between the two computations will be pointed out. For the sake of self-sufficiency we shall report the computation in detail for the one-dimensional case.

It is natural to assume that, if the system is acted upon only by conservative forces, its dynamics is described by a *Lagrangian density* $\mathcal{L}$, relative to the solid reference configuration space volume, depending on the space variable $X_s$ and on time through (in principle) $\dot{\chi}_s$, $\dot{\phi}$, $\chi_s''$, $\phi''$, $\chi'_s$, $\phi'$, $\chi_s$, and $\phi$. The Lagrangian density is equal to the *kinetic energy density* minus the *overall potential energy density* accounting for both the internal and the external conservative forces. The equations of motion for the two fields $\chi_s$ and $\phi$ can be derived assuming that the possible motions of the system in an interval of time $(t_1, t_2) \subset \mathbb{R}$ are those such that the fields $\chi_s$ and $\phi$ are extremals for the *action functional*

$$A(\chi_s, \phi) := \int_{t_1}^{t_2} dt \int_{B_s} dX_s \mathcal{L}(\dot{\chi}_s(X_s, t), \dots, \phi(X_s, t)) \quad (6)$$

in correspondence of the independent variations of the two fields $\chi_s$ and $\phi$ on $B_s \times (t_1, t_2)$. In other words any possible motion of the system in the considered interval is a solution of the Euler-Lagrange equations associated to the variational principle $\delta A = 0$.

A different variational principle is needed if the fluid component of the system is acted upon by dissipative forces; the virtual work made by these forces must be taken into account. Consider the independent variations $\delta\chi_s$ and $\delta\phi$ of the two fields $\chi_s$ and $\phi$ and denote by $\delta W$ the corresponding elementary *virtual work*
made by the dissipative forces acting on the fluid component. The possible motions of the system, see for instance [Blanchard et al. 1992, Chapter 5], in an interval of time $(t_1, t_2) \subset \mathbb{R}$ are those such that the fields $\chi_S$ and $\phi$ satisfies the variational principle

$$ \delta A = - \int_{t_1}^{t_2} \delta W dt \qquad (7) $$

namely, the variation of the the action integral in correspondence of a possible motion is equal to the integral over time of minus the virtual work of the dissipative forces corresponding to the considered variation of the fields.

Up to now the discussion has been conducted along very general lines. In the following sections we shall specify first the form of the virtual work of the drag forces and then that of the action. We shall obtain, finally, an explicit representation of the equations of motion (see (26) and (27)) with suitable boundary conditions.

## 2.3. Virtual work of the dissipation forces

The way in which dissipation has to be introduced in saturated porous media models is still under debate, see for instance [Nield 2007]. In particular according to the effectiveness of the hypothesis of separation of scales, between the local and macroscopic level, Darcy's or Stokes' effects are accounted for.

We first consider the so-called Darcy effect, i.e., the dissipation due to forces proportional to the velocity of the fluid component measured with respect to the solid. By following phenomenological arguments [Bear 1972] or by developing suitable homogenization schemes [Ene et al. 1975] it is seen that the fluid flow is, in this case, controlled by the permeability of the porous material. A natural expression for the virtual work of the dissipation forces acting on the fluid component and taking into account the Darcy effect is

$$ \delta W_D := - \int_{B_t} D[v_f(x, t) - v_S(x, t)][\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t)] dx \quad (8) $$

where $D > 0$ is a constant proportional to the inverse of permeability and $\delta\chi_f$ is the variation of the field $\chi_f$ induced by the independent variations $\delta\chi_S$ and $\delta\phi$. The quantity $-D[v_f - v_S]$ is the dissipative force density (relative to the current configuration space volume) which depends on the kinematic fields only through the velocity of the fluid component relative to the solid. Apparently equation (8) for the virtual work of Darcy dissipative forces is consistent with the classical expression of Darcy dissipation associated with the viscous flow through a porous continuum, see [Coussy 2004].

Following the recipe described above we have to express the virtual work in terms of the field $\chi_S$ and $\phi$. We first remark [Sciarra et al. 2008, Appendix C] that

$$ \delta\chi_f(\phi(X_S, t), t) - \delta\chi_S(X_S, t) = - \frac{\chi'_S(X_S, t)}{\phi'(X_S, t)} \delta\phi(X_S, t) \qquad (9) $$

for any $X_S \in B_S$ and $t \in \mathbb{R}$. Moreover, noted that $u_S(X_S, t) = \dot{\chi}_S(X_S, t)$, by using the definition of the Eulerian velocities, from (5) we get

$$ v_f(x, t) - v_S(x, t) = - \frac{\dot{\phi}(\chi_S^{-1}(x, t), t)}{\phi'(\chi_S^{-1}(x, t), t)} \chi'_S(\chi_S^{-1}(x, t), t) \qquad (10) $$
made by the dissipative forces acting on the fluid component. The possible motions of the system, see for instance [Blanchard et al. 1992, Chapter 5], in an interval of time $(t_1, t_2) \subset \mathbb{R}$ are those such that the fields $\chi_S$ and $\phi$ satisfies the variational principle

$$ \delta A = - \int_{t_1}^{t_2} \delta W dt \qquad (7) $$

namely, the variation of the the action integral in correspondence of a possible motion is equal to the integral over time of minus the virtual work of the dissipative forces corresponding to the considered variation of the fields.

Up to now the discussion has been conducted along very general lines. In the following sections we shall specify first the form of the virtual work of the drag forces and then that of the action. We shall obtain, finally, an explicit representation of the equations of motion (see (26) and (27)) with suitable boundary conditions.

## 2.3. Virtual work of the dissipation forces

The way in which dissipation has to be introduced in saturated porous media models is still under debate, see for instance [Nield 2007]. In particular according to the effectiveness of the hypothesis of separation of scales, between the local and macroscopic level, Darcy's or Stokes' effects are accounted for.

We first consider the so-called Darcy effect, i.e., the dissipation due to forces proportional to the velocity of the fluid component measured with respect to the solid. By following phenomenological arguments [Bear 1972] or by developing suitable homogenization schemes [Ene et al. 1975] it is seen that the fluid flow is, in this case, controlled by the permeability of the porous material. A natural expression for the virtual work of the dissipation forces acting on the fluid component and taking into account the Darcy effect is

$$ \delta W_D := - \int_{B_t} D[v_f(x, t) - v_S(x, t)][\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t)] dx \quad (8) $$

where $D > 0$ is a constant proportional to the inverse of permeability and $\delta\chi_f$ is the variation of the field $\chi_f$ induced by the independent variations $\delta\chi_S$ and $\delta\phi$. The quantity $-D[v_f - v_S]$ is the dissipative force density (relative to the current configuration space volume) which depends on the kinematic fields only through the velocity of the fluid component relative to the solid. Apparently equation (8) for the virtual work of Darcy dissipative forces is consistent with the classical expression of Darcy dissipation associated with the viscous flow through a porous continuum, see [Coussy 2004].

Following the recipe described above we have to express the virtual work in terms of the field $\chi_S$ and $\phi$. We first remark [Sciarra et al. 2008, Appendix C] that

$$ \delta\chi_f(\phi(X_S, t), t) - \delta\chi_S(X_S, t) = - \frac{\chi'_S(X_S, t)}{\phi'(X_S, t)} \delta\phi(X_S, t) \qquad (9) $$

for any $X_S \in B_S$ and $t \in \mathbb{R}$. Moreover, noted that $u_S(X_S, t) = \dot{\chi}_S(X_S, t)$, by using the definition of the Eulerian velocities, from (5) we get

$$ v_f(x, t) - v_S(x, t) = - \frac{\dot{\phi}(\chi_S^{-1}(x, t), t)}{\phi'(\chi_S^{-1}(x, t), t)} \chi'_S(\chi_S^{-1}(x, t), t) \qquad (10) $$
In order to get an expression of the virtual work comparable to that of the action, we have to rewrite the integral as an integral extended to the solid reference configuration. Given $t$, by performing the change of variables $x = \chi_S(X_S, t)$ in the integral in (8) and by using equations (9) and (10) we get

$$ \delta W_D = -D \int_{B_S} \frac{\dot{\phi}(X_S, t)}{\phi'(X_S, t)} \chi_S'(X_S, t) \frac{\chi_S'(X_S, t)}{\phi'(X_S, t)} \delta\phi(X_S, t) \chi_S'(X_S, t) dX_S \quad (11) $$

For the virtual work of the drag forces it is often considered also the contribution of the Stokes effect [Nield 2007], namely, the virtual work of dissipative forces controlled by the second derivative of the velocity of the fluid component relative to the solid. This contribution is typically the leading term in the case when separation of scales between the local and the macroscopic level is no more valid. In these cases dissipative stresses must be taken into account at macroscopic level too [Allaire 1997] and, thus, the second derivative of the fluid velocity appears in the governing equations. A naive choice for the virtual work of the dissipation forces would be

$$ -\int_{B_t} S[v_f(x, t) - v_S(x, t)]'' [\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t)] dx \quad (12) $$

with $S$ a positive constant; note that the second derivative is taken with respect to the space variable $x$ in the current configuration $B_t$. Following the same thermodynamic argument as that developed for Darcy dissipation, to the expression above of the virtual work corresponds the effective power

$$ -\int_{B_t} S[v_f(x, t) - v_S(x, t)]'' [v_f(x, t) - v_S(x, t)] dx \quad (13) $$

Since it is not negative definite, the proposed expression of the virtual work is not physically acceptable. This problem has been nicely solved in [Sciarra et al. 2008] where the following expression of the virtual work

$$ \delta W_S := -\int_{B_t} S[v_f(x, t) - v_S(x, t)]' [\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t)]' dx \quad (14) $$

with $S > 0$, has been proposed. Equation (14) for the virtual work of dissipative forces is indeed definitely consistent with the expression of the fluid dissipation when the fluid velocity satisfies the Navier-Stokes equations. As a matter of facts replacing the virtual displacement with velocities in the above expression yields the effective power of the Stokes drag forces which results to be negative definite. Note that, by performing an integration by part of (14), a term similar to (12) is found plus a contribution on the boundary of the actual domain.

We rewrite, now, the expression of the Stokes virtual work as we have done for the Darcy contribution and we shall find an equation similar to (11). In other words we shall express all the functions appearing in (14) in terms of the fields $\chi_S$ and $\phi$ and, via the changing of variable $x = \chi_S(X_S, t)$ for any fixed $t$, extend the integral to the domain $B_S$. The computation is performed in the Appendix A; the result is

$$ \delta W_S = -S \int_{B_S} \frac{1}{(\phi')^3} (\dot{\phi}\phi'\chi_S'' + \dot{\phi}'\phi'\chi_S' - \dot{\phi}\phi''\chi_S') [\delta\phi' + \frac{1}{\chi_S'\phi'}(\chi_S''\phi' - \phi''\chi_S')] \delta\phi] dX_S \quad (15) $$
In order to get an expression of the virtual work comparable to that of the action, we have to rewrite the integral as an integral extended to the solid reference configuration. Given $t$, by performing the change of variables $x = \chi_S(X_S, t)$ in the integral in (8) and by using equations (9) and (10) we get

$$ \delta W_D = -D \int_{B_S} \frac{\dot{\phi}(X_S, t)}{\phi'(X_S, t)} \chi_S'(X_S, t) \frac{\chi_S'(X_S, t)}{\phi'(X_S, t)} \delta\phi(X_S, t) \chi_S'(X_S, t) dX_S \quad (11) $$

For the virtual work of the drag forces it is often considered also the contribution of the Stokes effect [Nield 2007], namely, the virtual work of dissipative forces controlled by the second derivative of the velocity of the fluid component relative to the solid. This contribution is typically the leading term in the case when separation of scales between the local and the macroscopic level is no more valid. In these cases dissipative stresses must be taken into account at macroscopic level too [Allaire 1997] and, thus, the second derivative of the fluid velocity appears in the governing equations. A naive choice for the virtual work of the dissipation forces would be

$$ -\int_{B_t} S[v_f(x, t) - v_S(x, t)]'' [\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t)] dx \quad (12) $$

with $S$ a positive constant; note that the second derivative is taken with respect to the space variable $x$ in the current configuration $B_t$. Following the same thermodynamic argument as that developed for Darcy dissipation, to the expression above of the virtual work corresponds the effective power

$$ -\int_{B_t} S[v_f(x, t) - v_S(x, t)]'' [v_f(x, t) - v_S(x, t)] dx \quad (13) $$

Since it is not negative definite, the proposed expression of the virtual work is not physically acceptable. This problem has been nicely solved in [Sciarra et al. 2008] where the following expression of the virtual work

$$ \delta W_S := -\int_{B_t} S[v_f(x, t) - v_S(x, t)]' [\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t)]' dx \quad (14) $$

with $S > 0$, has been proposed. Equation (14) for the virtual work of dissipative forces is indeed definitely consistent with the expression of the fluid dissipation when the fluid velocity satisfies the Navier-Stokes equations. As a matter of facts replacing the virtual displacement with velocities in the above expression yields the effective power of the Stokes drag forces which results to be negative definite. Note that, by performing an integration by part of (14), a term similar to (12) is found plus a contribution on the boundary of the actual domain.

We rewrite, now, the expression of the Stokes virtual work as we have done for the Darcy contribution and we shall find an equation similar to (11). In other words we shall express all the functions appearing in (14) in terms of the fields $\chi_S$ and $\phi$ and, via the changing of variable $x = \chi_S(X_S, t)$ for any fixed $t$, extend the integral to the domain $B_S$. The computation is performed in the Appendix A; the result is

$$ \delta W_S = -S \int_{B_S} \frac{1}{(\phi')^3} (\dot{\phi}\phi'\chi_S'' + \dot{\phi}'\phi'\chi_S' - \dot{\phi}\phi''\chi_S') [\delta\phi' + \frac{1}{\chi_S'\phi'}(\chi_S''\phi' - \phi''\chi_S')] \delta\phi] dX_S \quad (15) $$
where we have omitted, for simplicity, to write explicitly that all the functions are evaluated in $(X_S, t)$.

From (11) and (15) we finally have the following expression for the virtual work of the drag forces obtained by taking into account both the Darcy and the Stokes effects

$$ \delta W := \delta W_D + \delta W_S = - \int_{B_S} (R\delta\phi + Q\delta\phi') \, dX_S \quad (16) $$

where we have set

$$ R := D \frac{1}{(\phi')^2} \dot{\phi}(\chi_S')^3 + S \frac{1}{\chi_S'(\phi')^4} (\dot{\phi}\phi'\chi_S'' + \dot{\phi}'\phi'\chi_S' - \dot{\phi}\phi''\chi_S')(\chi_S''\phi' - \phi''\chi_S') \quad (17) $$

and

$$ Q = S \frac{1}{(\phi')^3} (\dot{\phi}\phi'\chi_S'' + \dot{\phi}'\phi'\chi_S' - \dot{\phi}\phi''\chi_S') \quad (18) $$

Since we need an expression of the virtual work in terms of the variations $\delta\chi_S$ and $\delta\phi$ of the basic fields, we integrate (16) by parts and get

$$ \delta W = -(Q\delta\phi)_{\ell_1}^{\ell_2} - \int_{B_S} (R - Q')\delta\phi \, dX_S \quad (19) $$

As already noticed only in the Brinkman limit Darcy and Stokes dissipative effects can be accounted for together; on the other hand, in the general case, they are treated separately. This is indeed what we shall do in Section 3.

## 2.4. Variation of the action

In order to write explicitly the variation of the action we specify, now, the form of the Lagrangian density. In the sequel we shall not consider the inertial effects, so that, the Lagrangian density will be the opposite of the potential energy density associated to both the internal and external conservative forces. As it has been shown in [Sciarra et al. 2008] (see equation (18) therein) it is reasonable to assume that the potential energy density depends on the space and time variable only via two physically relevant functions: the strain of the solid and a properly normalized fluid mass density [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010, Sciarra et al. 2008].

More precisely, consider the Jacobian $\chi_S'(X_S, t)$ of the solid placement map, which measures the ratio between current and reference volumes of the solid component, and let

$$ \varepsilon(X_S, t) := [(\chi_S'(X_S, t))^2 - 1]/2 \quad (20) $$

be the strain field. Let $\varrho_{0,f}: B_f \to \mathbb{R}$ be the fluid reference density; we define the fluid mass density field

$$ m_f(X_S, t) := \varrho_{0,f}(\phi(X_S, t))\phi'(X_S, t) \quad (21) $$

Assuming that the mass is conserved, it is proven [Sciarra et al. 2008] that the field $m_f$ can be interpreted as the fluid mass density measured with respect to the solid reference volume.

We assume the total potential energy density $\Phi$ to depend on the fields $m_f$ and $\varepsilon$ and on their space derivative $m'_f$ and $\varepsilon'$. Since $m_f = \varrho_{0,f}(\phi)\phi'$, $m'_f =$
where we have omitted, for simplicity, to write explicitly that all the functions are evaluated in $(X_S, t)$.

From (11) and (15) we finally have the following expression for the virtual work of the drag forces obtained by taking into account both the Darcy and the Stokes effects

$$ \delta W := \delta W_D + \delta W_S = - \int_{B_S} (R\delta\phi + Q\delta\phi') \, dX_S \quad (16) $$

where we have set

$$ R := D \frac{1}{(\phi')^2} \dot{\phi}(\chi_S')^3 + S \frac{1}{\chi_S'(\phi')^4} (\dot{\phi}\phi'\chi_S'' + \dot{\phi}'\phi'\chi_S' - \dot{\phi}\phi''\chi_S') (\chi_S''\phi' - \phi''\chi_S') \quad (17) $$

and

$$ Q = S \frac{1}{(\phi')^3} (\dot{\phi}\phi'\chi_S'' + \dot{\phi}'\phi'\chi_S' - \dot{\phi}\phi''\chi_S') \quad (18) $$

Since we need an expression of the virtual work in terms of the variations $\delta\chi_S$ and $\delta\phi$ of the basic fields, we integrate (16) by parts and get

$$ \delta W = -(Q\delta\phi)_{\ell_1}^{\ell_2} - \int_{B_S} (R - Q')\delta\phi \, dX_S \quad (19) $$

As already noticed only in the Brinkman limit Darcy and Stokes dissipative effects can be accounted for together; on the other hand, in the general case, they are treated separately. This is indeed what we shall do in Section 3.

## 2.4. Variation of the action

In order to write explicitly the variation of the action we specify, now, the form of the Lagrangian density. In the sequel we shall not consider the inertial effects, so that, the Lagrangian density will be the opposite of the potential energy density associated to both the internal and external conservative forces. As it has been shown in [Sciarra et al. 2008] (see equation (18) therein) it is reasonable to assume that the potential energy density depends on the space and time variable only via two physically relevant functions: the strain of the solid and a properly normalized fluid mass density [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010, Sciarra et al. 2008].

More precisely, consider the Jacobian $\chi_S'(X_S, t)$ of the solid placement map, which measures the ratio between current and reference volumes of the solid component, and let

$$ \varepsilon(X_S, t) := [(\chi_S'(X_S, t))^2 - 1]/2 \quad (20) $$

be the strain field. Let $\varrho_{0,f} : B_f \to \mathbb{R}$ be the fluid reference density; we define the fluid mass density field

$$ m_f(X_S, t) := \varrho_{0,f}(\phi(X_S, t))\phi'(X_S, t) \quad (21) $$

Assuming that the mass is conserved, it is proven [Sciarra et al. 2008] that the field $m_f$ can be interpreted as the fluid mass density measured with respect to the solid reference volume.

We assume the total potential energy density $\Phi$ to depend on the fields $m_f$ and $\varepsilon$ and on their space derivative $m'_f$ and $\varepsilon'$. Since $m_f = \varrho_{0,f}(\phi)\phi'$, $m'_f =$
$\varrho'_{0,f}(\phi)(\phi')^2 + \varrho_{0,f}(\phi)\phi'', \varepsilon = ((\chi'_S)^2 - 1)/2$, and $\varepsilon' = \chi'_S\chi''_S$, we have that the Lagrangian density $\mathcal{L} = -\Phi$ depends on the space and time variables through the fields $\phi, \phi', \phi'', \chi'_S$, and $\chi''_S$. We then have

$$ \delta A = \int_{t_1}^{t_2} dt \int_{B_S} dX_S \left( \frac{\partial \mathcal{L}}{\partial \phi} \delta \phi + \frac{\partial \mathcal{L}}{\partial \phi'} \delta \phi' + \frac{\partial \mathcal{L}}{\partial \phi''} \delta \phi'' + \frac{\partial \mathcal{L}}{\partial \chi'_S} \delta \chi'_S + \frac{\partial \mathcal{L}}{\partial \chi''_S} \delta \chi''_S \right) $$

By using the above expression for the Lagrangian density we have that

$$ \begin{aligned} \delta A = - \int_{t_1}^{t_2} dt \int_{B_S} dX_S \left[ \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi} \right) \delta \phi \right. \\ \left. + \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi'} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi'} \right) \delta \phi' + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \delta \phi'' \right. \\ \left. + \left( \frac{\partial \Phi}{\partial \varepsilon} \frac{\partial \varepsilon}{\partial \chi'_S} + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi'_S} \right) \delta \chi'_S + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \delta \chi''_S \right] \end{aligned} \quad (22) $$

To get rid of the terms depending on the variations of the derivatives of the basic fields $\chi_S$ and $\phi$ we integrate by parts. Via a standard computation we find

$$ \begin{aligned} \delta A = - \int_{t_1}^{t_2} dt \int_{B_S} dX_S \left\{ \left[ \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi} \right. \right. \\ \left. - \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi'} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi'} - \left( \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \right)' \right)' \right] \delta \phi \\ - \left( \frac{\partial \Phi}{\partial \varepsilon} \frac{\partial \varepsilon}{\partial \chi'_S} + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi'_S} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \right)' \right)' \delta \chi_S \right\} \\ - \int_{t_1}^{t_2} dt \left\{ \left( \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \delta \phi' \right)_{\ell_1}^{\ell_2} + \left( \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \delta \chi'_S \right)_{\ell_1}^{\ell_2} \right. \\ \left. + \left[ \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi'} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi'} - \left( \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \right)' \right) \delta \phi \right]_{\ell_1}^{\ell_2} \right. \\ \left. + \left[ \left( \frac{\partial \Phi}{\partial \varepsilon} \frac{\partial \varepsilon}{\partial \chi'_S} + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi'_S} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \right)' \right) \delta \chi_S \right]_{\ell_1}^{\ell_2} \right\} \end{aligned} \quad (23) $$

## 2.5. Equations of motion

By using the variational principle (7), the equation (16) for the virtual work of the drag forces, and the equation (23) for the variation of the action, we get the equations of motion

$$ \left[ \chi'_S \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \varrho_{0,f}(\phi) \left[ \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m'_f} \right)' \right]' = R - Q' \quad (24) $$

and the boundary conditions

$$ \left\{ \begin{aligned} & \frac{\partial \Phi}{\partial m'_f} \varrho_{0,f}(\phi) \delta \phi' + \left[ \left( \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m'_f} \right)' \right) \varrho_{0,f}(\phi) + \frac{\partial \Phi}{\partial m'_f} \varrho'_{0,f}(\phi) \phi' + Q \right] \delta \phi \\ & + \frac{\partial \Phi}{\partial \varepsilon} \chi'_S \delta \chi'_S + \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \chi'_S \delta \chi_S \end{aligned} \right\}_{\ell_1}^{\ell_2} = 0 \quad (25) $$
$\varrho'_{0,f}(\phi)(\phi')^2 + \varrho_{0,f}(\phi)\phi'', \varepsilon = ((\chi'_S)^2 - 1)/2$, and $\varepsilon' = \chi'_S\chi''_S$, we have that the Lagrangian density $\mathcal{L} = -\Phi$ depends on the space and time variables through the fields $\phi, \phi', \phi'', \chi'_S$, and $\chi''_S$. We then have

$$ \delta A = \int_{t_1}^{t_2} dt \int_{B_S} dX_S \left( \frac{\partial \mathcal{L}}{\partial \phi} \delta \phi + \frac{\partial \mathcal{L}}{\partial \phi'} \delta \phi' + \frac{\partial \mathcal{L}}{\partial \phi''} \delta \phi'' + \frac{\partial \mathcal{L}}{\partial \chi'_S} \delta \chi'_S + \frac{\partial \mathcal{L}}{\partial \chi''_S} \delta \chi''_S \right) $$

By using the above expression for the Lagrangian density we have that

$$ \begin{aligned} \delta A = - \int_{t_1}^{t_2} dt \int_{B_S} dX_S \left[ \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi} \right) \delta \phi \right. \\ \left. + \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi'} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi'} \right) \delta \phi' + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \delta \phi'' \right. \\ \left. + \left( \frac{\partial \Phi}{\partial \varepsilon} \frac{\partial \varepsilon}{\partial \chi'_S} + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi'_S} \right) \delta \chi'_S + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \delta \chi''_S \right] \end{aligned} \quad (22) $$

To get rid of the terms depending on the variations of the derivatives of the basic fields $\chi_S$ and $\phi$ we integrate by parts. Via a standard computation we find

$$ \begin{aligned} \delta A = - \int_{t_1}^{t_2} dt \int_{B_S} dX_S \left\{ \left[ \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi} \right. \right. \\ \left. - \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi'} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi'} - \left( \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \right)' \right)' \right] \delta \phi \\ - \left( \frac{\partial \Phi}{\partial \varepsilon} \frac{\partial \varepsilon}{\partial \chi'_S} + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi'_S} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \right)' \right)' \delta \chi_S \right\} \\ - \int_{t_1}^{t_2} dt \left\{ \left( \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \delta \phi' \right)_{\ell_1}^{\ell_2} + \left( \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \delta \chi'_S \right)_{\ell_1}^{\ell_2} \right. \\ \left. + \left[ \left( \frac{\partial \Phi}{\partial m_f} \frac{\partial m_f}{\partial \phi'} + \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi'} - \left( \frac{\partial \Phi}{\partial m'_f} \frac{\partial m'_f}{\partial \phi''} \right)' \right) \delta \phi \right]_{\ell_1}^{\ell_2} \right. \\ \left. + \left[ \left( \frac{\partial \Phi}{\partial \varepsilon} \frac{\partial \varepsilon}{\partial \chi'_S} + \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi'_S} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \frac{\partial \varepsilon'}{\partial \chi''_S} \right)' \right) \delta \chi_S \right]_{\ell_1}^{\ell_2} \right\} \end{aligned} \quad (23) $$

## 2.5. Equations of motion

By using the variational principle (7), the equation (16) for the virtual work of the drag forces, and the equation (23) for the variation of the action, we get the equations of motion

$$ \left[ \chi'_{S} \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \varrho_{0,f}(\phi) \left[ \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m'_f} \right)' \right]' = R - Q' \quad (24) $$

and the boundary conditions

$$ \left\{ \begin{aligned} & \left( \frac{\partial \Phi}{\partial m'_f} \varrho_{0,f}(\phi) \delta \phi' + \left[ \left( \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m'_f} \right)' \right) \varrho_{0,f}(\phi) + \frac{\partial \Phi}{\partial m'_f} \varrho'_{0,f}(\phi) \phi' + Q \right] \delta \phi \right. \\ & \left. + \frac{\partial \Phi}{\partial \varepsilon} \chi'_{S} \delta \chi'_{S} + \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \chi'_{S} \delta \chi_{S} \right\}_{\ell_1}^{\ell_2} = 0 \end{aligned} \right. \quad (25) $$
We remark that, though they are partially written in terms of the fields $m_f$ and $\varepsilon$, the equations (24) are evolution equations for the two fields $\chi_s$ and $\phi$. Only under a brute approximation, see the geometrical linearization discussed in Section 3, they will reduce to a set of evolutionary equations for the fields $m_f$ and $\varepsilon$.

Notice that those equations are not simply the one dimensional particularization of the equations (30)–(32) in the paper [Sciarra et al. 2008], since there the reference mass density $\rho_{0,f}$ was assumed to be constant. When infinite systems are considered not constant reference densities can be useful in order to have finite total fluid mass integral. In the sequel we shall always consider finite systems so that it will be natural to assume the fluid reference density $\rho_{0,f} : B_f \to \mathbb{R}$ to be constant. Under this hypothesis we get the following simplified expression

$$ \left[ \chi_s \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \rho_{0,f} \left[ \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m_f'} \right)' \right]' = R - Q' \quad (26) $$

for the equations of motion and

$$ \left\{ \frac{\partial \Phi}{\partial m_f'} \rho_{0,f} \delta \phi' + \left[ \left( \frac{\partial \Phi}{\partial m_f'} - \left( \frac{\partial \Phi}{\partial m_f'} \right)' \right) \rho_{0,f} + Q \right] \delta \phi + \frac{\partial \Phi}{\partial \varepsilon'} \chi_s' \delta \chi_s' + \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \chi_s' \delta \chi_s \right\}_{\ell_1}^{\ell_2} = 0 \quad (27) $$

for the boundary conditions. We stress that from now on the symbol $\rho_{0,f}$ denotes a real constant.

### 3. Allen–Cahn and Cahn–Hilliard–like equations

In this section we rewrite the equations of motion of the poromechanical system under the so called geometrical linearization assumption, namely, in the case when only small deformations are present in the system. We first introduce the displacement fields $u(X_s, t)$ and $w(X_s, t)$ by setting

$$ \chi_s(X_s, t) = X_s + u(X_s, t) \quad \text{and} \quad \phi(X_s, t) = X_s + w(X_s, t) \quad (28) $$

for any $X_s \in B_s$ and $t \in \mathbb{R}$. We then assume that $u$ and $w$ are small, together with their space and time derivatives, and write the equations of motion up to the first order in $u$, $w$, and derivatives. By using (20), (21), (17), and (18) we get

$$ m_f = \rho_{0,f}(1+w'), \quad m := m_f - \rho_{0,f} = \rho_{0,f}w', \quad \varepsilon \approx u', \quad R \approx D\dot{w}, \quad \text{and} \quad Q \approx S\dot{w}' \quad (29) $$

where $\approx$ means that all the terms of order larger than one have been neglected.

We have introduced above the field $m$. In the following we shall imagine $\Phi$ as a function of $m$ and $m'$ and the equations of motion and the boundary conditions will be written in terms of this field. This can be done easily; indeed, since $\rho_{0,f}$ is constant, we have that $\partial/\partial m_f = \partial/\partial m$ and $m'_f = m'$. From (26), (27), and (29), we get the equations of motion

$$ \left[ (1 + \varepsilon) \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \rho_{0,f} \left[ \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' \right]' = D\dot{w} - S\dot{w}'' \quad (30) $$
We remark that, though they are partially written in terms of the fields $m_f$ and $\varepsilon$, the equations (24) are evolution equations for the two fields $\chi_s$ and $\phi$. Only under a brute approximation, see the geometrical linearization discussed in Section 3, they will reduce to a set of evolutionary equations for the fields $m_f$ and $\varepsilon$.

Notice that those equations are not simply the one dimensional particularization of the equations (30)–(32) in the paper [Sciarra et al. 2008], since there the reference mass density $\varrho_{0,f}$ was assumed to be constant. When infinite systems are considered not constant reference densities can be useful in order to have finite total fluid mass integral. In the sequel we shall always consider finite systems so that it will be natural to assume the fluid reference density $\varrho_{0,f} : B_f \to \mathbb{R}$ to be constant. Under this hypothesis we get the following simplified expression

$$ \left[ \chi_s \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \varrho_{0,f} \left[ \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m_f'} \right)' \right]' = R - Q' \quad (26) $$

for the equations of motion and

$$ \left\{ \frac{\partial \Phi}{\partial m_f'} \varrho_{0,f} \delta \phi' + \left[ \left( \frac{\partial \Phi}{\partial m_f} - \left( \frac{\partial \Phi}{\partial m_f'} \right)' \right) \varrho_{0,f} + Q \right] \delta \phi + \frac{\partial \Phi}{\partial \varepsilon'} \chi_s' \delta \chi_s' + \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \chi_s' \delta \chi_s \right\}_{\ell_1}^{\ell_2} = 0 \quad (27) $$

for the boundary conditions. We stress that from now on the symbol $\varrho_{0,f}$ denotes a real constant.

### 3. Allen–Cahn and Cahn–Hilliard–like equations

In this section we rewrite the equations of motion of the poromechanical system under the so called geometrical linearization assumption, namely, in the case when only small deformations are present in the system. We first introduce the displacement fields $u(X_s, t)$ and $w(X_s, t)$ by setting

$$ \chi_s(X_s, t) = X_s + u(X_s, t) \quad \text{and} \quad \phi(X_s, t) = X_s + w(X_s, t) \quad (28) $$

for any $X_s \in B_s$ and $t \in \mathbb{R}$. We then assume that $u$ and $w$ are small, together with their space and time derivatives, and write the equations of motion up to the first order in $u$, $w$, and derivatives. By using (20), (21), (17), and (18) we get

$$ m_f = \varrho_{0,f}(1+w'), \quad m := m_f - \varrho_{0,f} = \varrho_{0,f}w', \quad \varepsilon \approx u', \quad R \approx D\dot{w}, \quad \text{and} \quad Q \approx S\dot{w}' \quad (29) $$

where $\approx$ means that all the terms of order larger than one have been neglected.

We have introduced above the field $m$. In the following we shall imagine $\Phi$ as a function of $m$ and $m'$ and the equations of motion and the boundary conditions will be written in terms of this field. This can be done easily; indeed, since $\varrho_{0,f}$ is constant, we have that $\partial/\partial m_f = \partial/\partial m$ and $m_f' = m'$. From (26), (27), and (29), we get the equations of motion

$$ \left[ (1 + \varepsilon) \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \varrho_{0,f} \left[ \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' \right]' = D\dot{w} - S\dot{w}'' \quad (30) $$
and the associated boundary conditions

$$ \left\{ \left( \frac{\partial \Phi}{\partial \varepsilon'} \right) \delta \varepsilon + \frac{\partial \Phi}{\partial m'} \delta m + \left[ \left( \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' \right) \varrho_{0,f} + \frac{S}{\varrho_{0,f}} \dot{m} \right] \delta w + \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \delta u \right\}_{\ell_1}^{\ell_2} = 0 \quad (31) $$

A set of boundary conditions implying that (31) are satisfied is

$$ \left( \frac{\partial \Phi}{\partial \varepsilon'} \delta \varepsilon + \frac{\partial \Phi}{\partial m'} \delta m \right)_{\ell_1, \ell_2} = \left[ \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' + \frac{S}{\varrho_{0,f}^2} \dot{m} \right]_{\ell_1, \ell_2} = \left[ \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right]_{\ell_1, \ell_2} = 0 \quad (32) $$

where the notation above means that the functions in brackets are evaluated both in $\ell_1$ and $\ell_2$. With this choice it is possible to fix the boundary conditions directly on fields $m$ and $\varepsilon$ (and derivatives). The equations (30) are evolution equations for the original kinematic fields, but, by treating separately the Darcy and the Stokes effects, we can deduce a system of equations for the fields $m$ and $\varepsilon$.

The equations (31) are an extension to the case of second gradient elasticity of classical natural and essential boundary conditions [dell'Isola et al. 2009]. The first equation (32) is the additional boundary condition due to the presence of the gradient terms in the potential energy density $\Phi$. This equation specifies essential boundary conditions on the derivatives of the displacement fields or natural boundary conditions on the so called double forces, see [Germain 1973]. The generalized essential boundary conditions can be read as a prescription on the derivative of the independent fields $\chi_S$ and $\phi$, see equations (20) and (21); whilst the extended natural boundary conditions prescribe, on one hand, the additional forces which the solid continuum is able to balance at the boundary and, on the other, the wetting properties of the fluid which fills the pores [Seppecher 1989].

The second and the third equations (32) provide natural boundary conditions prescribing the chemical potential of the fluid and the traction exerted on the overall porous solid, respectively. These conditions extend the classical ones (which can be easily found in the literature) for a porous solid suffering internal stresses due to applied tractions and for a saturating fluid with fixed chemical potential; see [Baek et al. 2004] for more details. In this case we use the words generalized traction and chemical potential because of the additional contribution to stress ($\partial\Phi/\partial\varepsilon$) and chemical potential ($\partial\Phi/\partial m$) provided by the spatial derivatives of the corresponding hyperstress fields, say ($\partial\Phi/\partial\varepsilon'$) and ($\partial\Phi/\partial m'$).

*Pure Darcy dissipation.* We consider the system of equations (30) for $S = 0$. Recalling that $m = \varrho_{0,f}w'$, by deriving the second between the equations of motion with respect to $X_S$, we get

$$ \left[ (1 + \varepsilon) \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \varrho_{0,f}^2 \left[ \left( \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' \right) \right]'' = D \dot{m} \quad (33) $$

Those equations, together with the boundary conditions (32), are a partial differential equation problem for the two fields $m$ and $\varepsilon$. Note that, by exploiting
and the associated boundary conditions

$$ \left\{ \left( \frac{\partial \Phi}{\partial \varepsilon'} \right) \delta \varepsilon + \frac{\partial \Phi}{\partial m'} \delta m + \left[ \left( \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' \right) \varrho_{0,f} + \frac{S}{\varrho_{0,f}} \dot{m} \right] \delta w + \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \delta u \right\}_{\ell_1}^{\ell_2} = 0 \quad (31) $$

A set of boundary conditions implying that (31) are satisfied is

$$ \left( \frac{\partial \Phi}{\partial \varepsilon'} \delta \varepsilon + \frac{\partial \Phi}{\partial m'} \delta m \right)_{\ell_1, \ell_2} = \left[ \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' + \frac{S}{\varrho_{0,f}^2} \dot{m} \right]_{\ell_1, \ell_2} = \left[ \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right]_{\ell_1, \ell_2} = 0 \quad (32) $$

where the notation above means that the functions in brackets are evaluated both in $\ell_1$ and $\ell_2$. With this choice it is possible to fix the boundary conditions directly on fields $m$ and $\varepsilon$ (and derivatives). The equations (30) are evolution equations for the original kinematic fields, but, by treating separately the Darcy and the Stokes effects, we can deduce a system of equations for the fields $m$ and $\varepsilon$.

The equations (31) are an extension to the case of second gradient elasticity of classical natural and essential boundary conditions [dell'Isola et al. 2009]. The first equation (32) is the additional boundary condition due to the presence of the gradient terms in the potential energy density $\Phi$. This equation specifies essential boundary conditions on the derivatives of the displacement fields or natural boundary conditions on the so called double forces, see [Germain 1973]. The generalized essential boundary conditions can be read as a prescription on the derivative of the independent fields $\chi_S$ and $\phi$, see equations (20) and (21); whilst the extended natural boundary conditions prescribe, on one hand, the additional forces which the solid continuum is able to balance at the boundary and, on the other, the wetting properties of the fluid which fills the pores [Seppecher 1989].

The second and the third equations (32) provide natural boundary conditions prescribing the chemical potential of the fluid and the traction exerted on the overall porous solid, respectively. These conditions extend the classical ones (which can be easily found in the literature) for a porous solid suffering internal stresses due to applied tractions and for a saturating fluid with fixed chemical potential; see [Baek et al. 2004] for more details. In this case we use the words generalized traction and chemical potential because of the additional contribution to stress ($\partial\Phi/\partial\varepsilon$) and chemical potential ($\partial\Phi/\partial m$) provided by the spatial derivatives of the corresponding hyperstress fields, say ($\partial\Phi/\partial\varepsilon'$) and ($\partial\Phi/\partial m'$).

*Pure Darcy dissipation.* We consider the system of equations (30) for $S = 0$. Recalling that $m = \varrho_{0,f}w'$, by deriving the second between the equations of motion with respect to $X_S$, we get

$$ \left[ (1 + \varepsilon) \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \varrho_{0,f}^2 \left[ \left( \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' \right)' \right]'' = D \dot{m} \quad (33) $$

Those equations, together with the boundary conditions (32), are a partial differential equation problem for the two fields $m$ and $\varepsilon$. Note that, by exploiting
the third boundary condition (32) in $\ell_1$, we get the PDE problem

$$ egin{cases} \frac{\partial \Phi}{\partial \varepsilon} - \left(\frac{\partial \Phi}{\partial \varepsilon'}\right)' = 0 \quad \text{and} \quad \varrho_{0,f}^2 \left[\frac{\partial \Phi}{\partial m} - \left(\frac{\partial \Phi}{\partial m'}\right)''\right]' = D\dot{m} \\ \left(\frac{\partial \Phi}{\partial \varepsilon'}\delta\varepsilon + \frac{\partial \Phi}{\partial m'}\delta m\right)_{\ell_1, \ell_2} = \left[\frac{\partial \Phi}{\partial m} - \left(\frac{\partial \Phi}{\partial m'}\right)\right]'_{\ell_1, \ell_2} = 0 \end{cases} \quad (34) $$

Pure Stokes dissipation. We consider the system of equations (30) for $D = 0$. Recalling that $m = \varrho_{0,f}w'$, the equations of motion (30) become

$$ \left[ (1 + \varepsilon) \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \left[ \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' + \frac{S}{\varrho_{0,f}^2} \dot{m} \right]' = 0 \quad (35) $$

Those equations, together with the boundary conditions (32), are a partial differential equation problem for the two fields $m$ and $\varepsilon$. Note that, by exploiting the second and the third boundary condition (32) in $\ell_1$, we get the PDE problem

$$ egin{cases} \frac{\partial \Phi}{\partial \varepsilon} - \left(\frac{\partial \Phi}{\partial \varepsilon'}\right)' = 0 \quad \text{and} \quad \frac{\partial \Phi}{\partial m} - \left(\frac{\partial \Phi}{\partial m'}\right)' + \frac{S}{\varrho_{0,f}^2} \dot{m} = 0 \\ \left(\frac{\partial \Phi}{\partial \varepsilon'}\delta\varepsilon + \frac{\partial \Phi}{\partial m'}\delta m\right)_{\ell_1, \ell_2} = 0 \end{cases} \quad (36) $$

It is very interesting to note that, provided the potential energy density $\Phi$ depends on the space derivatives $m'$ and $\varepsilon'$ of the fields $m$ and $\varepsilon$ as a quadratic form with constant coefficients, it follows that the second between the equations of motion (34) becomes a Cahn-Hilliard-like equation with driving field still depending parametrically on $\varepsilon$, while the second between the equations of motion (36) becomes an Allen-Cahn-like equation with driving field depending parametrically on $\varepsilon$. More precisely we specialize the model we are studying by choosing the second gradient part of the dimensionless potential energy, that is we assume

$$ \Phi(m', \varepsilon', m, \varepsilon) := \frac{1}{2} [k_1(\varepsilon')^2 + 2k_2\varepsilon'm' + k_3(m')^2] + \Psi(m, \varepsilon) \quad (37) $$

with $k_1, k_3 > 0$, $k_2 \in \mathbb{R}$ such that $k_1k_3 - k_2^2 \ge 0$. These parameters provide energy penalties for the formation of interfaces; they have the physical dimensions of squared lengths and, according with the above mentioned conditions, provide a well-grounded identification of the intrinsic characteristic lengths of the one-dimensional solid and the fluid.

With this choice of the second gradient potential energy, in the pure Darcy case the PDE problem (34) becomes

$$ egin{cases} \frac{\partial \Psi}{\partial \varepsilon} - (k_1 \varepsilon'' + k_2 m'') = 0 \quad \text{and} \quad D \dot{m} = -\varrho_{0,f}^2 (k_2 \varepsilon'' + k_3 m'' - \frac{\partial \Psi}{\partial m})'' \\ \left( (k_1 \varepsilon' + k_2 m') \delta \varepsilon + (k_2 \varepsilon' + k_3 m') \delta m \right)_{\ell_1, \ell_2} = \left( k_2 \varepsilon'' + k_3 m'' - \frac{\partial \Psi}{\partial m} \right)_{\ell_1, \ell_2} = 0 \end{cases} \quad (38) $$

while in the pure Stokes case the PDE problem (36) reads

$$ egin{cases} \frac{\partial \Psi}{\partial \varepsilon} - (k_1 \varepsilon'' + k_2 m'') = 0 \quad \text{and} \quad \frac{S}{\varrho_{0,f}^2} \dot{m} = k_2 \varepsilon'' + k_3 m'' - \frac{\partial \Psi}{\partial m} \\ \left( (k_1 \varepsilon' + k_2 m') \delta \varepsilon + (k_2 \varepsilon' + k_3 m') \delta m \right)_{\ell_1, \ell_2} = 0 \end{cases} \quad (39) $$
the third boundary condition (32) in $\ell_1$, we get the PDE problem

$$ egin{cases} 
\frac{\partial \Phi}{\partial \varepsilon} - \left(\frac{\partial \Phi}{\partial \varepsilon'}\right)' = 0 \quad \text{and} \quad \varrho_{0,f}^2 \left[\frac{\partial \Phi}{\partial m} - \left(\frac{\partial \Phi}{\partial m'}\right)'\right]'' = D\dot{m} \\ 
\left(\frac{\partial \Phi}{\partial \varepsilon'}\delta\varepsilon + \frac{\partial \Phi}{\partial m'}\delta m\right)_{\ell_1, \ell_2} = \left[\frac{\partial \Phi}{\partial m} - \left(\frac{\partial \Phi}{\partial m'}\right)'\right]_{\ell_1, \ell_2} = 0 
\end{cases} 
\quad (34) $$

Pure Stokes dissipation. We consider the system of equations (30) for $D = 0$. Recalling that $m = \varrho_{0,f}w'$, the equations of motion (30) become

$$ \left[ (1 + \varepsilon) \left( \frac{\partial \Phi}{\partial \varepsilon} - \left( \frac{\partial \Phi}{\partial \varepsilon'} \right)' \right) \right]' = 0 \quad \text{and} \quad \left[ \frac{\partial \Phi}{\partial m} - \left( \frac{\partial \Phi}{\partial m'} \right)' + \frac{S}{\varrho_{0,f}^2} \dot{m} \right]' = 0 \quad (35) $$

Those equations, together with the boundary conditions (32), are a partial differential equation problem for the two fields $m$ and $\varepsilon$. Note that, by exploiting the second and the third boundary condition (32) in $\ell_1$, we get the PDE problem

$$ egin{cases} 
\frac{\partial \Phi}{\partial \varepsilon} - \left(\frac{\partial \Phi}{\partial \varepsilon'}\right)' = 0 \quad \text{and} \quad \frac{\partial \Phi}{\partial m} - \left(\frac{\partial \Phi}{\partial m'}\right)' + \frac{S}{\varrho_{0,f}^2} \dot{m} = 0 \\ 
\left(\frac{\partial \Phi}{\partial \varepsilon'}\delta\varepsilon + \frac{\partial \Phi}{\partial m'}\delta m\right)_{\ell_1, \ell_2} = 0 
\end{cases} 
\quad (36) $$

It is very interesting to note that, provided the potential energy density $\Phi$ depends on the space derivatives $m'$ and $\varepsilon'$ of the fields $m$ and $\varepsilon$ as a quadratic form with constant coefficients, it follows that the second between the equations of motion (34) becomes a Cahn-Hilliard-like equation with driving field still depending parametrically on $\varepsilon$, while the second between the equations of motion (36) becomes an Allen-Cahn-like equation with driving field depending parametrically on $\varepsilon$. More precisely we specialize the model we are studying by choosing the second gradient part of the dimensionless potential energy, that is we assume

$$ \Phi(m', \varepsilon', m, \varepsilon) := \frac{1}{2} [k_1(\varepsilon')^2 + 2k_2\varepsilon'm' + k_3(m')^2] + \Psi(m, \varepsilon) \quad (37) $$

with $k_1, k_3 > 0$, $k_2 \in \mathbb{R}$ such that $k_1k_3 - k_2^2 \ge 0$. These parameters provide energy penalties for the formation of interfaces; they have the physical dimensions of squared lengths and, according with the above mentioned conditions, provide a well-grounded identification of the intrinsic characteristic lengths of the one-dimensional solid and the fluid.

With this choice of the second gradient potential energy, in the pure Darcy case the PDE problem (34) becomes

$$ egin{cases} 
\frac{\partial \Psi}{\partial \varepsilon} - (k_1\varepsilon'' + k_2m'') = 0 \quad \text{and} \quad D\dot{m} = -\varrho_{0,f}^2(k_2\varepsilon'' + k_3m'' - \frac{\partial \Psi}{\partial m})'' \\ 
\left( (k_1\varepsilon' + k_2m')\delta\varepsilon + (k_2\varepsilon' + k_3m')\delta m \right)_{\ell_1, \ell_2} = \left( k_2\varepsilon'' + k_3m'' - \frac{\partial \Psi}{\partial m} \right)_{\ell_1, \ell_2} = 0 
\end{cases} 
\quad (38) $$

while in the pure Stokes case the PDE problem (36) reads

$$ egin{cases} 
\frac{\partial \Psi}{\partial \varepsilon} - (k_1\varepsilon'' + k_2m'') = 0 \quad \text{and} \quad \frac{S}{\varrho_{0,f}^2} \dot{m} = k_2\varepsilon'' + k_3m'' - \frac{\partial \Psi}{\partial m} \\ 
\left( (k_1\varepsilon' + k_2m')\delta\varepsilon + (k_2\varepsilon' + k_3m')\delta m \right)_{\ell_1, \ell_2} = 0 
\end{cases} 
\quad (39) $$
We discuss, now, some general properties of the PDE problems (38) and (39). First of all we recall that an harmonic function of a single real variable equal to zero at the extremes of an interval is necessarily equal to zero in the whole interval. We then have that for both the two problems (38) and (39) the stationary solutions are the solutions of the stationary problem

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \Psi}{\partial \varepsilon} - (k_1 \varepsilon'' + k_2 m'') = 0 \quad \text{and} \quad k_2 \varepsilon'' + k_3 m'' - \frac{\partial \Psi}{\partial m} = 0 \\[1.5em] \displaystyle \left( (k_1 \varepsilon' + k_2 m') \delta \varepsilon + (k_2 \varepsilon' + k_3 m') \delta m \right)_{\ell_1, \ell_2} = 0 \end{array} \right. \quad (40) $$

widely studied in [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010].

The dissipative character of the physical problem we are studying reflects in the existence of the not increasing energy functional

$$ \begin{aligned} F(\chi_S, \phi) &:= \int_{\ell_1}^{\ell_2} \Phi(m', \varepsilon', m, \varepsilon) \, dX_S \\ &= \int_{\ell_1}^{\ell_2} \left[ \frac{1}{2} [k_1(\varepsilon')^2 + 2k_2\varepsilon'm' + k_3(m')^2] + \Psi(m, \varepsilon) \right] dX_S \end{aligned} \quad (41) $$

for both the systems (38) and (39). To prove this we compute the time derivative of the functional F evaluated on the fields $\chi_S(X_S, t)$ and $\phi(X_S, t)$. We first get

$$ \frac{dF}{dt} = \int_{\ell_1}^{\ell_2} \left[ (k_1\varepsilon' + k_2m')\dot{\varepsilon}' + (k_2\varepsilon' + k_3m')\dot{m}' + \frac{\partial\Psi}{\partial m}\dot{m} + \frac{\partial\Psi}{\partial\varepsilon}\dot{\varepsilon} \right] dX_S $$

and, by integrating by parts, we obtain

$$ \begin{aligned} \frac{dF}{dt} = & \left[ (k_1\varepsilon' + k_2m')\dot{\varepsilon} + (k_2\varepsilon' + k_3m')\dot{m} \right]_{\ell_1}^{\ell_2} \\ & + \int_{\ell_1}^{\ell_2} \left[ -(k_1\varepsilon'' + k_2m'')\dot{\varepsilon} - (k_2\varepsilon'' + k_3m'')\dot{m} + \frac{\partial\Psi}{\partial m}\dot{m} + \frac{\partial\Psi}{\partial\varepsilon}\dot{\varepsilon} \right] dX_S \end{aligned} \quad (42) $$

Assume, first, that the fields $m$ and $\varepsilon$ are solution of the Allen–Cahn–like evolution equations (39), i.e., the PDE problems describing the evolution of the system when Stokes dissipation is considered. Moreover, assume that the boundary conditions have been chosen on the fields or on their first space derivatives so that those in (39) are fulfilled, we have that

$$ \frac{dF}{dt} = - \frac{S}{\rho_{0,f}^2} \int_{\ell_1}^{\ell_2} \dot{m}^2 dX_S \quad (43) $$

which proves that the functional F in not increasing along the motions of the Allen–Cahn–like system (39).

Assume, now, that the fields $m$ and $\varepsilon$ are solution of the Cahn–Hilliard–like evolution equations (38), i.e., the PDE problems describing the evolution of the system when Darcy dissipation is considered. Moreover, assume that the boundary conditions have been chosen on the fields or on their first space derivatives so that those in (38) are fulfilled, from (42) we have that

$$ \frac{dF}{dt} = \frac{\rho_{0,f}^2}{D} \int_{\ell_1}^{\ell_2} \left[ -(k_2\varepsilon'' + k_3m'') + \frac{\partial\Psi}{\partial m} \right] \left[ -(k_2\varepsilon'' + k_3m'') + \frac{\partial\Psi}{\partial\varepsilon} \right]'' dX_S $$
We discuss, now, some general properties of the PDE problems (38) and (39). First of all we recall that an harmonic function of a single real variable equal to zero at the extremes of an interval is necessarily equal to zero in the whole interval. We then have that for both the two problems (38) and (39) the stationary solutions are the solutions of the stationary problem

$$ egin{cases} 
\frac{\partial \Psi}{\partial \varepsilon} - (k_1 \varepsilon'' + k_2 m'') = 0 \quad \text{and} \quad k_2 \varepsilon'' + k_3 m'' - \frac{\partial \Psi}{\partial m} = 0 \\ 
\left( (k_1 \varepsilon' + k_2 m') \delta \varepsilon + (k_2 \varepsilon' + k_3 m') \delta m \right)_{\ell_1, \ell_2} = 0 
\end{cases} 
\quad (40) $$

widely studied in [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010].

The dissipative character of the physical problem we are studying reflects in the existence of the not increasing energy functional

$$ 
\begin{aligned} 
F(\chi_S, \phi) &:= \int_{\ell_1}^{\ell_2} \Phi(m', \varepsilon', m, \varepsilon) \, dX_S \\ 
&= \int_{\ell_1}^{\ell_2} \left[ \frac{1}{2} [k_1(\varepsilon')^2 + 2k_2\varepsilon'm' + k_3(m')^2] + \Psi(m, \varepsilon) \right] dX_S 
\end{aligned} 
\quad (41) $$

for both the systems (38) and (39). To prove this we compute the time derivative of the functional F evaluated on the fields $\chi_S(X_S, t)$ and $\phi(X_S, t)$. We first get

$$ \frac{dF}{dt} = \int_{\ell_1}^{\ell_2} \left[ (k_1\varepsilon' + k_2m')\dot{\varepsilon}' + (k_2\varepsilon' + k_3m')\dot{m}' + \frac{\partial\Psi}{\partial m}\dot{m} + \frac{\partial\Psi}{\partial\varepsilon}\dot{\varepsilon} \right] dX_S $$

and, by integrating by parts, we obtain

$$ 
\begin{aligned} 
\frac{dF}{dt} = & \left[ (k_1\varepsilon' + k_2m')\dot{\varepsilon} + (k_2\varepsilon' + k_3m')\dot{m} \right]_{\ell_1}^{\ell_2} \\ 
& + \int_{\ell_1}^{\ell_2} \left[ -(k_1\varepsilon'' + k_2m'')\dot{\varepsilon} - (k_2\varepsilon'' + k_3m'')\dot{m} + \frac{\partial\Psi}{\partial m}\dot{m} + \frac{\partial\Psi}{\partial\varepsilon}\dot{\varepsilon} \right] dX_S 
\end{aligned} 
\quad (42) $$

Assume, first, that the fields $m$ and $\varepsilon$ are solution of the Allen–Cahn–like evolution equations (39), i.e., the PDE problems describing the evolution of the system when Stokes dissipation is considered. Moreover, assume that the boundary conditions have been chosen on the fields or on their first space derivatives so that those in (39) are fulfilled, we have that

$$ \frac{dF}{dt} = - \frac{S}{\rho_{0,f}^2} \int_{\ell_1}^{\ell_2} \dot{m}^2 dX_S \quad (43) $$

which proves that the functional F in not increasing along the motions of the Allen–Cahn–like system (39).

Assume, now, that the fields $m$ and $\varepsilon$ are solution of the Cahn–Hilliard–like evolution equations (38), i.e., the PDE problems describing the evolution of the system when Darcy dissipation is considered. Moreover, assume that the boundary conditions have been chosen on the fields or on their first space derivatives so that those in (38) are fulfilled, from (42) we have that

$$ \frac{dF}{dt} = \frac{\rho_{0,f}^2}{D} \int_{\ell_1}^{\ell_2} \left[ -(k_2\varepsilon'' + k_3m'') + \frac{\partial\Psi}{\partial m} \right] \left[ -(k_2\varepsilon'' + k_3m'') + \frac{\partial\Psi}{\partial\varepsilon} \right]'' dX_S $$
By integrating by parts and exploiting the second boundary condition in (38) we get

$$ \frac{dF}{dt} = -\frac{\rho_{0,f}^2}{D} \int_{\ell_1}^{\ell_2} \left\{ \left[ -(k_2\varepsilon'' + k_3m'') + \frac{\partial\Psi}{\partial m} \right]' \right\}^2 dX_s \quad (44) $$

which proves that the functional F in not increasing along the motions of the Cahn-Hilliard-like system (38).

In conclusion, in this section, we have written the equations of motions (26) and the associated boundary conditions (27) under the geometrical linearization assumption. We have proven that, provided the Darcy and the Stokes effects are treated separately, those equations reduce respectively to the Cahn-Hilliard-like (38) and to the Allen-Cahn-like (39) PDE problems. Moreover, we have remarked that the stationary solutions of those problems are given by the same system of equations (40). Finally, we have proven that the energy functional (41) does not increase along the solutions of both the two systems (38) and (39), which suggests that in both cases the motions will tend asymptotically to the stationary profiles.

## 4. Pressure driven phase transition

We apply, now, the theory developed above to the special model introduced in [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] and whose stationary behavior has been widely discussed in those papers. We consider the following expression for the total potential energy density in the perspective of describing the transition between a fluid-poor and a fluid-rich phase

$$ \Psi(m, \varepsilon) := \frac{\alpha}{12} m^2 (3m^2 - 8b\varepsilon m + 6b^2\varepsilon^2) + \Psi_B(m, \varepsilon) \quad (45) $$

where

$$ \Psi_B(m, \varepsilon) := p\varepsilon + \frac{1}{2}\varepsilon^2 + \frac{1}{2}a(m - b\varepsilon)^2 \quad (46) $$

is the Biot potential energy density [Biot 1941], $a > 0$ is the ratio between the fluid and the solid rigidity, $b > 0$ is a coupling between the fluid and the solid component, $p > 0$ is the external pressure, and $\alpha > 0$ is a material parameter responsible for the showing up of an additional equilibrium.

In the papers [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] we have studied the stationary solutions of the equations (38) and (39), that is to say we have studied the problem (40). In this section we give a brief account of those results and then we will study numerically the not stationary solutions and, in particular, discuss how the stationary ones are approached. The two cases, pure Darcy and pure Stokes dissipation, will be discussed separately.

First of all we write explicitly the stationary problem corresponding to the potential energy (45). By (40) we get

$$ \left\{ \begin{array}{l} -(2/3)\alpha bm^3 + \alpha b^2 m^2 \varepsilon + p + \varepsilon - ab(m - b\varepsilon) - k_1 \varepsilon'' - k_2 m'' = 0 \\ \alpha m^3 - 2\alpha bm^2 \varepsilon + \alpha b^2 m \varepsilon^2 + a(m - b\varepsilon) - k_2 \varepsilon'' - k_3 m'' = 0 \\ \left( (k_1 \varepsilon' + k_2 m') \delta \varepsilon + (k_2 \varepsilon' + k_3 m') \delta m \right)_{\ell_1, \ell_2} = 0 \end{array} \right. \quad (47) $$

where the last line is the boundary condition.
By integrating by parts and exploiting the second boundary condition in (38) we get

$$ \frac{dF}{dt} = - \frac{\theta_{0,f}^2}{D} \int_{\ell_1}^{\ell_2} \left\{ \left[ -(k_2\epsilon'' + k_3m'') + \frac{\partial\Psi}{\partial m} \right]' \right\}^2 dX_s \quad (44) $$

which proves that the functional F in not increasing along the motions of the Cahn-Hilliard-like system (38).

In conclusion, in this section, we have written the equations of motions (26) and the associated boundary conditions (27) under the geometrical linearization assumption. We have proven that, provided the Darcy and the Stokes effects are treated separately, those equations reduce respectively to the Cahn-Hilliard-like (38) and to the Allen-Cahn-like (39) PDE problems. Moreover, we have remarked that the stationary solutions of those problems are given by the same system of equations (40). Finally, we have proven that the energy functional (41) does not increase along the solutions of both the two systems (38) and (39), which suggests that in both cases the motions will tend asymptotically to the stationary profiles.

## 4. Pressure driven phase transition

We apply, now, the theory developed above to the special model introduced in [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] and whose stationary behavior has been widely discussed in those papers. We consider the following expression for the total potential energy density in the perspective of describing the transition between a fluid-poor and a fluid-rich phase

$$ \Psi(m, \varepsilon) := \frac{\alpha}{12} m^2 (3m^2 - 8b\varepsilon m + 6b^2\varepsilon^2) + \Psi_B(m, \varepsilon) \quad (45) $$

where

$$ \Psi_B(m, \varepsilon) := p\varepsilon + \frac{1}{2}\varepsilon^2 + \frac{1}{2}a(m - b\varepsilon)^2 \quad (46) $$

is the Biot potential energy density [Biot 1941], $a > 0$ is the ratio between the fluid and the solid rigidity, $b > 0$ is a coupling between the fluid and the solid component, $p > 0$ is the external pressure, and $\alpha > 0$ is a material parameter responsible for the showing up of an additional equilibrium.

In the papers [Cirillo et al. 2011, Cirillo et al. 2009, Cirillo et al. 2010] we have studied the stationary solutions of the equations (38) and (39), that is to say we have studied the problem (40). In this section we give a brief account of those results and then we will study numerically the not stationary solutions and, in particular, discuss how the stationary ones are approached. The two cases, pure Darcy and pure Stokes dissipation, will be discussed separately.

First of all we write explicitly the stationary problem corresponding to the potential energy (45). By (40) we get

$$ \left\{ \begin{array}{l} -(2/3)\alpha bm^3 + \alpha b^2 m^2 \varepsilon + p + \varepsilon - ab(m - b\varepsilon) - k_1 \varepsilon'' - k_2 m'' = 0 \\ \alpha m^3 - 2\alpha bm^2 \varepsilon + \alpha b^2 m \varepsilon^2 + a(m - b\varepsilon) - k_2 \varepsilon'' - k_3 m'' = 0 \\ \left( (k_1 \varepsilon' + k_2 m') \delta \varepsilon + (k_2 \varepsilon' + k_3 m') \delta m \right)_{\ell_1, \ell_2} = 0 \end{array} \right. \quad (47) $$

where the last line is the boundary condition.
### 4.1. Phases: constant stationary solutions

In [Cirillo et al. 2011, Cirillo et al. 2009] we have studied the constant solutions of (47) which are called phases of the system. We have proven that there exists a pressure $p_c$, called critical pressure, such that for any $p \in [0, p_c)$ there exists a single phase ($m_S(p), \varepsilon_S(p)$), called the standard phase, which is very similar to the usual solution of the Biot model. For $p > p_c$ a second phase ($m_f(p), \varepsilon_f(p)$), richer in fluid with respect to the standard phase and hence called fluid-rich phase, appears.

We have shown that the standard phase ($m_S(p), \varepsilon_S(p)$) is the solution of the two equations $m = b\varepsilon$ and $p = f_1(\varepsilon)$, for any $p > 0$, where $f_1(\varepsilon) := -\varepsilon - \alpha b^4 \varepsilon^3 / 3$. On the other hand the fluid-rich phase ($m_f(p), \varepsilon_f(p)$) is the solution, with the smallest value of $\varepsilon$ (recall $\varepsilon \in (-1/2, 0)$, so that the smallest value has indeed largest modulus), of the two equations $m = m_+(\varepsilon)$ and $p = f_+(\varepsilon)$, where

$$m_+(\varepsilon) = \frac{b}{2} \left[ \varepsilon + \sqrt{\varepsilon^2 - \frac{4a}{\alpha b^2}} \right]$$

and

$$f_+(\varepsilon) := -\varepsilon + ab[m_+(\varepsilon) - b\varepsilon] - \alpha b^2 \varepsilon m_+^2(\varepsilon) + \frac{2}{3} \alpha b m_+^3(\varepsilon)$$

For $\varepsilon \le -2/(b\sqrt{\alpha/a})$ the function $f_+(\varepsilon)$ is positive, diverging to $+\infty$ for $\varepsilon \to -\infty$, and has a minimum at $\varepsilon_c$ such that $f_+(\varepsilon_c) = p_c$; this explains why the fluid-rich phase is seen only for $p > p_c$.

Moreover it has been shown that for any $p > 0$ the point ($m_S(p), \varepsilon_S(p)$) is a minimum of the two variable potential energy $\Psi(m, \varepsilon)$ with $p$ fixed, while ($m_f(p), \varepsilon_f(p)$) is a minimum for $p > p_c$ and a saddle point for $p = p_c$. For more details we refer to [Cirillo et al. 2011].

### 4.2. Profiles: not constant stationary solutions

In [Cirillo et al. 2010] it has been proven that there exists a unique value $p_{c0}$ of the pressure, called coexistence pressure, such that the potential energy of the two phases is equal. More precisely, it has been proven that the equation $\Psi(m_S(p), \varepsilon_S(p)) = \Psi(m_f(p), \varepsilon_f(p))$ has the single solution $p_{c0}$.

The behavior of the system at the coexistence pressure is particularly interesting; from now on we shall always consider $p = p_{c0}$ and, for this reason, we shall drop $p$ from the notation. When the external pressure is equal to $p_{c0}$, none of the two above phases is favored and we ask if profiles connecting one phase to the other exist. More precisely, in [Cirillo et al. 2010] we have addressed the problem of the existence of a connection between the two phases, that is, a solution of the stationary problem (47) on $\mathbb{R}$ tending to the standard phase for $X_S \to -\infty$ and to the fluid-rich one for $X_S \to +\infty$. Using results in [Alikakos et al. 2008] we have proven that such a connection does exist provided $k_1 k_3 - k_2^2 > 0$.

We have also shown that, for $k_1 k_3 - k_2^2 = 0$ (degenerate case), the problem of finding a solution of the stationary problem can be reduced to the computation of a definite integral. Indeed, in such a case one performs the rotation of the Cartesian reference system

$$x := \frac{m + k\varepsilon}{\sqrt{1 + k^2}} \quad \text{and} \quad y := \frac{-km + \varepsilon}{\sqrt{1 + k^2}} \qquad (48)$$
### 4.1. Phases: constant stationary solutions

In [Cirillo et al. 2011, Cirillo et al. 2009] we have studied the constant solutions of (47) which are called phases of the system. We have proven that there exists a pressure $p_c$, called critical pressure, such that for any $p \in [0, p_c)$ there exists a single phase ($m_S(p), \varepsilon_S(p)$), called the standard phase, which is very similar to the usual solution of the Biot model. For $p > p_c$ a second phase ($m_f(p), \varepsilon_f(p)$), richer in fluid with respect to the standard phase and hence called fluid-rich phase, appears.

We have shown that the standard phase ($m_S(p), \varepsilon_S(p)$) is the solution of the two equations $m = b\varepsilon$ and $p = f_1(\varepsilon)$, for any $p > 0$, where $f_1(\varepsilon) := -\varepsilon - \alpha b^4 \varepsilon^3 / 3$. On the other hand the fluid-rich phase ($m_f(p), \varepsilon_f(p)$) is the solution, with the smallest value of $\varepsilon$ (recall $\varepsilon \in (-1/2, 0)$, so that the smallest value has indeed largest modulus), of the two equations $m = m_+(\varepsilon)$ and $p = f_+(\varepsilon)$, where

$$m_+(\varepsilon) = \frac{b}{2} \left[ \varepsilon + \sqrt{\varepsilon^2 - \frac{4a}{\alpha b^2}} \right]$$

and

$$f_+(\varepsilon) := -\varepsilon + ab[m_+(\varepsilon) - b\varepsilon] - \alpha b^2 \varepsilon m_+^2(\varepsilon) + \frac{2}{3} \alpha b m_+^3(\varepsilon)$$

For $\varepsilon \le -2/(b\sqrt{\alpha/a})$ the function $f_+(\varepsilon)$ is positive, diverging to $+\infty$ for $\varepsilon \to -\infty$, and has a minimum at $\varepsilon_c$ such that $f_+(\varepsilon_c) = p_c$; this explains why the fluid-rich phase is seen only for $p > p_c$.

Moreover it has been shown that for any $p > 0$ the point ($m_S(p), \varepsilon_S(p)$) is a minimum of the two variable potential energy $\Psi(m, \varepsilon)$ with $p$ fixed, while ($m_f(p), \varepsilon_f(p)$) is a minimum for $p > p_c$ and a saddle point for $p = p_c$. For more details we refer to [Cirillo et al. 2011].

### 4.2. Profiles: not constant stationary solutions

In [Cirillo et al. 2010] it has been proven that there exists a unique value $p_{co}$ of the pressure, called coexistence pressure, such that the potential energy of the two phases is equal. More precisely, it has been proven that the equation $\Psi(m_S(p), \varepsilon_S(p)) = \Psi(m_f(p), \varepsilon_f(p))$ has the single solution $p_{co}$.

The behavior of the system at the coexistence pressure is particularly interesting; from now on we shall always consider $p = p_{co}$ and, for this reason, we shall drop $p$ from the notation. When the external pressure is equal to $p_{co}$, none of the two above phases is favored and we ask if profiles connecting one phase to the other exist. More precisely, in [Cirillo et al. 2010] we have addressed the problem of the existence of a connection between the two phases, that is, a solution of the stationary problem (47) on $\mathbb{R}$ tending to the standard phase for $X_S \to -\infty$ and to the fluid-rich one for $X_S \to +\infty$. Using results in [Alikakos et al. 2008] we have proven that such a connection does exist provided $k_1 k_3 - k_2^2 > 0$.

We have also shown that, for $k_1 k_3 - k_2^2 = 0$ (degenerate case), the problem of finding a solution of the stationary problem can be reduced to the computation of a definite integral. Indeed, in such a case one performs the rotation of the Cartesian reference system

$$x := \frac{m + k\varepsilon}{\sqrt{1 + k^2}} \quad \text{and} \quad y := \frac{-km + \varepsilon}{\sqrt{1 + k^2}} \qquad (48)$$
in the plane $m-\epsilon$, where $k := k_2/k_3 = \pm\sqrt{k_1/k_3}$, and defines

$$V(x, y) = -\Psi(m(x, y), \epsilon(x, y)) \qquad (49)$$

Then one shows that the two fields $m(X_S)$ and $\epsilon(X_S)$ are solutions of the two equations (47) if and only if the corresponding fields $x(X_S)$ and $y(X_S)$ satisfy

$$k_3(1 + k^2)x'' = -\frac{\partial V}{\partial x}(x, y) \quad \text{and} \quad \frac{\partial V}{\partial y}(x, y) = 0 \qquad (50)$$

The root locus of the constraint curve $\partial V(x, y)/\partial y = 0$ is made of a certain number of maximal components such that each of them is the graph of a function $x \in \mathbb{R} \to y(x) \in \mathbb{R}$; for each of them the first between the two equations (50) becomes a standard one dimensional conservative mechanical system with $X_S$ interpreted as time, kinetic energy $(1 + k^2)k_3(x')^2/2$, and potential energy $V(x, y(x))$.

Since the function $V$ has been obtained by flipping the sign of the function $\Psi$ and rotating the coordinate axes, then at the coexistence pressure it has the two absolute maximum points $(x_S, y_S)$ and $(x_f, y_f)$ corresponding, respectively, to the standard and to the fluid-rich phases. Since $(m_S, \epsilon_S)$ and $(m_f, \epsilon_f)$ satisfy the equations $\Psi_m(m, \epsilon) = 0$ and $\Psi_\epsilon(m, \epsilon) = 0$, we have that the two points $(x_S, y_S)$ and $(x_f, y_f)$ are solutions of the constraint equation $\partial V(x, y)/\partial y = 0$ and hence they belong to the constraint curve.

In [Cirillo et al. 2010] we have seen that there exist values of the second gradient parameters $k_1$, $k_2$, and $k_3$ such that the two points above fall on the same maximal component of the constraint equation. Since, in this case, the function $V$ has two isolated absolute maximum points which, by hypothesis, belong to the same maximal component of the constraint curve, we have that the function $V(x, y(x))$ of $x$ has two absolute isolated maxima in $x_S$ and $x_f$. The motion of the equivalent one dimensional conservative mechanical system corresponding to the energy level $V_{\max} := V(x_S, y_S)$ yields the heteroclinic orbit, that is to say it yields the above mentioned connection on $\mathbb{R}$ between the standard and the fluid-rich phases.

With a similar approach one can find solutions of the system (47) with Dirichlet boundary conditions on a finite interval, say $\ell_1 = 0$ and $\ell_2 = 1$. For instance, assume $x_S < x_f$ and consider the boundary condition $x(0) = x_S$ and $x(1) = x_f$. Motions of the equivalent mechanical system started in $x_S$ at time 0 with positive velocity $x'(0) = \bar{v} > 0$ (note that if it were $x_S > x_f$ one should consider a negative velocity at time 0) will reach $x_f$ in a finite time, which is a decreasing function of the initial velocity. It is then possible to choose properly the initial velocity so that

$$1 = \int_{x_S}^{x_f} \sqrt{\frac{k_3(1 + k^2)}{2[\bar{E} - V(x, y(x))]}} dx$$

where $\bar{E} = (1+k^2)k_3\bar{v}^2/2+V_{\max}$ is the total energy of the motion of the equivalent one-dimensional conservative system. Once $\bar{E}$ has been found, the solution of the stationary problem (47) with Dirichlet boundary conditions $m(0) = m_S$, $\epsilon(0) = \epsilon_S$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ is implicitly given, in terms of the variable $x$, by the definite integral

$$X_S = \int_{x_S}^{x} \sqrt{\frac{k_3(1 + k^2)}{2[\bar{E} - V(x, y(x))]}} dx$$
in the plane $m-\varepsilon$, where $k := k_2/k_3 = \pm\sqrt{k_1/k_3}$, and defines

$$V(x, y) = -\Psi(m(x, y), \varepsilon(x, y)) \qquad (49)$$

Then one shows that the two fields $m(X_S)$ and $\varepsilon(X_S)$ are solutions of the two equations (47) if and only if the corresponding fields $x(X_S)$ and $y(X_S)$ satisfy

$$k_3(1 + k^2)x'' = -\frac{\partial V}{\partial x}(x, y) \quad \text{and} \quad \frac{\partial V}{\partial y}(x, y) = 0 \qquad (50)$$

The root locus of the constraint curve $\partial V(x, y)/\partial y = 0$ is made of a certain number of maximal components such that each of them is the graph of a function $x \in \mathbb{R} \to y(x) \in \mathbb{R}$; for each of them the first between the two equations (50) becomes a standard one dimensional conservative mechanical system with $X_S$ interpreted as time, kinetic energy $(1 + k^2)k_3(x')^2/2$, and potential energy $V(x, y(x))$.

Since the function $V$ has been obtained by flipping the sign of the function $\Psi$ and rotating the coordinate axes, then at the coexistence pressure it has the two absolute maximum points $(x_S, y_S)$ and $(x_f, y_f)$ corresponding, respectively, to the standard and to the fluid-rich phases. Since $(m_S, \varepsilon_S)$ and $(m_f, \varepsilon_f)$ satisfy the equations $\Psi_m(m, \varepsilon) = 0$ and $\Psi_\varepsilon(m, \varepsilon) = 0$, we have that the two points $(x_S, y_S)$ and $(x_f, y_f)$ are solutions of the constraint equation $\partial V(x, y)/\partial y = 0$ and hence they belong to the constraint curve.

In [Cirillo et al. 2010] we have seen that there exist values of the second gradient parameters $k_1$, $k_2$, and $k_3$ such that the two points above fall on the same maximal component of the constraint equation. Since, in this case, the function $V$ has two isolated absolute maximum points which, by hypothesis, belong to the same maximal component of the constraint curve, we have that the function $V(x, y(x))$ of $x$ has two absolute isolated maxima in $x_S$ and $x_f$. The motion of the equivalent one dimensional conservative mechanical system corresponding to the energy level $V_{\max} := V(x_S, y_S)$ yields the heteroclinic orbit, that is to say it yields the above mentioned connection on $\mathbb{R}$ between the standard and the fluid-rich phases.

With a similar approach one can find solutions of the system (47) with Dirichlet boundary conditions on a finite interval, say $\ell_1 = 0$ and $\ell_2 = 1$. For instance, assume $x_S < x_f$ and consider the boundary condition $x(0) = x_S$ and $x(1) = x_f$. Motions of the equivalent mechanical system started in $x_S$ at time 0 with positive velocity $x'(0) = \bar{v} > 0$ (note that if it were $x_S > x_f$ one should consider a negative velocity at time 0) will reach $x_f$ in a finite time, which is a decreasing function of the initial velocity. It is then possible to choose properly the initial velocity so that

$$1 = \int_{x_S}^{x_f} \sqrt{\frac{k_3(1 + k^2)}{2[\bar{E} - V(x, y(x))]}} dx$$

where $\bar{E} = (1+k^2)k_3\bar{v}^2/2+V_{\max}$ is the total energy of the motion of the equivalent one-dimensional conservative system. Once $\bar{E}$ has been found, the solution of the stationary problem (47) with Dirichlet boundary conditions $m(0) = m_S$, $\varepsilon(0) = \varepsilon_S$, $m(1) = m_f$, and $\varepsilon(1) = \varepsilon_f$ is implicitly given, in terms of the variable $x$, by the definite integral

$$X_S = \int_{x_S}^{x} \sqrt{\frac{k_3(1 + k^2)}{2[\bar{E} - V(x, y(x))]}} dx$$
The solutions found in such a way for $k_1 = 10^{-3}$, $k_3 = 10^{-3}$, $k_2 = \pm 10^{-3}$, $a = 0.5$, $b = 1$, $\alpha = 100$, and at coexistence pressure, have been depicted in figures 1 and 2.

As we have noted above, in the not degenerate case, i.e., $k_1 k_3 - k_2^2 > 0$, in [Cirillo et al. 2010] we could prove the existence of a connection, but we were not able to find it explicitly. Studying in more detail the behavior of the solutions of the system (47) is a difficult task since, in this case, the problem cannot be rewritten as a one dimensional conservative mechanical system. We shall then study the behavior of the solutions close to the two phases by linearizing the stationary equations, while far from the points representing the phases we shall use a finite difference numerical approach.

The first question that naturally arises looking at the solutions obtained in the degenerate case, see figure 1, is the nature of the bump. In the degenerate case the bump is due to the fact that the solution of the stationary problem, once projected onto the $m-\epsilon$ plane, has to lie on the constraint curve, so that its shape depends on that of the constraint itself. In the not degenerate case no constraint curve exists, in other words it does not exist any a priori privileged manifold on the plane $m-\epsilon$. It is then possible (in principle) that the bump is associated to solutions oscillating close to the phases.

In order to discuss the existence of such oscillatory solutions we linearize the problem close to the phases. Let $(\bar{m}, \bar{\epsilon})$ be one of the two points $(m_S, \epsilon_S)$ and $(m_f, \epsilon_f)$ representing, respectively, the standard and the fluid-rich phase. By expanding the first gradient part of the stationary equations (40) around $(\bar{m}, \bar{\epsilon})$ to the first order, we get the linearized equations

$$ \begin{cases} k_1\epsilon'' + k_2m'' - \Psi_{\epsilon m}(\bar{m}, \bar{\epsilon})(m - \bar{m}) - \Psi_{\epsilon\epsilon}(\bar{m}, \bar{\epsilon})(\epsilon - \bar{\epsilon}) = 0 \\ k_2\epsilon'' + k_3m'' - \Psi_{m\epsilon}(\bar{m}, \bar{\epsilon})(m - \bar{m}) - \Psi_{m\bar{\epsilon}}(\bar{m}, \bar{\epsilon})(\epsilon - \bar{\epsilon}) = 0 \end{cases} $$

If we let $q_1 = \epsilon - \bar{\epsilon}$ and $q_2 = m - \bar{m}$, the above equations become

$$ \mathbf{T}q'' - \mathbf{\Psi}q = 0 \qquad (51) $$

where we consider the column vector $q \in \mathbb{R}^2$ and set

$$ \mathbf{T} = \begin{pmatrix} k_1 & k_2 \\ k_2 & k_3 \end{pmatrix} \quad \text{and} \quad \mathbf{\Psi} = \begin{pmatrix} \Psi_{\epsilon\epsilon}(\bar{m}, \bar{\epsilon}) & \Psi_{\epsilon m}(\bar{m}, \bar{\epsilon}) \\ \Psi_{m\epsilon}(\bar{m}, \bar{\epsilon}) & \Psi_{m\bar{\epsilon}}(\bar{m}, \bar{\epsilon}) \end{pmatrix} $$

Note that the matrix $\mathbf{T}$ is trivially real and symmetric; moreover, in the not degenerate case $k_1 k_3 - k_2^2 > 0$ it is also positive definite. The matrix $\mathbf{\Psi}$, on the other hand, is real and, by Schwartz theorem, symmetric. Moreover, as we proved in [Cirillo et al. 2011], it is positive definite at any pressure greater than the critical one $p_c$; in particular it is positive definite at the coexistence pressure.

Looking for solutions of (51) in the form $q = c \exp\{\lambda X_S\}$, with $c \in \mathbb{R}^2$ and $\lambda \in \mathbb{C}$, we get that the constant vector $c$ has to be a solution of $(\lambda^2 \mathbf{T} - \mathbf{\Psi})c = 0$. Set $\mu = \lambda^2$; in order to have not vanishing $c$ solving the above equation, $\mu$ has to be a solution of the secular equation $\det(\mu \mathbf{T} - \mathbf{\Psi}) = 0$. Since $\mathbf{T}$ is real, symmetric, and positive definite and $\mathbf{\Psi}$ is real and symmetric, it follows that the solutions $\mu_1$ and $\mu_2$ of the secular equations are real. Moreover, since $\mathbf{\Psi}$ is positive definite we have that those solutions are positive themselves (see, for instance, [Gantmacher 1975, Section 40, Chapter 6]). Thus the secular
The solutions found in such a way for $k_1 = 10^{-3}$, $k_3 = 10^{-3}$, $k_2 = \pm 10^{-3}$, $a = 0.5$, $b = 1$, $\alpha = 100$, and at coexistence pressure, have been depicted in figures 1 and 2.

As we have noted above, in the not degenerate case, i.e., $k_1 k_3 - k_2^2 > 0$, in [Cirillo et al. 2010] we could prove the existence of a connection, but we were not able to find it explicitly. Studying in more detail the behavior of the solutions of the system (47) is a difficult task since, in this case, the problem cannot be rewritten as a one dimensional conservative mechanical system. We shall then study the behavior of the solutions close to the two phases by linearizing the stationary equations, while far from the points representing the phases we shall use a finite difference numerical approach.

The first question that naturally arises looking at the solutions obtained in the degenerate case, see figure 1, is the nature of the bump. In the degenerate case the bump is due to the fact that the solution of the stationary problem, once projected onto the $m-\epsilon$ plane, has to lie on the constraint curve, so that its shape depends on that of the constraint itself. In the not degenerate case no constraint curve exists, in other words it does not exist any a priori privileged manifold on the plane $m-\epsilon$. It is then possible (in principle) that the bump is associated to solutions oscillating close to the phases.

In order to discuss the existence of such oscillatory solutions we linearize the problem close to the phases. Let $(\bar{m}, \bar{\epsilon})$ be one of the two points $(m_S, \epsilon_S)$ and $(m_f, \epsilon_f)$ representing, respectively, the standard and the fluid-rich phase. By expanding the first gradient part of the stationary equations (40) around $(\bar{m}, \bar{\epsilon})$ to the first order, we get the linearized equations

$$ \begin{cases} k_1\epsilon'' + k_2m'' - \Psi_{\epsilon m}(\bar{m}, \bar{\epsilon})(m - \bar{m}) - \Psi_{\epsilon\epsilon}(\bar{m}, \bar{\epsilon})(\epsilon - \bar{\epsilon}) = 0 \\ k_2\epsilon'' + k_3m'' - \Psi_{m\epsilon}(\bar{m}, \bar{\epsilon})(m - \bar{m}) - \Psi_{m\bar{\epsilon}}(\bar{m}, \bar{\epsilon})(\epsilon - \bar{\epsilon}) = 0 \end{cases} $$

If we let $q_1 = \epsilon - \bar{\epsilon}$ and $q_2 = m - \bar{m}$, the above equations become

$$ Tq'' - \Psi q = 0 \qquad (51) $$

where we consider the column vector $q \in \mathbb{R}^2$ and set

$$ T = \begin{pmatrix} k_1 & k_2 \\ k_2 & k_3 \end{pmatrix} \quad \text{and} \quad \Psi = \begin{pmatrix} \Psi_{\epsilon\epsilon}(\bar{m}, \bar{\epsilon}) & \Psi_{\epsilon m}(\bar{m}, \bar{\epsilon}) \\ \Psi_{m\epsilon}(\bar{m}, \bar{\epsilon}) & \Psi_{m\bar{\epsilon}}(\bar{m}, \bar{\epsilon}) \end{pmatrix} $$

Note that the matrix $T$ is trivially real and symmetric; moreover, in the not degenerate case $k_1 k_3 - k_2^2 > 0$ it is also positive definite. The matrix $\Psi$, on the other hand, is real and, by Schwartz theorem, symmetric. Moreover, as we proved in [Cirillo et al. 2011], it is positive definite at any pressure greater than the critical one $p_c$; in particular it is positive definite at the coexistence pressure.

Looking for solutions of (51) in the form $q = c \exp\{\lambda X_S\}$, with $c \in \mathbb{R}^2$ and $\lambda \in \mathbb{C}$, we get that the constant vector $c$ has to be a solution of $(\lambda^2 T - \Psi)c = 0$. Set $\mu = \lambda^2$; in order to have not vanishing $c$ solving the above equation, $\mu$ has to be a solution of the secular equation $\det(\mu T - \Psi) = 0$. Since $T$ is real, symmetric, and positive definite and $\Psi$ is real and symmetric, it follows that the solutions $\mu_1$ and $\mu_2$ of the secular equations are real. Moreover, since $\Psi$ is positive definite we have that those solutions are positive themselves (see, for instance, [Gantmacher 1975, Section 40, Chapter 6]). Thus the secular
equations have four real pairwise opposite eigenvalues. This implies that the two fixed points are unstable, as it also follows from the existence of the connection, and that no oscillation around the fixed points is expected.

We come, finally, to the numerical study of the stationary problem (47) in the not degenerate case on a finite interval (we shall always use the interval $[0, 1]$). We use the finite difference method powered with the Newton-Raphson algorithm. The substitution rules we adopted are standard and have been reported in appendix Appendix B. We solved the stationary problem (47) with the Dirichlet boundary conditions $m(0) = m_S$, $\varepsilon(0) = \varepsilon_S$, $m(1) = m_f$, $\varepsilon(1) = \varepsilon_f$, at the coexistence pressure for $a = 0.5$, $b = 1$, and $\alpha = 100$. For the second gradient coefficient we have considered two sets; the first set is $k_1 = 10^{-3}$, $k_3 = 10^{-3}$, and $k_2 = -0.4 \times 10^{-3}$, $0.2 \times 10^{-3}$, $0.8 \times 10^{-3}$ and the related numerical results have been depicted in the bottom row of figure 1 and in figure 2. The second set is $k_1 = 10^{-2}$, $k_3 = 10^{-2}$, and $k_2 = -0.4 \times 10^{-2}$, $0.2 \times 10^{-2}$, $0.8 \times 10^{-2}$ and the related numerical results have been reported in the top row of figure 1.

On the bottom (resp. top) row of figure 1 we have depicted the $m$ and $\varepsilon$-profiles as functions of the space variable $X_S$; solid lines refer to the degenerate case $k_2 = \pm 10^{-3}$ (resp. $k_2 = \pm 10^{-2}$), while dotted lines refer to the not degenerate cases $k_2 = -0.4 \times 10^{-3}$, $0.2 \times 10^{-3}$, $0.8 \times 10^{-3}$ (resp. $k_2 = -0.4 \times 10^{-2}$, $0.2 \times 10^{-2}$, $0.8 \times 10^{-2}$). In figure 2 the same solutions have been depicted on the plane $m-\varepsilon$. The behavior of the solutions of the stationary problem in the degenerate and in the not degenerate case is similar. The existence of the bump in the $m$ and $\varepsilon$-profiles is essentially due to the shape of the manifold of the $m-\varepsilon$ plane on which the solutions lie.

It is worth noting that, even in the degenerate case, due to the asymmetry of the potential energy $V$, the position of the interface in the stationary profile depends on the value of the second gradient constants $k_1$, $k_2$, and $k_3$. However in the limit for small values of $k_1$, $k_2$, and $k_3$ it is observed that the interface tends to a definite position [Cirillo et al. 2012].

### 4.3. Evolution under pure Stokes dissipation: Allen-Cahn-like equations

The equations of motion describing the evolution of the system subjected to a pure Darcy drag force are the equations (39) which, for the special model (45), read

$$ \left\{ \begin{array}{l} -(2/3)\alpha bm^3 + \alpha b^2 m^2 \varepsilon + p + \varepsilon - ab(m - b\varepsilon) - k_1 \varepsilon'' - k_2 m'' = 0 \\ -\rho_{0,f}^2 [\alpha m^3 - 2\alpha bm^2 \varepsilon + \alpha b^2 m \varepsilon^2 + a(m - b\varepsilon) - k_2 \varepsilon'' - k_3 m''] = S \dot{m} \end{array} \right. \quad (52) $$

with the boundary conditions

$$ ((k_2 m' + k_1 \varepsilon') \delta \varepsilon + (k_3 m' + k_2 \varepsilon') \delta m)_{\ell_1, \ell_2} = 0 \quad (53) $$

We have studied the above PDE problem with the same numerical approach used for the stationary problem and discussed in Section 4.2. That is we have used a finite difference method powered with the Newton-Raphson algorithm and with the standard substitution rules reported in Appendix B.

In the figures 3-5 we have depicted the solutions of the problem (52) with Dirichlet boundary conditions $m(0) = m_S$, $\varepsilon(0) = \varepsilon_S$, $m(1) = m_f$, and $\varepsilon(1) = \varepsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. All the graphs refer to the same spatially
equations have four real pairwise opposite eigenvalues. This implies that the two fixed points are unstable, as it also follows from the existence of the connection, and that no oscillation around the fixed points is expected.

We come, finally, to the numerical study of the stationary problem (47) in the not degenerate case on a finite interval (we shall always use the interval $[0, 1]$). We use the finite difference method powered with the Newton-Raphson algorithm. The substitution rules we adopted are standard and have been reported in appendix Appendix B. We solved the stationary problem (47) with the Dirichlet boundary conditions $m(0) = m_S$, $\varepsilon(0) = \varepsilon_S$, $m(1) = m_f$, $\varepsilon(1) = \varepsilon_f$, at the coexistence pressure for $a = 0.5$, $b = 1$, and $\alpha = 100$. For the second gradient coefficient we have considered two sets; the first set is $k_1 = 10^{-3}$, $k_3 = 10^{-3}$, and $k_2 = -0.4 \times 10^{-3}$, $0.2 \times 10^{-3}$, $0.8 \times 10^{-3}$ and the related numerical results have been depicted in the bottom row of figure 1 and in figure 2. The second set is $k_1 = 10^{-2}$, $k_3 = 10^{-2}$, and $k_2 = -0.4 \times 10^{-2}$, $0.2 \times 10^{-2}$, $0.8 \times 10^{-2}$ and the related numerical results have been reported in the top row of figure 1.

On the bottom (resp. top) row of figure 1 we have depicted the $m$ and $\varepsilon$-profiles as functions of the space variable $X_S$; solid lines refer to the degenerate case $k_2 = \pm 10^{-3}$ (resp. $k_2 = \pm 10^{-2}$), while dotted lines refer to the not degenerate cases $k_2 = -0.4 \times 10^{-3}$, $0.2 \times 10^{-3}$, $0.8 \times 10^{-3}$ (resp. $k_2 = -0.4 \times 10^{-2}$, $0.2 \times 10^{-2}$, $0.8 \times 10^{-2}$). In figure 2 the same solutions have been depicted on the plane $m-\varepsilon$. The behavior of the solutions of the stationary problem in the degenerate and in the not degenerate case is similar. The existence of the bump in the $m$ and $\varepsilon$-profiles is essentially due to the shape of the manifold of the $m-\varepsilon$ plane on which the solutions lie.

It is worth noting that, even in the degenerate case, due to the asymmetry of the potential energy $V$, the position of the interface in the stationary profile depends on the value of the second gradient constants $k_1$, $k_2$, and $k_3$. However in the limit for small values of $k_1$, $k_2$, and $k_3$ it is observed that the interface tends to a definite position [Cirillo et al. 2012].

### 4.3. Evolution under pure Stokes dissipation: Allen-Cahn-like equations

The equations of motion describing the evolution of the system subjected to a pure Darcy drag force are the equations (39) which, for the special model (45), read

$$ \left\{ \begin{array}{l} -(2/3)\alpha bm^3 + \alpha b^2 m^2 \varepsilon + p + \varepsilon - ab(m - b\varepsilon) - k_1 \varepsilon'' - k_2 m'' = 0 \\ -\rho_{0,f}^2 [\alpha m^3 - 2\alpha bm^2 \varepsilon + \alpha b^2 m \varepsilon^2 + a(m - b\varepsilon) - k_2 \varepsilon'' - k_3 m''] = S \dot{m} \end{array} \right. \quad (52) $$

with the boundary conditions

$$ ((k_2 m' + k_1 \varepsilon') \delta \varepsilon + (k_3 m' + k_2 \varepsilon') \delta m)_{\ell_1, \ell_2} = 0 \quad (53) $$

We have studied the above PDE problem with the same numerical approach used for the stationary problem and discussed in Section 4.2. That is we have used a finite difference method powered with the Newton-Raphson algorithm and with the standard substitution rules reported in Appendix B.

In the figures 3-5 we have depicted the solutions of the problem (52) with Dirichlet boundary conditions $m(0) = m_S$, $\varepsilon(0) = \varepsilon_S$, $m(1) = m_f$, and $\varepsilon(1) = \varepsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. All the graphs refer to the same spatially
random initial condition on the interval $[0, 1]$. Details on the depicted times can be found in the caption of the related figures.

Figures 3 and 4 show how the limiting stationary profile is formed. At times of order one a deviation from the initial random profile is noticed and a bumped profile with roughness due to the randomness of the initial condition is formed. The profile is smeared out in a time of order $10^1$ and an interface, far from the position of the stationary profile, is seen. The interface then undergoes a metastable, purely dissipative, evolution approaching the stationary state in times greater than $10^4$. This coarsening process has been widely studied for the standard Allen-Cahn equation, see for instance [Carr et al. 1989, Fusco et al. 1989, Ward 1994].

Figure 5 shows that the solution is soon attracted by the manifold in the plane $m-\epsilon$ defined by the stationary solution. Afterwards the projection of the profiles onto the $m-\epsilon$ plane keeps close to such a curve. We recall that, when the parameters $k_1, k_2$, and $k_3$, as in the cases illustrated in the figures 3–5, are chosen so that $k_1k_3 - k_2^2 = 0$, the algebraic equation of the limiting manifold is the second above the two equations (50).

We have measured the way in which the stationary solution is approached considering the difference between the value of the energy functional (41) at time $t$ and that of the corresponding stationary solution. Results are reported in figure 6. We note that the energy functional is a decreasing function of time as proved in Section 3 (see equation (43)).

Some of the features described above are peculiar of the initial random condition. We have also studied the same Allen-Cahn problem, with the same parameters and the same boundary conditions, with a deterministic stratified initial state. Results for the $\epsilon$-profile are depicted at times in the figure 7; we have not reported the $m$-profiles and the $m-\epsilon$-plane view since they do not add any more information.

At times of order $10^{-1}$ a deviation from the initial stratified state is detected. The profile deforms into a new one with a droplet-like shape, at times of order $10^1$, finally an interface profile nucleates at times of order $10^2$, as in the previous case corresponding to spatially random initial data. Similarly to this case the stationary interface is approached with the same velocity. Even in this case we have checked the way in which the solution tends to the stationary state via the energy functional (41). Results are reported in figure 8.

### 4.4. Evolution under pure Darcy dissipation: Cahn-Hilliard-like equations

The equations of motion describing the evolution of the system subjected to a pure Darcy drag force are the equations (38) which, for the special model (45), read

$$ \begin{cases} -(2/3)\alpha bm^3 + \alpha b^2 m^2 \epsilon + p + \epsilon - ab(m - b\epsilon) - k_1 \epsilon'' - k_2 m'' = 0 \\ \varrho_{0,f}^2 (\alpha m^3 - 2\alpha bm^2 \epsilon + \alpha b^2 m \epsilon^2 + a(m - b\epsilon) - k_2 \epsilon'' - k_3 m'')'' = D\dot{m} \end{cases} \quad (54) $$

with the boundary conditions

$$ \begin{cases} ((k_2 m' + k_1 \epsilon') \delta \epsilon + (k_3 m' + k_2 \epsilon') \delta m)_{\ell_1, \ell_2} = 0 \\ (\alpha m^3 - 2\alpha bm^2 \epsilon + \alpha b^2 m \epsilon^2 + a(m - b\epsilon) - k_2 \epsilon'' - k_3 m'')_{\ell_1, \ell_2} = 0 \end{cases} \quad (55) $$

We have studied the above PDE problem with the same numerical approach used for the stationary problem and for the Allen-Cahn-like system. We have
random initial condition on the interval $[0, 1]$. Details on the depicted times can be found in the caption of the related figures.

Figures 3 and 4 show how the limiting stationary profile is formed. At times of order one a deviation from the initial random profile is noticed and a bumped profile with roughness due to the randomness of the initial condition is formed. The profile is smeared out in a time of order $10^1$ and an interface, far from the position of the stationary profile, is seen. The interface then undergoes a metastable, purely dissipative, evolution approaching the stationary state in times greater than $10^4$. This coarsening process has been widely studied for the standard Allen-Cahn equation, see for instance [Carr et al. 1989, Fusco et al. 1989, Ward 1994].

Figure 5 shows that the solution is soon attracted by the manifold in the plane $m-\epsilon$ defined by the stationary solution. Afterwards the projection of the profiles onto the $m-\epsilon$ plane keeps close to such a curve. We recall that, when the parameters $k_1, k_2$, and $k_3$, as in the cases illustrated in the figures 3–5, are chosen so that $k_1k_3 - k_2^2 = 0$, the algebraic equation of the limiting manifold is the second above the two equations (50).

We have measured the way in which the stationary solution is approached considering the difference between the value of the energy functional (41) at time $t$ and that of the corresponding stationary solution. Results are reported in figure 6. We note that the energy functional is a decreasing function of time as proved in Section 3 (see equation (43)).

Some of the features described above are peculiar of the initial random condition. We have also studied the same Allen-Cahn problem, with the same parameters and the same boundary conditions, with a deterministic stratified initial state. Results for the $\epsilon$-profile are depicted at times in the figure 7; we have not reported the $m$-profiles and the $m-\epsilon$-plane view since they do not add any more information.

At times of order $10^{-1}$ a deviation from the initial stratified state is detected. The profile deforms into a new one with a droplet-like shape, at times of order $10^1$, finally an interface profile nucleates at times of order $10^2$, as in the previous case corresponding to spatially random initial data. Similarly to this case the stationary interface is approached with the same velocity. Even in this case we have checked the way in which the solution tends to the stationary state via the energy functional (41). Results are reported in figure 8.

### 4.4. Evolution under pure Darcy dissipation: Cahn-Hilliard-like equations

The equations of motion describing the evolution of the system subjected to a pure Darcy drag force are the equations (38) which, for the special model (45), read

$$ \begin{cases} -(2/3)\alpha bm^3 + \alpha b^2 m^2 \epsilon + p + \epsilon - ab(m - b\epsilon) - k_1 \epsilon'' - k_2 m'' = 0 \\ \varrho_{0,f}^2 (\alpha m^3 - 2\alpha bm^2 \epsilon + \alpha b^2 m \epsilon^2 + a(m - b\epsilon) - k_2 \epsilon'' - k_3 m'')'' = D\dot{m} \end{cases} \quad (54) $$

with the boundary conditions

$$ \begin{cases} ((k_2 m' + k_1 \epsilon') \delta \epsilon + (k_3 m' + k_2 \epsilon') \delta m)_{\ell_1, \ell_2} = 0 \\ (\alpha m^3 - 2\alpha bm^2 \epsilon + \alpha b^2 m \epsilon^2 + a(m - b\epsilon) - k_2 \epsilon'' - k_3 m'')_{\ell_1, \ell_2} = 0 \end{cases} \quad (55) $$

We have studied the above PDE problem with the same numerical approach used for the stationary problem and for the Allen-Cahn-like system. We have
repeated the same sequence of numerical computations, that is we have used
the same parameters and the same initial and boundary conditions as for the
above Allen–Cahn–like case.

In the figures 9–11 we have depicted the solutions of the problem (54) with Dirichlet boundary conditions $m(0) = m_s$, $\varepsilon(0) = \varepsilon_s$, $m(1) = m_f$, and $\varepsilon(1) = \varepsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. All the graphs refer to the same initial condition chosen randomly on the interval $[0, 1]$ as that considered for the solution of the Allen–Cahn–like equations. Details on the depicted times can be found in the caption of the related figures.

Figures 9 and 10 show that, starting from a random spatial distribution of strain and fluid mass, micro-structured paths appear at times of order $10^{-2}$. They progressively evolve towards a droplet-like profile which is formed at times of order $10^{-1}$. The droplet undergoes a metastable evolution at times of order $10^1$ up to the formation of an interface. This coarsening process has been well studied for the standard Cahn-Hilliard equation, see for instance [Ward 1998, Alikakos et al. 1991, Bates et al. 1994, Bates et al. 1995]. The last part of the evolution is a part of the motion of the interface approaching the stationary state at times of order $10^2$.

Figure 11 shows that even the solution of the Cahn-Hilliard-like system (54), as that of the Allen-Cahn-like system (52), is attracted immediately by the manifold in the plane $m-\varepsilon$ defined by the stationary solution and, afterwards, the projection of the profiles onto the $m-\varepsilon$ plane keeps close to such a manifold. Note that the profile temporarily leaves the stationary manifold at time 3.7 in correspondence of the shrinking of the droplet present in the $\varepsilon$-profile.

As for the Allen–Cahn–like case, also in the Cahn–Hilliard–like case we have studied the way in which the stationary state is approached by using the energy–like functional (41) (see figure 12).

Finally we have also studied the evolution of the system starting from a deterministic stratified initial datum; also in this case we used the same initial condition as the one used for the (52) problem, see figure 7. Results have been plotted in figures 13 and 14.

### 4.5. *Discussion of the numerical results*

Comparing the profiles of the strain and fluid mass density due to a pure Stokes (Cahn-Hilliard evolution) and a pure Darcy (Allen-Cahn evolution) dissipation, for spatially random or stratified initial conditions, some remarks naturally arise concerning short-time and long-time behavior.

Monitoring the evolution process described by the Cahn-Hilliard partial differential equations and considering spatially random initial data, short-time micro-scale oscillations can be detected at times of order $10^{-2}$; for these characteristic times the profile is smeared out and progressively deforms firstly into a metastable drop-like shape, which is approached at times of order $10^{-1}$, and lastly into an interface-like graph which is reached at times of order $10^1$. Analogously the short-time evolution process, described by the Allen-Cahn differential equations and stemming from the same randomly distributed initial condition, still exhibits a kind of intermediate feature and finally provides an interface-like graph which propagates towards the stationary state. However this last is approached at times several order of magnitude greater than those characteristic of the Cahn-Hilliard evolution process (in particular $\gg 10^4$).
repeated the same sequence of numerical computations, that is we have used
the same parameters and the same initial and boundary conditions as for the
above Allen–Cahn–like case.

In the figures 9–11 we have depicted the solutions of the problem (54) with Dirichlet boundary conditions $m(0) = m_s$, $\varepsilon(0) = \varepsilon_s$, $m(1) = m_f$, and $\varepsilon(1) = \varepsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. All the graphs refer to the same initial condition chosen randomly on the interval $[0, 1]$ as that considered for the solution of the Allen–Cahn–like equations. Details on the depicted times can be found in the caption of the related figures.

Figures 9 and 10 show that, starting from a random spatial distribution of strain and fluid mass, micro-structured paths appear at times of order $10^{-2}$. They progressively evolve towards a droplet-like profile which is formed at times of order $10^{-1}$. The droplet undergoes a metastable evolution at times of order $10^1$ up to the formation of an interface. This coarsening process has been well studied for the standard Cahn-Hilliard equation, see for instance [Ward 1998, Alikakos et al. 1991, Bates et al. 1994, Bates et al. 1995]. The last part of the evolution is a the motion of the interface approaching the stationary state at times of order $10^2$.

Figure 11 shows that even the solution of the Cahn-Hilliard-like system (54), as that of the Allen-Cahn-like system (52), is attracted immediately by the manifold in the plane $m-\varepsilon$ defined by the stationary solution and, afterwards, the projection of the profiles onto the $m-\varepsilon$ plane keeps close to such a manifold. Note that the profile temporarily leaves the stationary manifold at time 3.7 in correspondence of the shrinking of the droplet present in the $\varepsilon$-profile.

As for the Allen–Cahn–like case, also in the Cahn–Hilliard–like case we have studied the way in which the stationary state is approached by using the energy–like functional (41) (see figure 12).

Finally we have also studied the evolution of the system starting from a deterministic stratified initial datum; also in this case we used the same initial condition as the one used for the (52) problem, see figure 7. Results have been plotted in figures 13 and 14.

### 4.5. *Discussion of the numerical results*

Comparing the profiles of the strain and fluid mass density due to a pure Stokes (Cahn-Hilliard evolution) and a pure Darcy (Allen-Cahn evolution) dissipation, for spatially random or stratified initial conditions, some remarks naturally arise concerning short-time and long-time behavior.

Monitoring the evolution process described by the Cahn-Hilliard partial differential equations and considering spatially random initial data, short-time micro-scale oscillations can be detected at times of order $10^{-2}$; for these characteristic times the profile is smeared out and progressively deforms firstly into a metastable drop-like shape, which is approached at times of order $10^{-1}$, and lastly into an interface-like graph which is reached at times of order $10^1$. Analogously the short-time evolution process, described by the Allen-Cahn differential equations and stemming from the same randomly distributed initial condition, still exhibits a kind of intermediate feature and finally provides an interface-like graph which propagates towards the stationary state. However this last is approached at times several order of magnitude greater than those characteristic of the Cahn-Hilliard evolution process (in particular $\gg 10^4$).
The same difference is appreciated when considering stratified initial data; again the short-time behavior provides for both evolutions a deformation of the initial profile into a drop-like new one at different time scales. Looking at the profiles related to different characteristic times one can assess that the evolution process governed by the Allen-Cahn-like set of partial differential equations is affected by initial conditions during a characteristic time interval much larger than the corresponding one typical of the Cahn-Hilliard-like equations, see figures 3c and 9c, which correspond to the formation of the drop-like shape of the profile.

It is worth to notice that both the Allen-Cahn and the Cahn-Hilliard evolution processes separately exhibit a consistent behavior when the initial data is changed. In particular, starting from spatially random or stratified initial conditions, both equations provide, after a transitory regime, a quasi-static evolution of an interface-like profile which nucleates and propagates at the same time and with the same velocity. On the other hand the numerical simulations suggest that the quasi-static evolution process starts at different characteristic times and develops with different velocities when the Darcy and Stokes dissipation mechanisms are considered.

As already mentioned these two mechanisms correspond respectively to the presence (Darcy & Cahn-Hilliard regime) or absence (Stokes & Allen-Cahn regime) of scale separation and, thus, to the cases in which the fluid dissipative flow is parametrised either by the permeability of the solid or by the viscosity of the fluid. In these two situations the micro-structure of the corresponding porous media is quite different: materials with scale separation are constituted by solid grains which do not hinder the fluid flow strongly; the grains can be regarded as undeformable and impermeable, suffering an equivalent drag force due to the average fluid velocity. On the other hand materials with no scale separation are micro-porous and the fluid flow is much more hindered by the micro-structure of the skeleton. This could provide a possible explanation for the very high discrepancy between the characteristic times of evolution for the Cahn-Hilliard and the Allen-Cahn-like system of equation: when dealing with a Cahn-Hilliard dissipative process the time necessary to approach the stationary state is lower than that necessary to reach the same state through an Allen-Cahn evolution, see figures 6, 8, 12, and 14.

Comparing the distance between the current and the stationary strain energy for the Cahn-Hilliard and Allen-Cahn evolution process, three or two different regimes can be appreciated. These indeed correspond to the smearing out of the random oscillation (if any), to the formation of the intermediate drop-like solution, and to the nucleation of the interface-like profile.

Deeper investigations will be developed in the future in order to better explain such phenomenon, also accounting for different boundary conditions. Particular attention should be devoted to the study of initial boundary value problems involving Neumann boundary conditions, which naturally allow to describe the behavior of granular media under consolidation.
The same difference is appreciated when considering stratified initial data; again the short-time behavior provides for both evolutions a deformation of the initial profile into a drop-like new one at different time scales. Looking at the profiles related to different characteristic times one can assess that the evolution process governed by the Allen-Cahn-like set of partial differential equations is affected by initial conditions during a characteristic time interval much larger than the corresponding one typical of the Cahn-Hilliard-like equations, see figures 3c and 9c, which correspond to the formation of the drop-like shape of the profile.

It is worth to notice that both the Allen-Cahn and the Cahn-Hilliard evolution processes separately exhibit a consistent behavior when the initial data is changed. In particular, starting from spatially random or stratified initial conditions, both equations provide, after a transitory regime, a quasi-static evolution of an interface-like profile which nucleates and propagates at the same time and with the same velocity. On the other hand the numerical simulations suggest that the quasi-static evolution process starts at different characteristic times and develops with different velocities when the Darcy and Stokes dissipation mechanisms are considered.

As already mentioned these two mechanisms correspond respectively to the presence (Darcy & Cahn-Hilliard regime) or absence (Stokes & Allen-Cahn regime) of scale separation and, thus, to the cases in which the fluid dissipative flow is parametrised either by the permeability of the solid or by the viscosity of the fluid. In these two situations the micro-structure of the corresponding porous media is quite different: materials with scale separation are constituted by solid grains which do not hinder the fluid flow strongly; the grains can be regarded as undeformable and impermeable, suffering an equivalent drag force due to the average fluid velocity. On the other hand materials with no scale separation are micro-porous and the fluid flow is much more hindered by the micro-structure of the skeleton. This could provide a possible explanation for the very high discrepancy between the characteristic times of evolution for the Cahn-Hilliard and the Allen-Cahn-like system of equation: when dealing with a Cahn-Hilliard dissipative process the time necessary to approach the stationary state is lower than that necessary to reach the same state through an Allen-Cahn evolution, see figures 6, 8, 12, and 14.

Comparing the distance between the current and the stationary strain energy for the Cahn-Hilliard and Allen-Cahn evolution process, three or two different regimes can be appreciated. These indeed correspond to the smearing out of the random oscillation (if any), to the formation of the intermediate drop-like solution, and to the nucleation of the interface-like profile.

Deeper investigations will be developed in the future in order to better explain such phenomenon, also accounting for different boundary conditions. Particular attention should be devoted to the study of initial boundary value problems involving Neumann boundary conditions, which naturally allow to describe the behavior of granular media under consolidation.
<figure><img src="image_1.png" /><figcaption>Figure 1: Bottom: solutions ($\epsilon(X_S)$ on the left and $m(X_S)$ on the right) of the stationary problem (47) with the Dirichelet boundary conditions $m(0) = m_S$, $\epsilon(0) = \epsilon_S$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = 10^{-3}$, $k_3 = 10^{-3}$, $k_2 = \pm 10^{-3}$ (solid lines), and $k_2 = -0.4 \times 10^{-3}$, $0.2 \times 10^{-3}$, $0.8 \times 10^{-3}$ (dotted lines). Top: the same for $k_1 = 10^{-2}$, $k_3 = 10^{-2}$, $k_2 = \pm 10^{-2}$ (solid lines), and $k_2 = -0.4 \times 10^{-2}$, $0.2 \times 10^{-2}$, $0.8 \times 10^{-2}$ (dotted lines).</figcaption></figure>

## Appendix A. Stokes virtual dissipation in terms of the basic fields

Derivation of the equation (15) for the virtual work of the dissipation due to the Stokes effect. We first note that from the equality (10) we have that

$$[v_f(x, t) - v_S(x, t)]' = -\left(\left(\frac{\dot{\phi}}{\phi'}\chi_S'\right)'\right)_{\chi_S^{-1}(x,t)} (\chi^{-1})'_S(x, t) = \left(\frac{1}{\chi_S'}\left(\frac{\dot{\phi}}{\phi'}\chi_S'\right)'\right)_{\chi_S^{-1}(x,t)} \quad (A.1)$$

where in the first step we used the chain rule for the derivative of a composed function and in the second step the expression of the derivative of the inverse function. In the second term of the equation above the notation means that the derivative of $\phi\chi_S'/\phi'$ is computed with respect to $X_S$ and the result is then evaluated in $X_S = \chi_S^{-1}(x, t)$.

Now, we note that by evaluating the expression (9) for the variation of the field $\chi_f$ in $X_S = \chi_S^{-1}(x, t)$, for any $x \in B_t$ and $t \in \mathbb{R}$, we get

$$\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t) = -\frac{\chi_S'(\chi_S^{-1}(x, t), t)}{\phi'(\chi_S^{-1}(x, t), t)} \delta\phi(\chi_S^{-1}(x, t), t) \quad (A.2)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Bottom: solutions ($\epsilon(X_S)$ on the left and $m(X_S)$ on the right) of the stationary problem (47) with the Dirichelet boundary conditions $m(0) = m_S$, $\epsilon(0) = \epsilon_S$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = 10^{-3}$, $k_3 = 10^{-3}$, $k_2 = \pm 10^{-3}$ (solid lines), and $k_2 = -0.4 \times 10^{-3}$, $0.2 \times 10^{-3}$, $0.8 \times 10^{-3}$ (dotted lines). Top: the same for $k_1 = 10^{-2}$, $k_3 = 10^{-2}$, $k_2 = \pm 10^{-2}$ (solid lines), and $k_2 = -0.4 \times 10^{-2}$, $0.2 \times 10^{-2}$, $0.8 \times 10^{-2}$ (dotted lines).</figcaption></figure>

## Appendix A. Stokes virtual dissipation in terms of the basic fields

Derivation of the equation (15) for the virtual work of the dissipation due to the Stokes effect. We first note that from the equality (10) we have that

$$[v_f(x, t) - v_S(x, t)]' = -\left(\left(\frac{\dot{\phi}}{\phi'}\chi_S'\right)'\right)_{\chi_S^{-1}(x,t)} (\chi^{-1})'_S(x, t) = \left(\frac{1}{\chi_S'}\left(\frac{\dot{\phi}}{\phi'}\chi_S'\right)'\right)_{\chi_S^{-1}(x,t)} \quad (A.1)$$

where in the first step we used the chain rule for the derivative of a composed function and in the second step the expression of the derivative of the inverse function. In the second term of the equation above the notation means that the derivative of $\phi\chi_S'/\phi'$ is computed with respect to $X_S$ and the result is then evaluated in $X_S = \chi_S^{-1}(x, t)$.

Now, we note that by evaluating the expression (9) for the variation of the field $\chi_f$ in $X_S = \chi_S^{-1}(x, t)$, for any $x \in B_t$ and $t \in \mathbb{R}$, we get

$$\delta\chi_f(\chi_f^{-1}(x, t), t) - \delta\chi_S(\chi_S^{-1}(x, t), t) = -\frac{\chi_S'(\chi_S^{-1}(x, t), t)}{\phi'(\chi_S^{-1}(x, t), t)} \delta\phi(\chi_S^{-1}(x, t), t) \quad (A.2)$$
<figure><img src="image_2.png" /><figcaption>Figure 2: Plot on the $m-\epsilon$ plane of the solutions of the stationary problem (47) already depicted in the bottom row of the figure 1.</figcaption></figure>

As before, from (A.2) we get

$$[\delta\chi_f(\chi_f^{-1}(x,t),t) - \delta\chi_s(\chi_s^{-1}(x,t),t)]' = \left(\frac{1}{\chi_s'}\left(\frac{\chi_s'}{\phi'}\delta\phi\right)\right)'_{\chi_s^{-1}(x,t)} \quad (A.3)$$

Finally, by inserting (A.1) and (A.3) in (14) and by performing the change of variables $x = \chi_s(X_s, t)$ we get (15).

## Appendix B. Finite difference approximations

In this appendix we collect the finite difference substitution rules that we have adopted in our numerical computations. Let $n$ a positive integer number. Let $\sigma = 1/n$ and $\tau \in \mathbb{R}$ be respectively the space and the time increments. The space interval $[0, 1]$ is subdivided into $n$ small intervals of length $\sigma$. Given a field $h(X_s, t)$, for any $i = 1, \dots, n-1$ and $j \in \mathbb{N}$, we set

$$h'(i\sigma, j\tau) \approx \frac{1}{2\sigma} [h((i+1)\sigma, j\tau) - h((i-1)\sigma, j\tau)]$$

For the second space derivative we set

$$h''(i\sigma, j\tau) \approx \frac{1}{\sigma^2} [h((i+1)\sigma, j\tau) - 2h(i\sigma, j\tau) + h((i-1)\sigma, j\tau)]$$

for $i = 1, \dots, n-1$ and $j \in \mathbb{N}$. For the third space derivative we set

$$h'''(i\sigma, j\tau) \approx \frac{1}{2\sigma^3} [h((i+2)\sigma, j\tau) - 2h((i+1)\sigma, j\tau) \\ + 2h((i-1)\sigma, j\tau) - h((i-2)\sigma, j\tau)]$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Plot on the $m-\epsilon$ plane of the solutions of the stationary problem (47) already depicted in the bottom row of the figure 1.</figcaption></figure>

As before, from (A.2) we get

$$[\delta\chi_f(\chi_f^{-1}(x,t),t) - \delta\chi_s(\chi_s^{-1}(x,t),t)]' = \left(\frac{1}{\chi_s'}\left(\frac{\chi_s'}{\phi'}\delta\phi\right)\right)'_{\chi_s^{-1}(x,t)} \quad (A.3)$$

Finally, by inserting (A.1) and (A.3) in (14) and by performing the change of variables $x = \chi_s(X_s, t)$ we get (15).

## Appendix B. Finite difference approximations

In this appendix we collect the finite difference substitution rules that we have adopted in our numerical computations. Let $n$ a positive integer number. Let $\sigma = 1/n$ and $\tau \in \mathbb{R}$ be respectively the space and the time increments. The space interval $[0, 1]$ is subdivided into $n$ small intervals of length $\sigma$. Given a field $h(X_s, t)$, for any $i = 1, \dots, n-1$ and $j \in \mathbb{N}$, we set

$$h'(i\sigma, j\tau) \approx \frac{1}{2\sigma}[h((i+1)\sigma, j\tau) - h((i-1)\sigma, j\tau)]$$

For the second space derivative we set

$$h''(i\sigma, j\tau) \approx \frac{1}{\sigma^2}[h((i+1)\sigma, j\tau) - 2h(i\sigma, j\tau) + h((i-1)\sigma, j\tau)]$$

for $i = 1, \dots, n-1$ and $j \in \mathbb{N}$. For the third space derivative we set

$$h'''(i\sigma, j\tau) \approx \frac{1}{2\sigma^3}[h((i+2)\sigma, j\tau) - 2h((i+1)\sigma, j\tau) \\ + 2h((i-1)\sigma, j\tau) - h((i-2)\sigma, j\tau)]$$
for $i = 2, \dots, n-2$ and $j \in \mathbb{N}$. For the fourth space derivative we set

$$h'''(i\sigma, j\tau) \approx \frac{1}{\sigma^4} [h((i+2)\sigma, j\tau) - 4h((i+1)\sigma, j\tau) + 6h(i\sigma, j\tau) - 4h((i-1)\sigma, j\tau) + h((i-2)\sigma, j\tau)]$$

for $i = 2, \dots, n-2$ and $j \in \mathbb{N}$. For the time derivative we finally set

$$\dot{h}(i\sigma, j\tau) \approx \frac{1}{\tau} [h(i\sigma, j\tau) - h(i\sigma, (j-1)\tau)]$$

for $i = 0, \dots, n$ and $j = 1, \dots$.

## References

[Alikakos et al. 1991] N. Alikakos, P.W. Bates, G. Fusco, "Slow motion for the Cahn-Hilliard equation in one space dimension." *Journal of Differential Equations* **90**, 81-135 (1991).

[Alikakos et al. 2008] N. Alikakos, G. Fusco, "On the Connection Problem for Potentials with Several Global Minima." *Indiana University Mathematics Journal* **57**, 1871-1906 (2008).

[Allaire 1997] G. Allaire, "One-phase Newtonian flow." In "Homogenization and Porous Media," U. Hornung (Editor), Springer-Verlag, New York, 1997.

[Allen et al. 1979] S.M. Allen and J.W. Cahn, "A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening." *Acta Metallurgica* **27**, 1085-1095 (1979).

[Andrade et al. 2011] J.E. Andrade, C.F. Avila, S.A. Hall, N. Lenoir, G. Viggiani "Multiscale modeling and characterization of granular matter: from grain kinematics to continuum mechanics," *Journal of the Mechanics and Physics of Solids* **59**, 237-250 (2011).

[Auriault et al. 2005] J.L. Auriault, C. Geindreau, C. Boutin "Filtration Law in Porous Media with Poor Separation of Scales." *Transport in Porous Media* **60**, 89-108 (2005).

[Baek et al. 2004] S. Baek, A.R. Srinivasa, "Diffusion of a fluid through an elastic solid undergoing large deformation," *International Journal of Non-linear Mechanics* **39**, 201-218 (2004).

[Bates et al. 1994] P.W. Bates, J. Xun, "Metastable patterns for the Cahn-Hilliard equation: part I," *Journal of Differential Equations* **111**, 421-457 (1994).

[Bates et al. 1995] P.W. Bates, J. Xun, "Metastable patterns for the Cahn-Hilliard equation: part II," *Journal of Differential Equations* **117**, 165-216 (1995).

[Bear 1972] J. Bear, "Dynamics of fluids in porous media," Elsevier, 1972.
for $i = 2, \dots, n-2$ and $j \in \mathbb{N}$. For the fourth space derivative we set

$$h'''(i\sigma, j\tau) \approx \frac{1}{\sigma^4} [h((i+2)\sigma, j\tau) - 4h((i+1)\sigma, j\tau) + 6h(i\sigma, j\tau) - 4h((i-1)\sigma, j\tau) + h((i-2)\sigma, j\tau)]$$

for $i = 2, \dots, n-2$ and $j \in \mathbb{N}$. For the time derivative we finally set

$$\dot{h}(i\sigma, j\tau) \approx \frac{1}{\tau} [h(i\sigma, j\tau) - h(i\sigma, (j-1)\tau)]$$

for $i = 0, \dots, n$ and $j = 1, \dots$.

## References

[Alikakos et al. 1991] N. Alikakos, P.W. Bates, G. Fusco, "Slow motion for the Cahn-Hilliard equation in one space dimension." *Journal of Differential Equations* **90**, 81-135 (1991).

[Alikakos et al. 2008] N. Alikakos, G. Fusco, "On the Connection Problem for Potentials with Several Global Minima." *Indiana University Mathematics Journal* **57**, 1871-1906 (2008).

[Allaire 1997] G. Allaire, "One-phase Newtonian flow." In "Homogenization and Porous Media," U. Hornung (Editor), Springer-Verlag, New York, 1997.

[Allen et al. 1979] S.M. Allen and J.W. Cahn, "A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening." *Acta Metallurgica* **27**, 1085-1095 (1979).

[Andrade et al. 2011] J.E. Andrade, C.F. Avila, S.A. Hall, N. Lenoir, G. Viggiani "Multiscale modeling and characterization of granular matter: from grain kinematics to continuum mechanics," *Journal of the Mechanics and Physics of Solids* **59**, 237-250 (2011).

[Auriault et al. 2005] J.L. Auriault, C. Geindreau, C. Boutin "Filtration Law in Porous Media with Poor Separation of Scales." *Transport in Porous Media* **60**, 89-108 (2005).

[Baek et al. 2004] S. Baek, A.R. Srinivasa, "Diffusion of a fluid through an elastic solid undergoing large deformation," *International Journal of Non-linear Mechanics* **39**, 201-218 (2004).

[Bates et al. 1994] P.W. Bates, J. Xun, "Metastable patterns for the Cahn-Hilliard equation: part I," *Journal of Differential Equations* **111**, 421-457 (1994).

[Bates et al. 1995] P.W. Bates, J. Xun, "Metastable patterns for the Cahn-Hilliard equation: part II," *Journal of Differential Equations* **117**, 165-216 (1995).

[Bear 1972] J. Bear, "Dynamics of fluids in porous media," Elsevier, 1972.
[Besuelle 2001] P. Besuelle, "Compacting and dilating shear bands in porous rock: Theoretical and experimental conditions," *Journal of Geophysical Research*, **106**, 13,435–13,442 (2001).

[Biot 1941] M.A. Biot, "General Theory of Three-Dimensional Consolidation." *Journal of Applied Physics* **12**, 155 (1941).

[Blanchard et al. 1992] P. Blanchard, E. Brüning, "Variational methods in mathematical physics." Springer, New York, US, 1992.

[Bray 1994] A.J. Bray, "Theory of phase-ordering kinetics." *Advances in Physics* **43**, 357–459 (1994).

[Cahn 1961] J.W. Cahn, "On spinodal decomposition." *Acta Metallurgica* **9**, 795–801 (1961).

[Carr et al. 1989] J. Carr, B. Pego, "Metastable patterns is solutions of $u_t = \epsilon^2 u_{xx} - f(x)$," *Communications on Pure and Applied Mathematics* **42**, 523–576 (1989).

[Chambon et al. 2001] R. Chambon, D. Caillerie, T. Matschima, "Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies," *International Journal of Solids and Structures* **38**, 8503–8527 (2001).

[Christoffersen et al. 1981] J. Christoffersen, M.M. Mehrabadi, S. Nemat-Nasser, "A micromechanical description of granular material behavior," *Journal of Applied Mechanics* **48**, 339–344 (1981).

[Cirillo et al. 2009] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Solid-fluid segregation in saturated porous media." *Poromechanics IV: Proceedings of the Fourth BIOT Conference on Poromechanics, DEStech Publications, Inc.*, June 2009.

[Cirillo et al. 2010] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Phase coexistence in consolidating porous media." *Physical Review E* **81**, 061121 (2010).

[Cirillo et al. 2011] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Phase transition in saturated porous-media: pore-fluid segregation in consolidation." *Physica D* **240**, 1345–1351 (2011).

[Cirillo et al. 2012] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Kink Localization under Asymmetric Double-Well Potential." preprint arXiv:1207.5727v1 (2012).

[Collin et al. 2006] F. Collin, R. Chambon, R. Charlier, "A finite element method for poro mechanical modelling of geotechnical problems using local second gradient models," *International Journal for Numerical Methods in Engineering* **65**, 1749–1772 (2006).

[Convertino et al. 2002] A. Convertino, A. Valentini, A. Bassi, N. Cioffi, L. Torsi, and E.N.M. Cirillo, "Effect of metal clusters on the swelling of gold-fluorocarbon-polymer composite films," *Applied Physics Letters* **80**, 1565–1567 (2002).
[Besuelle 2001] P. Besuelle, "Compacting and dilating shear bands in porous rock: Theoretical and experimental conditions," *Journal of Geophysical Research*, **106**, 13,435–13,442 (2001).

[Biot 1941] M.A. Biot, "General Theory of Three-Dimensional Consolidation." *Journal of Applied Physics* **12**, 155 (1941).

[Blanchard et al. 1992] P. Blanchard, E. Brüning, "Variational methods in mathematical physics." Springer, New York, US, 1992.

[Bray 1994] A.J. Bray, "Theory of phase-ordering kinetics." *Advances in Physics* **43**, 357–459 (1994).

[Cahn 1961] J.W. Cahn, "On spinodal decomposition." *Acta Metallurgica* **9**, 795–801 (1961).

[Carr et al. 1989] J. Carr, B. Pego, "Metastable patterns is solutions of $u_t = \epsilon^2 u_{xx} - f(x)$," *Communications on Pure and Applied Mathematics* **42**, 523–576 (1989).

[Chambon et al. 2001] R. Chambon, D. Caillerie, T. Matschima, "Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies," *International Journal of Solids and Structures* **38**, 8503–8527 (2001).

[Christoffersen et al. 1981] J. Christoffersen, M.M. Mehrabadi, S. Nemat-Nasser, "A micromechanical description of granular material behavior," *Journal of Applied Mechanics* **48**, 339–344 (1981).

[Cirillo et al. 2009] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Solid-fluid segregation in saturated porous media." *Poromechanics IV: Proceedings of the Fourth BIOT Conference on Poromechanics, DEStech Publications, Inc.*, June 2009.

[Cirillo et al. 2010] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Phase coexistence in consolidating porous media." *Physical Review E* **81**, 061121 (2010).

[Cirillo et al. 2011] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Phase transition in saturated porous-media: pore-fluid segregation in consolidation." *Physica D* **240**, 1345–1351 (2011).

[Cirillo et al. 2012] E.N.M. Cirillo, N. Ianiro, G. Sciarra, "Kink Localization under Asymmetric Double-Well Potential." preprint arXiv:1207.5727v1 (2012).

[Collin et al. 2006] F. Collin, R. Chambon, R. Charlier, "A finite element method for poro mechanical modelling of geotechnical problems using local second gradient models," *International Journal for Numerical Methods in Engineering* **65**, 1749–1772 (2006).

[Convertino et al. 2002] A. Convertino, A. Valentini, A. Bassi, N. Cioffi, L. Torsi, and E.N.M. Cirillo, "Effect of metal clusters on the swelling of gold-fluorocarbon-polymer composite films," *Applied Physics Letters* **80**, 1565–1567 (2002).
[Convertino et al. 2003] A. Convertino, A. Capobianchi, A. Valentini, and E.N.M. Cirillo, "A New Approach to Organic Solvent Detection: High reflectivity Bragg reflectors based on a gold nanoparticle/Teflon-like composite material," *Advanced Materials* **15**, 1103 (2003).

[Coussy 2004] O. Coussy, "Poromechanics." John Wiley & Sons, Ltd, 2004.

[Coussy 2010] O. Coussy, "Mechanics and Physics of porous solids." John Wiley and Sons, 2010.

[Cryer 1963] C.W. Cryer, "A comparison of the three-dimensional consolidation theories of Biot and Terzaghi," *The Quarterly Journal of Mechanics and Applied Mathematics* **16**, 401–412 (1963).

[dell'Isola et al. 2009] F. dell'Isola, G. Sciarra, S. Vidoli, "Generalized Hooke's law for isotropic second gradient materials," *Proceeding of the Royal Society A* **465**, 2177–2196 (2009).

[Ene et al. 1975] H.I. Ene, E. Sanchez-Palencia, "Equations et phénomènes de surface pour l'ècoulement dans un modèle de milieux poreux," *Journal de Mécanique* **14**, 73–108 (1975).

[Eyre 1993] D.J. Eyre, "Systems of Cahn–Hilliard equations." *SIAM Journal on Applied Mathematics* **53**, 1686–1712 (1993).

[Fife 2002] P.C. Fife, "Pattern formation in gradient systems," in "Handbook of Dynamical Systems" vol. 2, pages 677–722, 2002, ed. B. Fielder, Elsevier Science, Amsterdam, The Netherlands.

[Fleck et al. 1997] N.A. Fleck, J.W. Hutchinson, "Strain gradient plasticity," *Advances in Applied Mechanics* **33**, 295–361 (1997).

[Fleck et al. 2001] N.A. Fleck, J.W. Hutchinson, "A reformulation of strain gradient plasticity," *Journal of the Mechanics and Physics of Solids* **49**, 2245–2271 (2001).

[Fusco et al. 1989] G. Fusco, J.K. Hale, "Slow motion manifolds, dormant instability and singular perturbations," *Journal of Dynamics and Differential Equations* **1**, 75–94 (1989).

[Gantmacher 1975] F. Gantmacher, "Lectures in analytical mechanics." Mir Publishers, Moscow, 1975.

[Gelb et al. 1999] L.D. Gelb, K.E. Gubbins, R. Radhakrishnan, M. Sliwinska-Bartkoviak, "Phase separation in confined systems." *Reports on Progress in Physics* **62**, 1573–1659 (1999).

[Germain 1973] P. Germain, "La méthode des puissances virtuelles en mécanique des milieux continus. Première partie: Théorie du second gradient." *Journal de Mécanique* **12**, 235–274 (1973).

[Gurtin 1996] M.E. Gurtin, "Generalized Ginzburg–Landau and Cahn–Hilliard equations based on a microforce balance." *Physica D* **92**, 178–192 (1996).
[Convertino et al. 2003] A. Convertino, A. Capobianchi, A. Valentini, and E.N.M. Cirillo, "A New Approach to Organic Solvent Detection: High reflectivity Bragg reflectors based on a gold nanoparticle/Teflon-like composite material," *Advanced Materials* **15**, 1103 (2003).

[Coussy 2004] O. Coussy, "Poromechanics." John Wiley & Sons, Ltd, 2004.

[Coussy 2010] O. Coussy, "Mechanics and Physics of porous solids." John Wiley and Sons, 2010.

[Cryer 1963] C.W. Cryer, "A comparison of the three-dimensional consolidation theories of Biot and Terzaghi," *The Quarterly Journal of Mechanics and Applied Mathematics* **16**, 401–412 (1963).

[dell'Isola et al. 2009] F. dell'Isola, G. Sciarra, S. Vidoli, "Generalized Hooke's law for isotropic second gradient materials," *Proceeding of the Royal Society A* **465**, 2177–2196 (2009).

[Ene et al. 1975] H.I. Ene, E. Sanchez-Palencia, "Equations et phénomènes de surface pour l'ècoulement dans un modèle de milieux poreux," *Journal de Mécanique* **14**, 73–108 (1975).

[Eyre 1993] D.J. Eyre, "Systems of Cahn–Hilliard equations." *SIAM Journal on Applied Mathematics* **53**, 1686–1712 (1993).

[Fife 2002] P.C. Fife, "Pattern formation in gradient systems," in "Handbook of Dynamical Systems" vol. 2, pages 677–722, 2002, ed. B. Fielder, Elsevier Science, Amsterdam, The Netherlands.

[Fleck et al. 1997] N.A. Fleck, J.W. Hutchinson, "Strain gradient plasticity," *Advances in Applied Mechanics* **33**, 295–361 (1997).

[Fleck et al. 2001] N.A. Fleck, J.W. Hutchinson, "A reformulation of strain gradient plasticity," *Journal of the Mechanics and Physics of Solids* **49**, 2245–2271 (2001).

[Fusco et al. 1989] G. Fusco, J.K. Hale, "Slow motion manifolds, dormant instability and singular perturbations," *Journal of Dynamics and Differential Equations* **1**, 75–94 (1989).

[Gantmacher 1975] F. Gantmacher, "Lectures in analytical mechanics." Mir Publishers, Moscow, 1975.

[Gelb et al. 1999] L.D. Gelb, K.E. Gubbins, R. Radhakrishnan, M. Sliwinska-Bartkoviak, "Phase separation in confined systems." *Reports on Progress in Physics* **62**, 1573–1659 (1999).

[Germain 1973] P. Germain, "La méthode des puissances virtuelles en mécanique des milieux continus. Première partie: Théorie du second gradient." *Journal de Mécanique* **12**, 235–274 (1973).

[Gurtin 1996] M.E. Gurtin, "Generalized Ginzburg–Landau and Cahn–Hilliard equations based on a microforce balance." *Physica D* **92**, 178–192 (1996).
[Kolymbas 1994] D. Kolymbas, "Compaction waves as phase transitions," *Acta Mechanica* **107**, 171-181 (1994).

[Kolymbas 1998] D. Kolymbas, "Behavior of liquefied sand," *Phil. Trans. R. Soc. Lond. A* **356**, 2609-2622 (1998)

[Laforgue et al. 1998] J.G.L. Laforgue, R.E. O'Malley, M.J. Ward, "Metastable travelling-wave solutions of singularly-perturbed reaction-diffusion equations." *European Journal of Applied Mathematics* **9**, 397-416 (1998).

[Langer 1992] J.S. Langer, *Solids far from Equilibrium*, edited by C. Godrèche, Cambridge University Press, page 297, 1992.

[Lopatnikov et al. 2004] S.L. Lopatnikov, A.H.-D. Cheng, "Macroscopic Lagrangian formulation of poroelasticity with porosity dynamics." *Journal of the Mechanics and Physics of Solids* **52**, 2801-2839 (2004).

[Madeo et al. 2008] A. Madeo, F. Dell'Isola, N. Ianiro, and G. Sciarra, "A variational deduction of second gradient poroelasticity. Part II: an application to the consolidation problem." *Journal of Mechanics of Materials and Structures* **3**, 607-625 (2008).

[Mandel 1953] J. Mandel, "Consolidation des sols (étude mathématique)," *Géothèchnique* **3**, 287-299 (1953).

[Nemat-Nasser et al. 1993] S. Nemat-Nasser, M. Hori, "Micromechanics: Overall properties of heterogeneous materials." North Holland, 1993.

[Nichols et al. 1994] 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).

[Nield 2007] D.A. Nield, "The modeling of viscous dissipation in a saturated porous medium," *Journal of Heat Transfer-Transactions of the ASME* **129**, 1459-1463 (2007).

[Sciarra et al. 2008] 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**, 507-526 (2008).

[Seppecher 1989] P. Seppecher, "Etude des conditions aux limites en théorie du second gradient: cas de la capillarité." *C.R. Acad. Sci. Paris* **309** Serie II, 497-502 (1989).

[von Terzaghi 1946] K. von Terzaghi, *Theoretical Soil Mechanics*. John Wiley and Sons, 1946.

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

[Vardoulakis 2004b] I. Vardoulakis, "Fluidisation in artesian flow conditions: Hydromechanically unstable granular media," *Géotechnique* **54**, 165-177 (2004).
[Kolymbas 1994] D. Kolymbas, "Compaction waves as phase transitions," *Acta Mechanica* **107**, 171-181 (1994).

[Kolymbas 1998] D. Kolymbas, "Behavior of liquefied sand," *Phil. Trans. R. Soc. Lond. A* **356**, 2609-2622 (1998)

[Laforgue et al. 1998] J.G.L. Laforgue, R.E. O'Malley, M.J. Ward, "Metastable travelling-wave solutions of singularly-perturbed reaction-diffusion equations." *European Journal of Applied Mathematics* **9**, 397-416 (1998).

[Langer 1992] J.S. Langer, *Solids far from Equilibrium*, edited by C. Godrèche, Cambridge University Press, page 297, 1992.

[Lopatnikov et al. 2004] S.L. Lopatnikov, A.H.-D. Cheng, "Macroscopic Lagrangian formulation of poroelasticity with porosity dynamics." *Journal of the Mechanics and Physics of Solids* **52**, 2801-2839 (2004).

[Madeo et al. 2008] A. Madeo, F. Dell'Isola, N. Ianiro, and G. Sciarra, "A variational deduction of second gradient poroelasticity. Part II: an application to the consolidation problem." *Journal of Mechanics of Materials and Structures* **3**, 607-625 (2008).

[Mandel 1953] J. Mandel, "Consolidation des sols (étude mathématique)," *Géothecnique* **3**, 287-299 (1953).

[Nemat-Nasser et al. 1993] S. Nemat-Nasser, M. Hori, "Micromechanics: Overall properties of heterogeneous materials." North Holland, 1993.

[Nichols et al. 1994] 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).

[Nield 2007] D.A. Nield, "The modeling of viscous dissipation in a saturated porous medium," *Journal of Heat Transfer-Transactions of the ASME* **129**, 1459-1463 (2007).

[Sciarra et al. 2008] 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**, 507-526 (2008).

[Seppecher 1989] P. Seppecher, "Etude des conditions aux limites en théorie du second gradient: cas de la capillarité." *C.R. Acad. Sci. Paris* **309** Serie II, 497-502 (1989).

[von Terzaghi 1946] K. von Terzaghi, *Theoretical Soil Mechanics*. John Wiley and Sons, 1946.

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

[Vardoulakis 2004b] I. Vardoulakis, "Fluidisation in artesian flow conditions: Hydromechanically unstable granular media," *Géotechnique* **54**, 165-177 (2004).
[Ward 1994] M.J. Ward, “Metastable patterns, layer collapse, and coarsening for a one-dimensional Ginzburg–Landau equation,” *Studies in Applied Mathematics* **91**, 51–93 (1994).

[Ward 1998] M.J. Ward, “Dynamic Metastability and Singular Perturbations,” in “Boundaries, Interfaces and Transitions,” M.C. Delfour ed., CRM proceedings and lecture notes, vol. 13, American Mathematical Society, Providence, R.I., US, 1998.

[Zhang et al. 1999] H.W. Zhang, L. Sanavia and B.A. Schrefler, “An internal length scale in dynamic strain localization of multiphase porous media,” *Mechanics of Cohesive-Frictional Materials* **4**, 443–460 (1999).
[Ward 1994] M.J. Ward, “Metastable patterns, layer collapse, and coarsening for a one-dimensional Ginzburg–Landau equation,” *Studies in Applied Mathematics* **91**, 51–93 (1994).

[Ward 1998] M.J. Ward, “Dynamic Metastability and Singular Perturbations,” in “Boundaries, Interfaces and Transitions,” M.C. Delfour ed., CRM proceedings and lecture notes, vol. 13, American Mathematical Society, Providence, R.I., US, 1998.

[Zhang et al. 1999] H.W. Zhang, L. Sanavia and B.A. Schrefler, “An internal length scale in dynamic strain localization of multiphase porous media,” *Mechanics of Cohesive-Frictional Materials* **4**, 443–460 (1999).
<figure><img src="image_3.png" /><figcaption>Figure 3: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Allen-Cahn-like system (52) with a random initial condition and Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 2, 5.4, 7, 15, 200, 50000$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Allen-Cahn-like system (52) with a random initial condition and Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 2, 5.4, 7, 15, 200, 50000$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4: The same as in figure 3 for the $m(X_s, t)$ profiles.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: The same as in figure 3 for the $m(X_s, t)$ profiles.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: The solution of the same problem as in figure 3 is depicted on the plane $m-\epsilon$ at times $t = 5.4, 7, 15$ from the left to the right.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: For the same problem as in figure 3 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 2, 5.4, 7, 15, 200$ at which the $\epsilon$ and $m$-profiles are depicted in figures 3–5.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: The solution of the same problem as in figure 3 is depicted on the plane $m-\epsilon$ at times $t = 5.4, 7, 15$ from the left to the right.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: For the same problem as in figure 3 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 2, 5.4, 7, 15, 200$ at which the $\epsilon$ and $m$-profiles are depicted in figures 3–5.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Allen-Cahn-like system (52) with a deterministic stratified initial state and Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 0.07, 2, 10, 320, 340, 50000$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Allen-Cahn-like system (52) with a deterministic stratified initial state and Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 0.07, 2, 10, 320, 340, 50000$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: For the same problem as in figure 7 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 0.07, 2, 10, 320, 340$ at which the $\epsilon$ profiles are depicted in figures 7.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: For the same problem as in figure 7 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 0.07, 2, 10, 320, 340$ at which the $\epsilon$ profiles are depicted in figures 7.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 9: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Cahn-Hilliard-like system (54) with the same random initial state as the one used in figure 3. We used Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 0.009, 0.04, 0.2, 3.7, 3.9, 89$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Cahn-Hilliard-like system (54) with the same random initial state as the one used in figure 3. We used Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 0.009, 0.04, 0.2, 3.7, 3.9, 89$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 10: The same as in figure 9 for the profiles $m(X_s, t)$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: The same as in figure 9 for the profiles $m(X_s, t)$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: The solution of the same problem as in figure 9 is depicted on the plane $m-\epsilon$ at times $t = 0.04, 0.2, 3.7$ from the left to the right.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: For the same problem as in figure 9 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 0.009, 0.04, 0.2, 3.7, 3.9, 89$ at which the $\epsilon$ and $m$-profiles are depicted in figures 9–11.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: The solution of the same problem as in figure 9 is depicted on the plane $m-\epsilon$ at times $t = 0.04, 0.2, 3.7$ from the left to the right.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: For the same problem as in figure 9 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 0.009, 0.04, 0.2, 3.7, 3.9, 89$ at which the $\epsilon$ and $m$-profiles are depicted in figures 9–11.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 13: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Cahn-Hilliard-like system (54) with a deterministic stratified initial state (the same used in figure 7) and Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 0.001, 0.01, 10.3, 11, 11.1, 89$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Profiles (solid thin lines) $\epsilon(X_s, t)$ obtained by solving the Cahn-Hilliard-like system (54) with a deterministic stratified initial state (the same used in figure 7) and Dirichlet boundary conditions $m(0) = m_s$, $\epsilon(0) = \epsilon_s$, $m(1) = m_f$, and $\epsilon(1) = \epsilon_f$ on the finite interval $[0, 1]$, at the coexistence pressure for $a = 0.5$, $b = 1$, $\alpha = 100$, $k_1 = k_2 = k_3 = 10^{-3}$. The solid thick line is the corresponding stationary profile. Profiles at times $t = 0.001, 0.01, 10.3, 11, 11.1, 89$ are depicted in lexicographic order.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: For the same problem as in figure 13 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 0.001, 0.01, 10.3, 11, 11.1, 89$ at which the $\epsilon$ profiles are depicted in figures 13.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: For the same problem as in figure 13 the difference of the energy (41) at time $t$ and that corresponding to the stationary profile is reported as function of time. The vertical thin lines denote the times $t = 0.001, 0.01, 10.3, 11, 11.1, 89$ at which the $\epsilon$ profiles are depicted in figures 13.</figcaption></figure>