Fluid-fluid phase separation in a soft porous medium

Fluid-flow problems are typically considered in an Eulerian reference frame (fixed in space), whereas solid-mechanics problems are typically considered in a Lagrangian reference frame (fixed to the solid material). For a poromechanical system, however, the fluids and the solid co-exist and must be considered in the same reference frame. To do so, it is convenient to begin by defining both frames and the relationship between them.

The deformation of the solid skeleton is defined by comparing its current (deformed) state to an undeformed (relaxed) reference state. The Lagrangian coordinate $\mathbf{X}$ refers to the reference state and the Eulerian coordinate $\mathbf{x}$ refers to the current state. We denote spatial gradients with respect to the Lagrangian and Eulerian coordinates by $\nabla_{0}$ and $\nabla$, respectively. The Lagrangian displacement of the solid is then $\mathbf{U}_{s}\left(\mathbf{X},t\right)=\mathbf{x}\left(\mathbf{X},t\right)-\mathbf{X}$ and the Eulerian displacement is $\mathbf{u}_{s}\left(\mathbf{x},t\right)=\mathbf{x}-\mathbf{X}\left(\mathbf{x},t\right)$. The deformation gradient tensor, $\mathbb{F}$, describes the deformation of the solid skeleton relative to its reference state,

$$\mathbb{F}=\nabla_{0}\mathbf{x}.$$

(1)

Note that we use the convention $\left(\nabla\mathbf{A}\right)_{ij}=\frac{\partial A_{i}}{\partial x_{j}}$ for a general vector field $\mathbf{A}=A_{i}\mathbf{\hat{e}}_{i}$ and $\left(\nabla\cdot\mathbb{B}\right)_{i}=\frac{\partial B_{ij}}{\partial x_{j}}$ for a general tensor field $\mathbb{B}=B_{ij}\mathbf{\hat{e}}_{i}\mathbf{\hat{e}}_{j}$, where $\mathbf{\hat{e}}_{i}$ and $\mathbf{\hat{e}}_{j}$ are Cartesian unit vectors and where we have adopted the Einstein summation convention. Using the definitions above, $\mathbb{F}$ is related to $\mathbf{U}_{s}$ and $\mathbf{u}_{s}$ via

$$\mathbb{F}=\mathbb{I}+\nabla_{0}\mathbf{U}_{s}=\left[\mathbb{I}-\nabla\mathbf{u}_{s}\right]^{-1},$$

(2)

where $\mathbb{I}$ is the identity tensor.

The volume fraction of each phase can be measured relative to either the relaxed or deformed configurations. The nominal volume fractions, $\varphi_{\alpha}$, measure the volume of phase $\alpha$ per unit reference bulk volume, whereas the true volume fractions, $\phi_{\alpha}$, measure the volume of phase $\alpha$ per unit current bulk volume, where $\alpha=s,~{}l,~{}g$ indicates the solid, liquid, and gas phases, respectively. These volume fractions are related by $\varphi_{\alpha}=\phi_{\alpha}J$, where the Jacobian determinant, $J=\textrm{det}(\mathbb{F})$, measures the ratio of current (deformed) bulk volume to reference (relaxed) bulk volume and is thus equal to one in the reference state. The no-void condition can thus be expressed in two ways,

$$\sum_{\alpha}\varphi_{\alpha}=J\quad\textrm{and}\quad\sum_{\alpha}\phi_{\alpha}=1.$$

(3)

A result of the no-void condition is that any two of the three volume fractions fully determine the local composition.

The incompressibility of the solid grains requires that the nominal volume fraction of solid must remain unchanged during deformation, such that $\varphi_{s}\equiv\varphi_{s}^{0}\equiv\phi_{s}^{0}$, where $\varphi_{s}^{0}\equiv\phi_{s}^{0}$ is the solid fraction in the reference state, which we take to be uniform for simplicity. As a result,

$$J=\frac{\phi_{s}^{0}}{\phi_{s}}.$$

(4)

It is also convenient to define the true porosity, $\phi=\phi_{g}+\phi_{l}=1-\phi_{s}$, and the relaxed porosity, $\phi_{0}=1-\phi_{s}^{0}$. Below, we develop evolution equations for $\phi_{g}$ and $\phi$ and then evaluate $\phi_{l}=\phi-\phi_{g}$ and $\phi_{s}=1-\phi$ as needed.

We formulate evolution equations for $\phi_{g}$ and $\phi$ by first considering conservation of mass for each phase. In the Eulerian frame, conservation of mass can be written

$$\frac{\partial}{\partial t}\left(\rho_{\alpha}\phi_{\alpha}\right)+\nabla\cdot\left(\rho_{\alpha}\phi_{\alpha}\mathbf{v}_{\alpha}\right)=0,$$

(5)

where $\rho_{\alpha}$ and $\mathbf{v}_{\alpha}$ are the density and local velocity of phase $\alpha$, respectively. We assume for simplicity that all three phases are incompressible, meaning that $\rho_{\alpha}$ is a constant. The Eulerian volume flux, $\mathbf{w}_{\alpha}$, of phase $\alpha$ relative to the solid skeleton is

$$\mathbf{w}_{\alpha}=\phi_{\alpha}\left(\mathbf{v}_{\alpha}-\mathbf{v}_{s}\right).$$

(6)

We reduce Equations (5) to two evolution equations, one each for $\phi_{g}$ and $\phi$, and an elliptic equation for the total mixture flux, $\mathbf{q}=\sum\limits_{\alpha}\phi_{\alpha}\mathbf{v}_{\alpha}=\mathbf{v}_{s}+\mathbf{w}_{g}+\mathbf{w}_{l}$, by eliminating the phase velocities in favour of the relative fluxes. The result of these manipulations is

	$\displaystyle\frac{\partial\phi_{g}}{\partial t}+\nabla\cdot\left(\phi_{g}\mathbf{q}\right)+\nabla\cdot\left[\left(1-\phi_{g}\right)\mathbf{w}_{g}-\phi_{g}\mathbf{w}_{l}\right]$	$\displaystyle=0,$		(7a)
	$\displaystyle\frac{\partial\phi}{\partial t}+\nabla\cdot\left(\phi\mathbf{q}\right)+\nabla\cdot\left[\left(1-\phi\right)\left(\mathbf{w}_{g}+\mathbf{w}_{l}\right)\right]$	$\displaystyle=0,$		(7b)
	$\displaystyle\nabla\cdot\mathbf{q}$	$\displaystyle=0.$		(7c)

We derive thermodynamically consistent constitutive expressions for the fluxes $\mathbf{w}_{\alpha}$ in §II.5. We express mass conservation in Lagrangian form by using the Reynolds transport theorem to rewrite Equation (5) as

$$\frac{\partial}{\partial t}\left(\rho_{\alpha}\varphi_{\alpha}\right)+\nabla_{0}\cdot\left(\rho_{\alpha}\mathbf{W}_{\alpha}\right)=0,$$

(8)

where the nominal volume flux, $\mathbf{W}_{\alpha}$, is related to the true flux by $\mathbf{W}_{\alpha}=J\mathbb{F}^{-1}\mathbf{w}_{\alpha}$.

Further restrictions can be imposed on our system by considering conservation of momentum. The true (Cauchy) total stress, $\mathbb{T}$, measures the total internal force per unit current area within the mixture. In the absence of body forces and neglecting inertia, mechanical equilibrium requires that $\mathbb{T}$ must be divergence free,

$$\nabla\cdot\mathbb{T}=0.$$

(9)

In the Lagrangian frame, momentum balance is most naturally expressed in terms of the first Piola-Kirchhoff stress, $\mathbb{S}$,

$$\nabla_{0}\cdot\mathbb{S}=0,$$

(10)

where $\mathbb{S}=J\mathbb{T}\mathbb{F}^{-\intercal}$ measures the total internal force per unit reference area within the mixture.

The first law of thermodynamics requires that the rate of change of free energy in a particular control volume, via mass flux into the volume or working due to internal compositional changes, cannot exceed the rate of work done by external forces. Following Gurtin [39], we formulate this restriction by considering the change in Helmholtz free energy, $\psi$, of a bulk material element in the Lagrangian frame. The rate of change of energy per unit mass of phase $\alpha$ due to a net flux of phase $\alpha$ is measured by the chemical potential, $\mu_{\alpha}$. The key ingredient in Gurtin’s approach is that changes in energy resulting from changes in composition, as measured by changes in $\varphi_{\alpha}$, can be represented by the action of a nominal vector ‘microstress’, $\bm{\xi}_{\alpha}$, acting on the boundary of the material element. This microstress is intended to capture the net work done by subscale physical mechanisms such as capillary forces that would otherwise not be resolved at the continuum scale. Mechanical damage to the solid skeleton can be captured in the same way by introducing a damage parameter, $d$, and an associated microstress, $\bm{\xi}_{d}$ (see §II.4 below). The working of external forces is represented by the action of the total stress, $\mathbb{S}$, on the boundary of the material element. Summing these different contributions leads to a dissipation inequality,

$$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{V}_{0}}\psi\mathrm{d}V\leq-\int_{\partial\mathcal{V}_{0}}\rho_{g}\mu_{g}\mathbf{W}_{g}\cdot\hat{\mathbf{n}}\mathrm{d}A-\int_{\partial\mathcal{V}_{0}}\rho_{l}\mu_{l}\mathbf{W}_{l}\cdot\hat{\mathbf{n}}\mathrm{d}A+\int_{\partial\mathcal{V}_{0}}\left(\bm{\xi}_{g}\cdot\hat{\mathbf{n}}\right)\dot{\varphi_{g}}\mathrm{d}A\\ +\int_{\partial\mathcal{V}_{0}}\left(\bm{\xi}_{l}\cdot\hat{\mathbf{n}}\right)\dot{\varphi_{l}}\mathrm{d}A+\int_{\partial\mathcal{V}_{0}}\left(\bm{\xi}_{s}\cdot\hat{\mathbf{n}}\right)\dot{J}\mathrm{d}A+\int_{\partial\mathcal{V}_{0}}\left(\bm{\xi}_{d}\cdot\hat{\mathbf{n}}\right)\dot{d}\mathrm{d}A+\int_{\partial\mathcal{V}_{0}}\left(\mathbb{S}\cdot\hat{\mathbf{n}}\right)\cdot\dot{\mathbf{u}}_{s}\mathrm{d}A,$$

(11)

where $\mathcal{V}_{0}$ and $\partial\mathcal{V}_{0}$ represent the material element and its boundary, respectively, and $\dot{\square}~{}\equiv~{}\frac{\partial}{\partial t}\square$. Grouping the surface integrals in Equation (11) and then using the divergence theorem leads to a local form of this inequality,

$$\dot{\psi}+\nabla_{0}\cdot\left(\rho_{g}\mu_{g}\mathbf{W}_{g}+\rho_{l}\mu_{l}\mathbf{W}_{l}-\dot{\varphi_{g}}\bm{\xi}_{g}-\dot{\varphi_{l}}\bm{\xi}_{l}-\dot{J}\bm{\xi}_{s}-\dot{d}\bm{\xi}_{d}-\mathbb{S}\cdot\dot{\mathbf{u}}_{s}\right)\leq 0.$$

(12)

We must also ensure that the no-void condition (Equation 3) is satisfied. This constraint can be incorporated into the above thermodynamic restriction through the use of a Lagrange multiplier, $p$, which is the thermodynamic mixture pressure. To make this constraint compatible with the dimensions of the dissipation inequality, we differentiate Equation (3) with respect to time and use Jacobi’s formula to evaluate the derivative of $J$, arriving at

$$\dot{J}=J\mathbb{F}^{-\intercal}\colon\dot{\mathbb{F}}=\dot{\varphi_{g}}+\dot{\varphi_{l}}.$$

(13)

We combine the dissipation inequality (Equation 12) with the incompressibility constraint (Equation 13) by multiplying the latter by $p$ and summing the two.

Finally, we elucidate the consequences of this inequality by asserting that $\psi$ can, in general, be a function of composition via $\varphi_{g}$, $\varphi_{l}$, $d$, and their gradients, and of deformation via $\mathbb{F}$ and $\nabla_{0}J$, so that $\psi~{}=~{}\psi\left(\varphi_{g},\varphi_{l},d,\mathbb{F},\nabla_{0}\varphi_{g},\nabla_{0}\varphi_{l},\nabla_{0}d,\nabla_{0}J\right)$. Note that we could equivalently write $\psi$ as a function of true volume fractions and their Eulerian gradients, but it is convenient to use nominal quantities as the independent variables here. We then use the chain rule as well as conservation of mass (Equation 8) and conservation of momentum (Equation 10) to arrive at

$$\begin{split}&\left(\frac{\partial\psi}{\partial\varphi_{g}}-\rho_{g}\mu_{g}+p-\nabla_{0}\cdot\bm{\xi}_{g}\right)\dot{\varphi_{g}}+\left(\frac{\partial\psi}{\partial\varphi_{l}}-\rho_{l}\mu_{l}+p-\nabla_{0}\cdot\bm{\xi}_{l}\right)\dot{\varphi_{l}}\\ &+\left(\frac{\partial\psi}{\partial\nabla_{0}\varphi_{g}}-\bm{\xi}_{g}\right)\nabla_{0}\dot{\varphi_{g}}+\left(\frac{\partial\psi}{\partial\nabla_{0}\varphi_{l}}-\bm{\xi}_{l}\right)\nabla_{0}\dot{\varphi_{l}}+\left(\frac{\partial\psi}{\partial d}-\nabla_{0}\cdot\bm{\xi}_{d}\right)\dot{d}\\ &+\left(\frac{\partial\psi}{\partial\nabla_{0}d}-\bm{\xi}_{d}\right)\nabla_{0}\dot{d}+\left(\frac{\partial\psi}{\partial\mathbb{F}}-\mathbb{S}-pJ\mathbb{F}^{-\intercal}-J\left(\nabla_{0}\cdot\bm{\xi}_{s}\right)\mathbb{F}^{-\intercal}\right)\colon\dot{\mathbb{F}}\\ &+\left(\frac{\partial\psi}{\partial\nabla_{0}J}-\bm{\xi}_{s}\right)\nabla_{0}\dot{J}+\mathbf{W}_{g}\cdot\nabla_{0}\left(\rho_{g}\mu_{g}\right)+\mathbf{W}_{l}\cdot\nabla_{0}\left(\rho_{l}\mu_{l}\right)\leq 0.\end{split}$$

(14)

Due to the mutual independence of the arguments of $\psi$, this inequality is only satisfied for all possible configurations if each of the bracketed terms vanishes independently. This requirement then leads to a set of thermodynamic constraints,

	$$\rho_{g}\mu_{g}=p+\frac{\partial\psi}{\partial\varphi_{g}}-\nabla_{0}\cdot\frac{\partial\psi}{\partial\nabla_{0}\varphi_{g}},$$		(15a)
	$$\rho_{l}\mu_{l}=p+\frac{\partial\psi}{\partial\varphi_{l}}-\nabla_{0}\cdot\frac{\partial\psi}{\partial\nabla_{0}\varphi_{l}},$$		(15b)
	$$\frac{\partial\psi}{\partial d}-\nabla_{0}\cdot\frac{\partial\psi}{\partial\nabla_{0}d}=0,$$		(15c)
	$$\mathbb{S}=-pJ\mathbb{F}^{-\intercal}+\frac{\partial\psi}{\partial\mathbb{F}}-J\left(\nabla_{0}\cdot\frac{\partial\psi}{\partial\nabla_{0}J}\right)\mathbb{F}^{-\intercal},$$		(15d)

and the additional requirement that

$$\mathbf{W}_{g}\cdot\nabla_{0}\left(\rho_{g}\mu_{g}\right)+\mathbf{W}_{l}\cdot\nabla_{0}\left(\rho_{l}\mu_{l}\right)\leq 0.$$

(16)

The latter expression can be satisfied by choosing an appropriate form for the fluxes $\mathbf{W}_{g}$ and $\mathbf{W}_{l}$, as discussed in §II.5 below.

The evolution equations derived in §II.2 drive the system toward an equilibrium state defined by the minima of $\psi$. We must next construct a function $\psi$ that captures the different physical processes at play in our system. We assume that $\psi$ is, in general, an additive function of each physical contribution,

$$\psi=\psi_{\mathrm{bulk}}+\psi_{\mathrm{mix}}+\psi_{\mathrm{elastic}}+{\psi_{\mathrm{damage}}}+\psi_{\mathrm{interface}},$$

(17)

where $\psi_{\mathrm{bulk}}$, $\psi_{\mathrm{mix}}$, $\psi_{\mathrm{elastic}}$, ${{\psi_{\mathrm{damage}}}}$ and $\psi_{\mathrm{interface}}$ are the energetic contributions from changing the mass of gas or liquid in the material element, gas-liquid-solid interactions within the material element, elasticity of the solid skeleton, mechanical damage to the solid skeleton, and interfacial effects, respectively. We write these contributions in terms of the true volume fractions and their Eulerian gradients in order to provide greater physical insight; conversion to nominal quantities is straightforward.

The bulk free energy, $\psi_{\mathrm{bulk}}$, is the free energy associated with the amount of gas and liquid in the material element and is given by

$$\psi_{\mathrm{bulk}}\left(\phi_{\alpha},J\right)=\left(\mu_{l}^{0}\rho_{l}\phi_{l}+\mu_{g}^{0}\rho_{g}\phi_{g}\right)J,$$

(18)

where the reference chemical potential $\mu^{0}_{\alpha}$ is the energy associated with adding one unit mass of fluid $\alpha$ to the material element, neglecting interactions between phases. For immiscible materials, $\mu^{0}_{\alpha}$ is a constant.

The mixing (interaction) energy, $\psi_{\mathrm{mix}}$, depends on the pore-scale arrangement of the phases and on the wetting characteristics of the fluid-fluid-solid system. To construct a free energy that gives the desired phenomenological behaviour for our system, we draw inspiration from Flory-Huggins theory used in polymer physics [42]. In a polymer solution, the Flory-Huggins energy, $\psi_{\mathrm{FH}}$, for a multiphase mixture is

$$\psi_{\mathrm{FH}}=k_{B}TJ\left[\sum_{\alpha\neq s}\Omega_{\alpha}\phi_{\alpha}\log\phi_{\alpha}+\sum_{\alpha\neq\beta}\hat{\chi}_{\alpha\beta}\phi_{\alpha}\phi_{\beta}\right],$$

(19)

where $k_{B}$ is the Boltzmann constant, $T$ is temperature, $\Omega_{\alpha}$ is inversely proportional to the characteristic size of particles of phase $\alpha$, and $\hat{\chi}_{\alpha\beta}$ characterises the strength of unfavourable interactions between phases $\alpha$ and $\beta$. The first sum in equation (19) represents the entropy of mixing of the different phases. Due to the large size of solid grains compared to fluid molecules, it is assumed that $\Omega_{s}\ll\Omega_{f}$, where the subscript $s$ indicates a general solid phase and $f$ a general fluid phase, and so the contribution of solid to the entropy term is neglected. The entropy term promotes the mixing of different phases, and also constrains the relevant volume fractions to remain between zero and one. The second sum in equation (19) represents the enthalpy of interactions between different phases, and results in the formation of double-well potentials which promote phase separation. The Flory-Huggins free energy thus captures the key features we expect from gas-liquid-solid interactions in a granular system, namely a double-well structure in which volume fractions are constrained to remain between zero and one. As such, we use $\psi_{\mathrm{FH}}$ as motivation for constructing an appropriate form for $\psi_{\mathrm{mix}}$.

In our system, the most important interaction is that the gas is non-wetting to the solid relative to the liquid. We enforce this requirement by assuming that gas-liquid and gas-solid interactions are much more energetically expensive than liquid-solid interactions, meaning that $\hat{\chi}_{ls}\ll\hat{\chi}_{gs,gl}$. We therefore neglect the liquid-solid interaction. In general, the remaining coefficients, $\hat{\chi}_{gs}$ and $\hat{\chi}_{gl}$, depend on physical quantities such as surface energies and pore structure; we take them to be constant material properties here. For simplicity, we further assume that $\Omega_{g}=\Omega_{l}=\Omega$. The mixing energy is then given by

$$\psi_{\mathrm{mix}}\left(\phi_{\alpha},J\right)=\mathcal{E}_{\mathrm{mix}}\left(\phi_{g}\log{\phi_{g}}+\phi_{l}\log{\phi_{l}}+\chi_{gs}\phi_{g}\phi_{s}+\chi_{gl}\phi_{g}\phi_{l}\right)J,$$

(20)

where $\chi_{\alpha\beta}=\frac{\hat{\chi}_{\alpha\beta}}{\Omega}$ and $\mathcal{E}_{\mathrm{mix}}$ is a characteristic energy density associated with mixing, which we treat as a material property, and into which we have absorbed $\Omega$. Figure 2 shows a ternary plot of this interaction energy. The key feature of the energetic landscape is that there are local minima near $\phi_{g}=1$ and along $\phi_{g}=0$, corresponding to open gas cavities and to a liquid-saturated porous matrix, respectively.

The specific values of $\chi_{gs}$ and $\chi_{gl}$ control the relative depth of the two minima, but not the overall structure of the energy landscape.

The elastic energy of the solid matrix due to deformation is given by $\psi_{\mathrm{elastic}}$, and is usually expressed in terms of a strain-energy density function, $\mathcal{W}_{\mathrm{el}}\left(\mathbb{F}\right)$. For simplicity, we use a standard neo-Hookean strain-energy density,

$$\mathcal{W}_{\mathrm{el}}=\frac{1}{2}\mathcal{G}\left[\textrm{tr}\left(\mathbb{F}^{\intercal}\mathbb{F}\right)-3-2\log{J}\right]+\frac{1}{2}\mathcal{K}\left(\log{J}\right)^{2},$$

(21)

where $\mathcal{G}$ and $\mathcal{K}$ are the ‘drained’ shear and bulk moduli of the solid skeleton, respectively. Note that the derivation that follows is valid for any choice of $\mathcal{W}_{\mathrm{el}}$.

A notable characteristic of a solid skeleton with a granular microstructure is that tension does not result in the storage of elastic energy, since grains only interact when in contact. In order to specify a different elastic response between open cavities and a continuous solid packing, we define an additional order parameter, $d\left(\mathbf{x},t\right)$, which we identify as a Bourdin-style damage parameter [43]. The damage parameter represents a phase-field approximation of the fractures to the solid skeleton that are caused by the opening of gas cavities. Locally, a value of $d=0$ corresponds to an undamaged solid skeleton, in which the solid grains are in contact and hence provide a standard elastic response to compression. In contrast, a value of $d=1$ corresponds to a fully damaged solid skeleton that provides no elastic resistance to further compression or tension. The use of a phase-field damage parameter ensures a smooth transition between open cavities and the rest of the porous medium

To implement damage, we assume that the strain-energy density can be decomposed into distinct tensile and compressive parts, denoted $\mathcal{W}^{+}$ and $\mathcal{W}^{-}$, respectively [40, 41]. We demonstrate a suitable decomposition of $\mathcal{W}_{\mathrm{el}}$ into tensile and compressive parts for a non-linear material in Appendix A, following the approach of Tang et al. [41]. We can then assume that the tensile part of the elastic energy is degraded by $d$ [40], such that,

	$$\psi_{\mathrm{elastic}}\left(d,\mathbb{F},J\right)=g\left(d\right)\mathcal{W}^{+}+\mathcal{W}^{-},$$		(22a)
	$$g\left(d\right)=\left(1-\Theta\right)\left(1-d\right)^{2}+\Theta,$$		(22b)

where $\Theta\ll 1$ is a numerical smoothing parameter. In the undamaged limit ($d=0$), this formulation reduces to $\psi_{\mathrm{elastic}}=\mathcal{W}^{+}+\mathcal{W}^{-}=\mathcal{W}_{\mathrm{el}}$, as expected.

According to Griffith’s theory of fracture, the energy required to create a fracture in a material is given by the critical fracture energy, $\mathcal{G}_{c}$, integrated over the area of the fracture [44]. In a granular medium, the ‘fracture’ energy is associated with the energy required to cause the decohesion of neighbouring solid grains. In our phase-field approach, the energy required to form such fractures is approximated by a bulk energy, $\psi_{\mathrm{damage}}$, [40, 43], of the form,

$$\psi_{\mathrm{damage}}\left(d,\nabla_{0}d\right)=\mathcal{G}_{c}\left(\frac{d^{2}}{4l_{d}}+l_{d}|\nabla_{0}d|^{2}\right),$$

(23)

where $l_{d}$ controls the width of the transition between damaged and undamaged material. For a non-cohesive solid skeleton, $\mathcal{G}_{c}$ is very small, with the only resistance to decohesion coming from wetting liquid bridges between solid grains. Note also that damage of a non-cohesive solid skeleton is reversible: if a cavity closes, $d$ will return to zero and the skeleton will return to standard elastic behaviour.

Finally, $\psi_{\mathrm{interface}}$ measures the energy of diffuse interfaces between different phases. In a phase-field formulation, interfaces are represented implicitly as regions with large gradients in composition. The interfacial energy thus depends on gradients in volume fraction [25, 45] and the total interfacial energy is the sum of the energy of each of the possible types of interface [36]. The simplest such dependence can in general be written as

$$\Upsilon_{\textrm{int}}=\sum_{\alpha}\frac{\gamma_{\alpha}}{2}|\nabla\phi_{\alpha}|^{2}J,$$

(24)

where $\gamma_{\alpha}$ are the associated interfacial coefficients. For a two-phase system (say, liquid and gas), the no-void condition allows elimination of one of the two volume fractions (say, $\phi_{g}$) such that $\Upsilon_{\textrm{int}}=\frac{\gamma}{2}|\nabla\phi_{l}|^{2}J$, where $\gamma=\gamma_{l}+\gamma_{g}$, and $\gamma$ can be directly related to the surface energy of the sharp interface formed between these two phases. For our three-phase system, we use the no-void condition to eliminate $\phi_{s}$ and write the total interfacial energy solely in terms of the gas and liquid fractions,

$$\psi_{\mathrm{interface}}\left(J,\nabla\phi_{\alpha}\right)=\frac{\gamma_{1}}{2}|\nabla\phi_{g}|^{2}J+\frac{\gamma_{2}}{2}|\nabla\phi_{l}|^{2}J+\frac{\gamma_{3}}{2}\left(\nabla\phi_{g}\cdot\nabla\phi_{l}\right)J,$$

(25)

where $\gamma_{1}=\gamma_{g}+\gamma_{s}$, $\gamma_{2}=\gamma_{l}+\gamma_{s}$ and $\gamma_{3}=2\gamma_{s}$. Unlike for a two-phase system, the relationship between these quantities and the associated surface energies is not clear.

The free energy constructed in §II.4 can now be combined with the thermodynamic constraints resulting from the energy inequality (Equations 15) to find expressions for the constitutive behaviour of the system in terms of the different driving mechanisms. We use Equations (15a) and (15b) to decompose the chemical potential for the gas and liquid phases as the sum of the pressure and a capillary potential, $\Psi_{\alpha}$,

	$\displaystyle\rho_{\alpha}\mu_{\alpha}$	$\displaystyle=p+\Psi_{\alpha},$		(26a)
	$\displaystyle\Psi_{\alpha}$	$\displaystyle=\frac{\partial\psi}{\partial\varphi_{\alpha}}-\nabla_{0}\cdot\left(\frac{\partial\psi}{\partial\nabla_{0}\varphi_{\alpha}}\right)\equiv\frac{\delta\psi}{\delta\varphi_{\alpha}}.$		(26b)

The capillary potentials represent the thermodynamic forcing on the gas and liquid phases due to the interactions between them. Note that $\Psi_{\alpha}$ can be viewed as the variational derivative of $\psi$ with respect to $\varphi_{\alpha}$.

Equation (15c) leads to an elliptic equation for $d$ as a function of the tensile elastic energy. For our specified free energy, the result is

$$4l_{d}\left(1-\Theta\right)\left(1-d\right)\mathcal{W}^{+}+4\mathcal{G}_{c}l_{d}^{2}\nabla_{0}^{2}d=\mathcal{G}_{c}d.$$

(27)

In the limit $\mathcal{G}_{c}\to 0$ (no energetic resistance to decohesion), Equation (27) reduces to $\mathcal{W}^{+}\left(1-d\right)=0$, such that any non-zero tension results in immediate fracture of the solid skeleton. This case essentially reproduces the sharp-interface limit of open cavities within a poroelastic continuum.

Equation (15d) allows the total stress to be decomposed into three distinct parts, converting from the first Piola-Kirchhoff stress to the Cauchy stress to arrive at

$$\mathbb{T}=-p\mathbb{I}+\bm{\sigma}^{\prime}+\mathbb{K}.$$

(28)

The effective (elastic) stress, $\bm{\sigma}^{\prime}$, is

$$\bm{\sigma}^{\prime}=\frac{1}{J}\frac{\partial\psi_{\mathrm{elastic}}}{\partial\mathbb{F}}\mathbb{F}^{\intercal}=g\left(d\right)\frac{1}{J}\frac{\partial\mathcal{W}^{+}}{\partial\mathbb{F}}\mathbb{F}^{\intercal}+\frac{1}{J}\frac{\partial\mathcal{W}^{-}}{\partial\mathbb{F}}\mathbb{F}^{\intercal},$$

(29)

where we have used $\psi_{\mathrm{elastic}}$ from §II.4 and $g\left(d\right)$ is the degradation function defined in Equation (22b). The Korteweg stress is

$$\mathbb{K}=\frac{1}{J}\frac{\partial\psi_{\mathrm{interface}}}{\partial\mathbb{F}}\mathbb{F}^{\intercal}-\left(\nabla_{0}\cdot\frac{\partial\psi_{\mathrm{interface}}}{\partial\nabla_{0}J}\right)\mathbb{I}.$$

(30)

The Korteweg stress is the bulk representation of interfacial tension at diffuse, internal interfaces and thus resists gradients in composition. The stress balance derived in Equation (9) allows us to link the effective and Korteweg stresses with gradients in the chemical potential via the pressure field.

Recall that Equation (16) constrains the relationship between fluid chemical potentials and fluxes. The simplest way to satisfy this constraint is for the fluxes to be proportional to gradients in chemical potential, such that,

$$\mathbf{W}_{\alpha}=-\mathbb{M}_{\alpha}\nabla_{0}\left(\rho_{\alpha}\mu_{\alpha}\right),$$

(31)

where $\mathbb{M}_{\alpha}\left(\varphi_{\alpha},\mathbb{F}\right)$ is a positive definite mobility tensor. The associated Eulerian volume fluxes are given by

$$\mathbf{w}_{\alpha}=-\mathcal{M}_{\alpha}\nabla\left(\rho_{\alpha}\mu_{\alpha}\right),$$

(32)

where $\mathcal{M}_{\alpha}\left(\phi_{\alpha}\right)$ is the Eulerian mobility of fluid $\alpha$. In the Eulerian frame, the pore space of granular packings is typically fairly isotropic, so we assume a scalar mobility for simplicity. At the least, $\mathcal{M}_{\alpha}$ must be a linear function of the fluid volume fraction, so that if there is no fluid present, then there will be no associated flux. In the interest of developing a simple model that captures the leading-order behaviour of the system, we choose $\mathcal{M}_{\alpha}=\phi_{\alpha}\mathcal{M}_{\alpha}^{0}$, where $\mathcal{M}_{\alpha}^{0}$ is a constant for a fluid $\alpha$ within a particular porous medium; however, more complicated mobility laws could also be used. For example, $\mathcal{M}_{\alpha}$ could be decomposed into a product of viscosity, absolute permeability, and relative permeability. These relations between the fluid fluxes and their respective chemical potentials complete our model.

Explicit expressions for the effective stress, the capillary potential and the Korteweg stress for the particular form of the free energy specified in §II.4 are given in Appendices A, B and C, respectively. A full set of governing equations is presented in Appendix D.

For uniaxial flow and deformation, the model derived above can be reduced to one dimension. Although uniaxial behaviour is a strong constraint, studying such a system is significantly easier, both conceptually and computationally, and allows us to gain substantial insight into what is ultimately a complex model. Uniaxial behaviour implies that

	$$\mathbf{v}_{\alpha}=v_{\alpha}\left(x,t\right)\mathbf{\hat{e}}_{x},$$		(33a)
	$$\mathbf{u}_{s}=u_{s}\left(x,t\right)\mathbf{\hat{e}}_{x},$$		(33b)
	$$\phi_{\alpha}=\phi_{\alpha}\left(x,t\right),$$		(33c)
	$$d=d\left(x,t\right),$$		(33d)

as illustrated in Figure 3.

At this point, it is convenient to non-dimensionalise our model. We do so using the following scalings: $x=L\tilde{x}$, $q=q_{0}\tilde{q}$, and $t=\tau\tilde{t}$, where $L$ is a characteristic length scale, $q_{0}$ is the mixture flux at the boundary, and $\tau=\frac{L^{2}}{\mathcal{M}^{0}_{l}\mathcal{E}_{\textrm{mix}}}$ is the time scale associated with the transport of liquid across a distance $L$ by capillary effects. This scaling motivates us to introduce the following dimensionless groups, which characterise the relative strengths of the different physical processes at play:

$\displaystyle Q_{0}=\frac{q_{0}\tau}{L},~{}~{}\mathcal{M}_{r}=\frac{\mathcal{M}^{0}_{g}}{\mathcal{M}^{0}_{l}},\quad\mathcal{D}=\frac{\mathcal{E}_{\mathrm{mix}}}{\mathcal{G}},\quad\Gamma_{i}=\frac{\gamma_{i}}{\mathcal{E}_{\mathrm{mix}}L^{2}},\quad\kappa_{c}=\frac{\mathcal{G}_{c}}{\mathcal{E}_{\mathrm{mix}}l_{d}},\quad\Lambda_{d}=\frac{l_{d}^{2}}{L^{2}}.$

The strength of background fluid flow relative to capillary forces is measured by $Q_{0}$, while $\mathcal{M}_{r}$ is the ratio of gas mobility to liquid mobility, $\mathcal{D}$ is the characteristic ‘deformability’ of the solid, $\Gamma_{i}$ ($i=1,2,3$) are the rescaled interfacial coefficients, $\kappa_{c}$ is the strength of cohesion between solid grains, and $\Lambda_{d}$ characterises the sharpness of the transition between damaged and undamaged material. We work exclusively with dimensionless quantities going forward, dropping the tildes for convenience.

In 1D, Equation (7c), which enforces the incompressibility of the mixture, simplifies to

$$\frac{\partial q}{\partial x}=0.$$

(34)

For the scalings defined above, the non-dimensional boundary condition for the flux is $q\left(0,t\right)=Q_{0}$. The solution to Equation (34) is therefore given by $q=Q_{0}$. For uniaxial flow and deformation, the deformation gradient tensor, $\mathbb{F}$, is diagonal, which greatly simplifies the stresses: only the $xx$-components of $\bm{\sigma}^{\prime}$ and $\mathbb{K}$, denoted ${\sigma}^{\prime}_{xx}$ and $K_{xx}$, respectively, appear explicitly in the 1D model. We use momentum conservation (Equation 9) and the stress decomposition (Equation 28) to link $p$ with $\sigma^{\prime}_{xx}$ and $K_{xx}$:

$$\frac{\partial p}{\partial x}=\frac{\partial}{\partial x}\left({\sigma}^{\prime}_{xx}+K_{xx}\right).$$

(35)

We now replace the fluid fluxes in the 1D evolution equation for $\phi_{g}$, given by Equation (7a), with the constitutive behaviour outlined above in Equations (26) and (32). The gas fraction then evolves according to

$$\begin{split}\frac{\partial\phi_{g}}{\partial t}=-Q_{0}\frac{\partial\phi_{g}}{\partial x}+\mathcal{M}_{r}&\frac{\partial}{\partial x}\left[\left(1-\phi_{g}\right)\phi_{g}\frac{\partial}{\partial x}\left({\Psi}_{g}+{\sigma}^{\prime}_{xx}+{K}_{xx}\right)\right]\\ -&\frac{\partial}{\partial x}\left[\phi_{g}\left(\phi-\phi_{g}\right)\frac{\partial}{\partial x}\left({\Psi}_{l}+{\sigma}^{\prime}_{xx}+{K}_{xx}\right)\right],\end{split}$$

(36)

where $\Psi_{\alpha}$ is now the dimensionless capillary potential, which is specified below. Similarly, we combine the 1D evolution equation for $\phi$ (Equation 7b) with the relevant constitutive behaviour to give

$$\begin{split}\frac{\partial\phi}{\partial t}=-Q_{0}\frac{\partial\phi}{\partial x}+\mathcal{M}_{r}&\frac{\partial}{\partial x}\left[\left(1-\phi\right)\phi_{g}\frac{\partial}{\partial x}\left({\Psi}_{g}+{\sigma}^{\prime}_{xx}+{K}_{xx}\right)\right]\\ +&\frac{\partial}{\partial x}\left[\left(1-\phi\right)\left(\phi-\phi_{g}\right)\frac{\partial}{\partial x}\left({\Psi}_{l}+{\sigma}^{\prime}_{xx}+{K}_{xx}\right)\right].\end{split}$$

(37)

We close our model with the 1D version of Equation (27), which describes the distribution of damage,

$$4\left(1-\Theta\right)\left(1-d\right){\mathcal{W}}^{+}+4\kappa_{c}\Lambda_{d}\frac{\partial^{2}d}{\partial X^{2}}=\kappa_{c}d,$$

(38)

where ${\mathcal{W}}^{+}$ is the dimensionless tensile energy, which in 1D is solely a function of $\phi$. The model is thus reduced to two non-linear evolution equations (for $\phi_{g}$ and for $\phi$) and one elliptic equation (for $d$).

By non-dimensionalising the capillary potentials (Equation 26b) and undertaking the required differentiation (Appendix B), we find that

	$${\Psi}_{g}={\Pi}_{g}-\Gamma_{1}\frac{\partial^{2}\phi_{g}}{\partial x^{2}}-\frac{\Gamma_{3}}{2}\frac{\partial^{2}\phi_{l}}{\partial x^{2}},$$		(39a)
	$${\Psi}_{l}={\Pi}_{l}-\Gamma_{2}\frac{\partial^{2}\phi_{l}}{\partial x^{2}}-\frac{\Gamma_{3}}{2}\frac{\partial^{2}\phi_{g}}{\partial x^{2}},$$		(39b)

where the interaction potentials, ${\Pi}_{g}$ and ${\Pi}_{l}$, are

	$${\Pi}_{g}=\log{\phi_{g}}+\phi_{s}+\left(1-\phi_{g}\right)\left(\chi_{gl}\phi_{l}+\chi_{gs}\phi_{s}\right)$$		(40a)
and
	$${\Pi}_{l}=\log{\phi_{l}}+\phi_{s}+\phi_{g}\left(\chi_{gl}\phi_{g}-\chi_{gs}\phi_{s}\right).$$		(40b)

The interaction potentials result from the thermodynamic interactions between the phases, as described by $\psi_{\textrm{mix}}$, and are the drivers for phase separation.

For uniaxial deformation, tension can be defined solely volumetrically: the material is in tension if $J>1$ and in compression if $J\leq 1$. The non-dimensional effective stress (Equation 29) thus reduces to

$${\sigma}^{\prime}_{xx}=\begin{cases}{\sigma}^{\prime}_{0}\left(\phi\right)&J\leq 1,\\ {\sigma}^{\prime}_{0}\left(\phi\right)g\left(d\right)&J>1,\end{cases}$$

(41)

where ${\sigma}^{\prime}_{0}$ is the undamaged neo-Hookean stress given by

$${\sigma}^{\prime}_{0}=\frac{1}{\mathcal{D}}\left[J-\frac{1}{J}+\varsigma\frac{\log{J}}{J}\right],$$

(42)

and $\varsigma=\frac{\mathcal{K}}{\mathcal{G}}$. The derivation of this effective stress can be found in Appendix A. Similarly, the tensile energy used in Equation (38) is given by,

$${\mathcal{W}}^{+}=\begin{cases}0&\textrm{for}~{}J\leq 1,\\ {\mathcal{W}}_{el}^{0}&\textrm{for}~{}J>1,\end{cases}$$

(43)

where ${\mathcal{W}}_{el}^{0}=\frac{1}{2\mathcal{D}}\left[J^{2}-1-2\log{J}+\varsigma\left(\log{J}\right)^{2}\right]$. The incompressibility condition (Equation 4) links $J$ to the porosity field:

$$J=\frac{1-\phi_{0}}{1-\phi}.$$

(44)

Finally, $K_{xx}$ is derived in Appendix C and is given by

$$\begin{split}{K}_{xx}=\Gamma_{1}\left[\phi_{g}\frac{\partial^{2}\phi_{g}}{\partial x^{2}}-\frac{1}{2}\left(\frac{\partial\phi_{g}}{\partial x}\right)^{2}\right]&+\Gamma_{2}\left[\phi_{l}\frac{\partial^{2}\phi_{l}}{\partial x^{2}}-\frac{1}{2}\left(\frac{\partial\phi_{l}}{\partial x}\right)^{2}\right]\\ &+\frac{1}{2}\Gamma_{3}\left[\phi_{g}\frac{\partial^{2}\phi_{l}}{\partial x^{2}}+\phi_{l}\frac{\partial^{2}\phi_{g}}{\partial x^{2}}-\frac{\partial\phi_{l}}{\partial x}\frac{\partial\phi_{g}}{\partial x}\right].\end{split}$$

(45)

Equations (36)-(45) comprise a complete model for our system in 1D, and can be solved numerically given appropriate initial and boundary conditions. We do so below by discretising in space using finite differences on a staggered grid and integrating in time using MATLAB’s built-in solver for stiff differential equations, ode15s. For the remainder of this paper, we limit ourselves to the 1D case. We do not anticipate fundamentally different physical behaviour in 2D or 3D.

In the no-solid limit, $\phi\equiv\phi_{0}\equiv 1$, our model reduces to a gas-liquid system in which the two phases separate in an unconstrained domain (no porous medium), as illustrated in Figure 3a. This scenario is typically described by the Cahn-Hilliard equation for binary mixtures [25]. In this limit, the evolution equation for the porosity (Equation 37) and the elliptic equation for the damage (Equation 38) degenerate. The evolution equation for the gas fraction (Equation 36) simplifies greatly, most notably due to disappearance of the elastic terms, becoming

$$\frac{\partial\phi_{g}}{\partial t}-\mathcal{M}_{r}\frac{\partial}{\partial x}\left[\phi_{g}\frac{\partial\mu_{g}}{\partial x}\right]=0.$$

(46)

The gas chemical potential is now given by

$$\mu_{g}=\log{\phi_{g}}+\chi_{gl}\left(1-\phi_{g}\right)^{2}-\Gamma\left[\frac{1}{2}\left(\frac{\partial\phi_{g}}{\partial x}\right)^{2}+\left(1-\phi_{g}\right)\frac{\partial^{2}\phi_{g}}{\partial x^{2}}\right],$$

(47)

where we have defined $\Gamma=\Gamma_{1}+\Gamma_{2}-\Gamma_{3}=\Gamma_{g}+\Gamma_{l}$ as the characteristic interfacial coefficient between the fluid phases. Equations (46) and (47) resemble the standard Cahn-Hilliard equation for fluid-fluid phase separation [e.g., 31] but, in our case, also account for the mechanical stress generated across diffuse interfaces.

We now consider a test problem with periodic boundary conditions. For the initial condition, we use a homogeneous gas fraction superimposed with small, random perturbations. The subsequent evolution is standard fluid-fluid phase separation, also known as spinodal decomposition, in which the mixture separates spontaneously into distinct gas-rich and liquid-rich domains separated by diffuse interfaces (Figure 4).

Over time, these domains coarsen as smaller gas ‘bubbles’ merge with larger gas ‘bubbles’. Coarsening is driven by the evolution of the system toward minimum global energy (the smaller the number of bubbles, the smaller the total interfacial energy). In 1D, coarsening is very slow [46]; as such, we stop our simulations before reaching the eventual equilibrium of a single bubble.

In the no-gas limit, our model reduces to a diffuse-interface formulation of the well-known theory of single-phase large-deformation poroelasticity [47, 48]. In this case, there will be no thermodynamically driven phase separation because there is minimal energetic cost to liquid remaining within the pore space of the solid skeleton. However, flow can still drive phase separation, as demonstrated below. Setting $\phi_{g}\equiv 0$, Equation (36) for the evolution of the gas fraction degenerates and Equation (37) for the evolution of the porosity simplifies to

$$\frac{\partial\phi}{\partial t}+Q_{0}\frac{\partial\phi}{\partial x}-\frac{\partial}{\partial x}\left[\left(1-\phi\right)\mathcal{M}\left(\phi\right)\frac{\partial}{\partial x}\left({\sigma}^{\prime}_{xx}+{\Psi}\right)\right]=0,$$

(48)

with

$${\Psi}=\log{\phi}+1-\phi-\Gamma_{2}\left[\left(1-\phi\right)\frac{\partial^{2}\phi}{\partial x^{2}}+\frac{1}{2}\left(\frac{\partial\phi}{\partial x}\right)^{2}\right],$$

(49)

and a dimensionless mobility given by $\mathcal{M}\left(\phi\right)=\phi$. This equation is reminiscent of standard sharp-interface poroelasticity [e.g., 49], but with the addition of a potential, ${\Psi}$, that allows the formation of internal diffuse interfaces, such as those resulting from the fluid-driven opening of solid-free cavities within the porous skeleton. The effective stress here also depends on $d$ via Equation (38).

For a quantitative comparison of our model with sharp-interface poroelasticity, we consider the uniaxial deformation of a packing of soft grains due to net fluid flow from left to right (Figure 3b). We impose a constant liquid influx at the left boundary of non-dimensional flow rate $Q_{0}$. We take the right boundary to be permeable, allowing free outflow of the liquid but not of the solid. The result is that the packing will be compressed in the direction of the flow, with the left edge acting as an internal free boundary that moves downstream away from its initial position. Traditionally, this problem would be approached by solving the equations of large-deformation poroelasticity within the packing while explicitly tracking the position of the left edge of the packing as a moving boundary, assuming that the effective stress vanishes there. In our model, the left edge of the packing is an internal diffuse interface between a solid-free cavity ($\phi\sim 1$) and a porous domain ($0<\phi<1$).

For this problem, the sharp-interface formulation can be solved analytically at steady state, thus allowing for direct comparison with our phase-field approach. Following MacMinn et al. [49], we define the left edge of the packing to have a position $\delta\left(t\right)$, with a fluid-filled domain for $x<\delta$. According to the theory of large-deformation poroelasticity, the porosity in the domain $x~{}\in~{}[\delta\left(t\right),1]$ then evolves according to

$$\frac{\partial\phi}{\partial t}+Q_{0}\frac{\partial\phi}{\partial x}-\frac{\partial}{\partial x}\left[\left(1-\phi\right)\mathcal{M}\left(\phi\right)\frac{\partial{\sigma}^{\prime}_{0}}{\partial x}\right]=0.$$

(50)

We anticipate that our phase-field model would reduce to this formulation in the sharp-interface limit. Conservation of mass requires that the influx of liquid at the left boundary, $x=\delta$, must match the outflux at the right boundary, $x=1$. We use this fact, along with the requirement that the solid velocity vanishes at $x=1$, to enforce that the relative liquid flux at the right boundary is equal to the total mixture flux, $w_{l}=-\mathcal{M}\frac{\partial{\sigma}^{\prime}_{0}}{\partial x}=Q_{0}$. We also assume that the solid packing is relaxed at its left edge, so that the effective stress vanishes at $x=\delta$. The boundary conditions for Equation (50) are then given by

$${\sigma}^{\prime}_{0}\left(x=\delta\right)=0,\qquad\phi\frac{\partial{\sigma}^{\prime}_{0}}{\partial x}\Bigr{|}_{x=1}=-Q_{0}.$$

(51)

As discussed above, the steady-state porosity is piecewise:

$$\phi\left(x\right)=\begin{cases}1&x<\delta,\\ f\left(x\right)&\delta\leq x\leq 1.\end{cases}$$

(52)

Solving Equation (50) subject to the boundary conditions above (Equations 51) gives $f(x)$. An analytical solution for $f(x)$ is derived in Appendix E and plotted in Figure 5.

To compare our phase-field model to the sharp-interface result, we solve Equation (38) and Equation (48) numerically until reaching steady state. With regard to boundary conditions for these equations, we assume that the solid velocity vanishes at each end of the domain, so that the relative liquid flux, $w_{l}$, is equal to $Q_{0}$ at each end. To enforce the fact that the solution should not depend on anything exterior to the domain, we also impose that the gradients of porosity and damage vanish at both ends. The boundary conditions for Equation (48) are thus

$$w_{l}\left(x=0,1\right)=Q_{0},\quad\frac{\partial\phi}{\partial x}\Bigr{|}_{x=0,1}=0,\quad\frac{\partial d}{\partial x}\Bigr{|}_{x=0,1}=0.$$

(53)

With these boundary conditions, we solve Equations (38) and (48) numerically with initial condition $\phi\left(x,0\right)=\phi_{0}$.

Figure 5 compares the steady state of the phase-field simulations for various values of $\Gamma_{2}$ (solid lines), with the analytical solution to the sharp-interface model given in Appendix E (dashed line). In both cases, the solution comprises a solid-free region at the left and a compressed porous packing at the right. Within the latter, the porosity is non-linear with position. Away from the interface, the two solutions are in excellent qualitative and quantitative agreement. Near the interface, the discontinuity in the sharp-interface solution is approximated by a smooth interpolation in the phase-field model, the width of which is controlled by $\Gamma_{2}$. As $\Gamma_{2}$ decreases, the phase-field simulations converge toward the sharp-interface solution (Figure 5, inset).

To focus on a porous medium that is sufficiently large that boundary effects are negligible in the bulk of the system, we impose periodic boundary conditions. We also set $Q_{0}=0$ without loss of generality; a non-zero value of $Q_{0}$ corresponds to bulk translation at a fixed rate. For initial conditions, we set

	$\displaystyle\phi_{g}(x,0)$	$\displaystyle=S_{0}\phi_{0}+\eta(x),$		(54)
	$\displaystyle\phi(x,0)$	$\displaystyle=\phi_{0}-|\phi^{*}|,$		(55)

where $S_{0}$ is the initial gas saturation, $\eta$ is a field of small, random fluctuations, and $\phi^{*}\ll 1$ ensures an initial state of slight compression, so that the material is initially undamaged.

Evolving the system from this initial state, we see that perturbations in the gas fraction grow over time, deforming the solid skeleton (Figure 6).

As gas bubbles outgrow the available pore space, they expand to form cavities in the solid by displacing solid grains. As the cavities form and then grow, the rest of the solid skeleton is increasingly compressed into a smaller region. A local equilibrium state is reached once elastic stresses balance the thermodynamic forcing of the phase separation. In this equilibrium state, the compressed solid packing has a uniform volume fraction throughout the porous medium. The formation of cavities damages the solid skeleton such that there is negligible elastic stress within the cavities.

We investigate the impact of the deformability of the porous medium on phase separation by running numerical simulations over a range of values of $\mathcal{D}$ (Figure 7).

For a stiff porous medium ($\mathcal{D}\lesssim 1$), the elastic forces within the solid skeleton are too strong for the gas to deform the packing. The gas thus remains within the pore space and phase separation is suppressed. Conversely, for a sufficiently soft porous medium ($\mathcal{D}\gtrsim 1$), gas is able to open macroscopic cavities within the packing as described above. We investigate the impact of other parameters (such as $\chi_{gs}$ and $S_{0}$) in §IV.2 below.

To identify the parameters that control the onset of phase separation, we now perform a linear stability analysis for our model. We linearise the system about a base state that represents an undamaged, precompressed porous matrix of porosity $\phi_{0}^{\prime}=\phi_{0}~{}-~{}|\phi^{*}|$, whose pore space is occupied by a homogeneous distribution of gas with saturation $S_{0}$ and volume fraction $\phi_{g0}=S_{0}\phi_{0}$. We assume that perturbations of size $\epsilon\ll 1$ about this base state take the form of normal modes with growth rate $s$ and wavenumber $k$, such that the perturbed gas fraction and porosity are given by,

	$\displaystyle\phi_{g}$	$\displaystyle=\phi_{g0}+\epsilon\widehat{\phi}_{g}~{}\mathrm{exp}\left(st+ikx\right)+\mathcal{O}\left(\epsilon^{2}\right),$		(56a)
	$\displaystyle\phi$	$\displaystyle=\phi_{0}^{\prime}+\epsilon\widehat{\phi}~{}\mathrm{exp}\left(st+ikx\right)+\mathcal{O}\left(\epsilon^{2}\right),$		(56b)

where $\widehat{\phi}_{g}$ and $\widehat{\phi}$ characterise the $\mathcal{O}(\epsilon)$ contributions to $\phi_{g}$ and $\phi$, respectively.

Substituting this ansatz (Equations 56) into our model (Equations 36–45) and linearising in $\epsilon$ leads to a set of equations that describe the evolution of small perturbations. To $\mathcal{O}(\epsilon)$, our model thus reduces to the following algebraic equations for $\widehat{\phi}_{g}$ and $\widehat{\phi}$,

	$$s\widehat{\phi}_{g}=-\mathcal{M}_{r}\left(1-\phi_{g0}\right)\phi_{g0}\left[\left(h_{gg}k^{2}+\nu_{gg}k^{4}\right)\widehat{\phi}_{g}+\left(h_{gf}k^{2}+\nu_{gf}k^{4}\right)\widehat{\phi}\right]\\ +\left(\phi^{\prime}_{0}-\phi_{g0}\right)\phi_{g0}\left[\left(h_{lg}k^{2}+\nu_{lg}k^{4}\right)\widehat{\phi}_{g}+\left(h_{lf}k^{2}+\nu_{lf}k^{4}\right)\widehat{\phi}\right],$$		(57a)
	$$s\widehat{\phi}=-\mathcal{M}_{r}\left(1-\phi^{\prime}_{0}\right)\phi_{g0}\left[\left(h_{gg}k^{2}+\nu_{gg}k^{4}\right)\widehat{\phi}_{g}+\left(h_{gf}k^{2}+\nu_{gf}k^{4}\right)\widehat{\phi}\right]\\ -\left(1-\phi^{\prime}_{0}\right)\left(\phi^{\prime}_{0}-\phi_{g0}\right)\left[\left(h_{lg}k^{2}+\nu_{lg}k^{4}\right)\widehat{\phi}_{g}+\left(h_{lf}k^{2}+\nu_{lf}k^{4}\right)\widehat{\phi}\right].$$		(57b)

In the above, we have introduced the functions $h_{\alpha\beta}=\frac{\partial h_{\alpha}}{\partial\phi_{\beta}}|_{\phi^{\prime}_{0},\phi_{g0}}$, where $h_{\alpha}$ is the homogeneous part of the chemical potential of phase $\alpha$, and the parameters $\nu_{\alpha\beta}$, which are different combinations of interfacial coefficients. In general, $h_{\alpha\beta}$ and $\nu_{\alpha\beta}$ depend on the dimensionless groups defined in §II.6 and the base states $\phi_{g0}$ and $\phi^{\prime}_{0}$, but they are constant with respect to $k$ and $s$. Explicit expressions for $h_{\alpha\beta}$ and $\nu_{\alpha\beta}$ are presented in Appendix F. Collecting like terms and finding the eigenvalues of Equations (57) leads to the dispersion relation

$$s^{2}+\zeta_{1}\left(k\right)k^{2}s+\zeta_{2}\left(k\right)k^{4}=0,$$

(58)

where $\zeta_{1}\left(k\right)=a_{1}+a_{2}k^{2}$ and $\zeta_{2}\left(k\right)=b_{1}+b_{2}k^{2}+b_{3}k^{4}$. Explicit forms of $a_{i}$ and $b_{i}$ are given in Appendix F.