Fluid-fluid phase separation in a soft porous medium

Various biological and chemical processes lead to the nucleation and growth of non-wetting fluid bubbles within the pore space of a granular medium, such as the formation of gas bubbles in liquid-saturated lake-bed sediments. In sufficiently soft porous materials, the non-wetting nature of these bubbles can result in the formation of open cavities within the granular solid skeleton. Here, we consider this process through the lens of phase separation, where thermomechanics govern the separation of the non-wetting phase from a fluid-fluid-solid mixture. We construct a phase-field model informed by large-deformation poromechanics, in which two immiscible fluids interact with a poroelastic solid skeleton. Our model captures the competing effects of elasticity and fluid-fluid-solid interactions. We use a phase-field damage model to capture the mechanics of the granular solid. As a model problem, we consider an initial distribution of non-wetting fluid in the pore space that separates into multiple cavities. We use simulations and linear-stability analysis to identify the key parameters that control phase separation, the conditions that favour the formation of cavities, and the characteristic size of the resulting cavities.


I. INTRODUCTION
Multiphase flow through rigid porous materials has been studied extensively, but multiphase flow through deformable porous materials remains relatively poorly understood [1].
These systems are challenging because flows through deformable porous media inherently induce a strong two-way coupling: flow can deform the pore structure, which in turn affects the flow. Recent works have explored this coupling by considering the injection of an immiscible fluid into a saturated, deformable porous medium [2][3][4][5], highlighting how capillary forces can deform the structure of the pore space. Similar phenomena have also attracted attention from the perspective of pattern formation in multiphase frictional flows [6,7].
Deformation enables the emergence of completely new flow phenomena that do not have a rigid analogue, the most striking of which is the ability of the non-wetting phase to form open (solid-free) pathways and cavities [8,9]. Non-wetting cavities form due to capillary forces, which make it energetically costly for the non-wetting phase to invade narrow pore throats. If the solid skeleton is sufficiently soft, it may be energetically favourable for the non-wetting phase to instead displace solid and form macroscopic cavities within the porous skeleton [8,9]. The non-wetting phase can thus occupy one of two distinct states within a soft porous medium; either (1) invading the pore space by displacing the wetting phase or (2) forming open cavities by additionally displacing the solid (Figure 1). In a soft porous medium, the latter is resisted by elasticity and the competition between capillarity and elasticity determines which of these states is energetically favourable. The pore-scale physics of cavity formation are relatively well understood [10,11], but a continuum-scale model for this phenomenon is currently lacking. A quantitative understanding of the continuum-scale physics of cavity formation would have application to a variety of natural and industrial systems, including, for example, predicting the rate of gas venting from lake beds, sea beds and waste ponds [12][13][14].
The spontaneous formation of non-wetting cavities within a soft porous medium can be conceptualised as a process of thermodynamic phase separation. Recently, the problem of liquid-liquid phase separation within an elastic network has attracted substantial attention from the perspective of droplet formation within living cells [15][16][17][18][19][20]. Elastically modulated phase separation has also been studied within the context of hydrogels [21][22][23][24]. It is now well understood that the presence of an elastic network fundamentally alters the phase-separation process by introducing an energetic cost to the rearrangements that are needed for phase separation to occur. This energetic cost limits the final size of phase-separated domains [19], decreases the rate of post-separation coarsening [21], and in some cases can suppress phase separation entirely [19,20]. Fluid-fluid phase separation in a poroelastic medium has many similarities to these systems, as demonstrated and discussed below.
Phase-field models provide a powerful framework for modelling phase-separation processes. These models were first developed within the materials-science community to describe spinodal decomposition during the solidification of alloys [25][26][27]. The resulting Cahn-Hilliard equation has since been generalised for use in a range of problems in fluid dynamics [e.g., [28][29][30][31], including multiphase flow in rigid porous media [32,33]. Many Cahn-Hilliardstyle models focus on two phases, either one fluid and one solid [e.g., 23], two fluids [e.g., 34], or two solids [e.g., 35], but these models have also been generalised to an arbitrary number of phases [36]. Although most Cahn-Hilliard-style models are derived via an ad-hoc variational approach, Gurtin [37][38][39] later developed a thermodynamically consistent framework for these models by introducing the concept of subscale 'microforces' associated with changes in local composition. Gurtin's framework also allows for the natural inclusion of a deformable solid phase via conservative (elastic) and/or dissipative (plastic) contributions to the free energy [39], which are otherwise challenging to capture consistently in a traditional Cahn-Hilliard model.
Here, we develop a continuum model that captures the fundamental mechanics of the formation of non-wetting cavities within a soft porous medium. By conceptualising the transition between pore invasion and cavity formation as a phase-separation process, our phase-field model captures the competing physical processes in a thermodynamically consistent manner. In §II, we derive our model by extending the work of Hennessy et al. [23] to include an additional fluid phase and to allow for elastic degradation of the solid skeleton via a damage model [40,41]. Our derivation follows the approach of Gurtin [39], using balance laws to set out an energy inequality that imposes certain restrictions on the allowable constitutive behaviours. In §III, we simplify our model to one spatial dimension and benchmark it against known solutions for two different two-phase limits. In §IV, we consider a nearly uniform initial distribution of gas within the pore space of a slightly compressed solid skeleton, which may or may not lead to the spontaneous formation of multiple gas cavities. We then study the onset of phase separation through numerical simulations ( §IV A) and linear stability analysis ( §IV B). In §IV C, we investigate the characteristic size of the resulting cavities. Finally, in §V, we summarise our findings and highlight areas for future work.

II. MODEL DEVELOPMENT
We begin by considering a mixture of three phases, two fluids and a solid. At the pore scale, these phases exist within separate domains separated by sharp interfaces. At the continuum (Darcy) scale, however, the phases are mixed and the proportion of each at a particular point in space and time is measured by its volume fraction. The volume fractions thus act as phase-field order parameters that can be used to distinguish the physical characteristics of different regions, and to smoothly interpolate between them. Relative motion of the phases is governed by a thermodynamic framework that drives the system toward minima of its total free energy, as described in detail below.
We take the solid phase to be a porous skeleton comprising non-cohesive, incompressible solid 'grains'. The pore space between the grains is occupied by two immiscible, incompressible fluid phases. For clarity and without loss of generality, we refer to these fluids as 'liquid' and 'gas', where the gas is the non-wetting phase. The key feature of our model is that, under certain conditions, it is energetically favourable for the gas to separate from the other two phases by forming open cavities.
We next formulate our model. We begin with the kinematics of the mixture ( §II A).
We then formulate balance laws for mass and momentum ( §II B) and consider the thermodynamics of the mixture ( §II C). We formulate a free energy that captures the energetic characteristics of our system (elasticity, damage, and gas-liquid-solid interactions; §II D). We close the model by using thermodynamic arguments to determine constitutive relationships that link the various quantities in our balance laws ( §II E). Finally, we simplify our model to the case of a one dimensional (1D) system ( §II F).

A. Kinematics
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 X refers to the reference state and the Eulerian coordinate x refers to the current state. We denote spatial gradients with respect to the Lagrangian and Eulerian coordinates by ∇ 0 and ∇, respectively. The Lagrangian displacement of the solid is then U s (X, t) = x (X, t) − X and the Eulerian displacement is u s (x, t) = x − X (x, t). The deformation gradient tensor, F, describes the deformation of the solid skeleton relative to its reference state, Note that we use the convention (∇A) ij = ∂A i ∂x j for a general vector field A = A iêi and (∇ · B) i = ∂B ij ∂x j for a general tensor field B = B ijêiêj , whereê i andê j are Cartesian unit vectors and where we have adopted the Einstein summation convention. Using the definitions above, F is related to U s and u s via where 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, ϕ α , measure the volume of phase α per unit reference bulk volume, whereas the true volume fractions, φ α , measure the volume of phase α per unit current bulk volume, where α = s, l, g indicates the solid, liquid, and gas phases, respectively. These volume fractions are related by ϕ α = φ α J, where the Jacobian determinant, J = det(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, 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 ϕ s ≡ ϕ 0 is the solid fraction in the reference state, which we take to be uniform for simplicity. As a result, It is also convenient to define the true porosity, φ = φ g + φ l = 1 − φ s , and the relaxed porosity, φ 0 = 1 − φ 0 s . Below, we develop evolution equations for φ g and φ and then evaluate φ l = φ − φ g and φ s = 1 − φ as needed.

B. Balance Laws
We formulate evolution equations for φ g and φ by first considering conservation of mass for each phase. In the Eulerian frame, conservation of mass can be written where ρ α and v α are the density and local velocity of phase α, respectively. We assume for simplicity that all three phases are incompressible, meaning that ρ α is a constant. The Eulerian volume flux, w α , of phase α relative to the solid skeleton is We reduce Equations (5) to two evolution equations, one each for φ g and φ, and an elliptic equation for the total mixture flux, q = α φ α v α = v s + w g + w l , by eliminating the phase velocities in favour of the relative fluxes. The result of these manipulations is We derive thermodynamically consistent constitutive expressions for the fluxes w α in §II E.
We express mass conservation in Lagrangian form by using the Reynolds transport theorem to rewrite Equation (5) as where the nominal volume flux, W α , is related to the true flux by W α = JF −1 w α .
Further restrictions can be imposed on our system by considering conservation of momentum. The true (Cauchy) total stress, 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 T must be divergence free, In the Lagrangian frame, momentum balance is most naturally expressed in terms of the first Piola-Kirchhoff stress, S, where S = JTF − measures the total internal force per unit reference area within the mixture.

C. Thermodynamics
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, ψ, of a bulk material element in the Lagrangian frame. The rate of change of energy per unit mass of phase α due to a net flux of phase α is measured by the chemical potential, µ α .
The key ingredient in Gurtin's approach is that changes in energy resulting from changes in composition, as measured by changes in ϕ α , can be represented by the action of a nominal vector 'microstress', ξ α , 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, ξ d (see §II D below). The working of external forces is represented by the action of the total stress, S, on the boundary of the material element.
Summing these different contributions leads to a dissipation inequality, where V 0 and ∂V 0 represent the material element and its boundary, respectively, anḋ ≡ ∂ ∂t . Grouping the surface integrals in Equation (11) and then using the divergence theorem leads to a local form of this inequality, 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 aṫ 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 ψ can, in general, be a function of composition via ϕ g , ϕ l , d, and their gradients, and of deformation Note that we could equivalently write ψ 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 Due to the mutual independence of the arguments of ψ, 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, and the additional requirement that The latter expression can be satisfied by choosing an appropriate form for the fluxes W g and W l , as discussed in §II E below.

D. Free Energy
The evolution equations derived in §II B drive the system toward an equilibrium state defined by the minima of ψ. We must next construct a function ψ that captures the different physical processes at play in our system. We assume that ψ is, in general, an additive function of each physical contribution, where ψ bulk , ψ mix , ψ elastic , ψ damage and ψ 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, ψ bulk , is the free energy associated with the amount of gas and liquid in the material element and is given by where the reference chemical potential µ 0 α is the energy associated with adding one unit mass of fluid α to the material element, neglecting interactions between phases. For immiscible materials, µ 0 α is a constant. The mixing (interaction) energy, ψ 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, ψ FH , for a multiphase mixture is where k B is the Boltzmann constant, T is temperature, Ω α is inversely proportional to the characteristic size of particles of phase α, andχ αβ characterises the strength of unfavourable interactions between phases α and β. 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 Ω s Ω 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 ψ FH as motivation for constructing an appropriate form for ψ 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χ ls χ gs,gl . We therefore neglect the liquid-solid interaction. In general, the remaining coefficients,χ gs andχ 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 Ω g = Ω l = Ω. The mixing energy is then given by where χ αβ =χ αβ Ω and E mix is a characteristic energy density associated with mixing, which we treat as a material property, and into which we have absorbed Ω. Figure 2 shows a ternary plot of this interaction energy. The key feature of the energetic landscape is that there are local minima near φ g = 1 and along φ g = 0, corresponding to open gas cavities and to a liquid-saturated porous matrix, respectively. The specific values of χ gs and χ 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 ψ elastic , and is usually expressed in terms of a strain-energy density function, W el (F). For simplicity, we use a standard neo-Hookean strain-energy density, where G and 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 W 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 (x, t), which we identify as a Bourdinstyle damage parameter [43]. The damage parameter represents a phase-field approximation To implement damage, we assume that the strain-energy density can be decomposed into distinct tensile and compressive parts, denoted W + and W − , respectively [40,41]. We demonstrate a suitable decomposition of W el into tensile and compressive parts for a nonlinear 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, where Θ 1 is a numerical smoothing parameter. In the undamaged limit (d = 0), this formulation reduces to ψ elastic = W + + W − = W 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, 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, ψ damage , [40,43], of the form, where l d controls the width of the transition between damaged and undamaged material.
For a non-cohesive solid skeleton, 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 noncohesive solid skeleton is reversible: if a cavity closes, d will return to zero and the skeleton will return to standard elastic behaviour.
Finally, ψ 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 where γ α 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, φ g ) such that Υ int = γ 2 |∇φ l | 2 J, where γ = γ l + γ g , and γ 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 φ s and write the total interfacial energy solely in terms of the gas and liquid fractions, where γ 1 = γ g +γ s , γ 2 = γ l +γ s and γ 3 = 2γ s . Unlike for a two-phase system, the relationship between these quantities and the associated surface energies is not clear.

E. Constitutive Behaviour
The free energy constructed in §II D 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, Ψ α , The capillary potentials represent the thermodynamic forcing on the gas and liquid phases due to the interactions between them. Note that Ψ α can be viewed as the variational derivative of ψ with respect to ϕ α .
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 In the limit G c → 0 (no energetic resistance to decohesion), Equation (27)  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 The effective (elastic) stress, σ , is where we have used ψ elastic from §II D and g (d) is the degradation function defined in Equation (22b). The Korteweg stress is 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, where M α (ϕ α , F) is a positive definite mobility tensor. The associated Eulerian volume fluxes are given by where M α (φ α ) is the Eulerian mobility of fluid α. 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, M α 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 M α = φ α M 0 α , where M 0 α is a constant for a fluid α within a particular porous medium; however, more complicated mobility laws could also be used. For example, M α 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 D are given in Appendices A, B and C, respectively. A full set of governing equations is presented in Appendix D.

F. Reduction to 1D
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 as illustrated in Figure 3.
At this point, it is convenient to non-dimensionalise our model. We do so using the following scalings: x = Lx, q = q 0q , and t = τt, where L is a characteristic length scale, q 0 is the mixture flux at the boundary, and τ = L 2 M 0 l E 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: The strength of background fluid flow relative to capillary forces is measured by Q 0 , while M r is the ratio of gas mobility to liquid mobility, D is the characteristic 'deformability' of the solid, Γ i (i = 1, 2, 3) are the rescaled interfacial coefficients, κ c is the strength of cohesion between solid grains, and Λ 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 For the scalings defined above, the non-dimensional boundary condition for the flux is q (0, t) = Q 0 . The solution to Equation (34) is therefore given by q = Q 0 . For uniaxial flow and deformation, the deformation gradient tensor, F, is diagonal, which greatly simplifies the stresses: only the xx-components of σ and K, denoted σ 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 σ xx and K xx : We now replace the fluid fluxes in the 1D evolution equation for φ g , given by Equation (7a), with the constitutive behaviour outlined above in Equations (26) and (32). The gas fraction then evolves according to where Ψ α is now the dimensionless capillary potential, which is specified below. Similarly, we combine the 1D evolution equation for φ (Equation 7b) with the relevant constitutive behaviour to give We close our model with the 1D version of Equation (27), which describes the distribution of damage, where W + is the dimensionless tensile energy, which in 1D is solely a function of φ. The model is thus reduced to two non-linear evolution equations (for φ g and for φ) and one elliptic equation (for d).
By non-dimensionalising the capillary potentials (Equation 26b) and undertaking the required differentiation (Appendix B), we find that where the interaction potentials, Π g and Π l , are and The interaction potentials result from the thermodynamic interactions between the phases, as described by ψ 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 ≤ 1. The non-dimensional effective stress (Equation 29) thus reduces to where σ 0 is the undamaged neo-Hookean stress given by and ς = K G . The derivation of this effective stress can be found in Appendix A. Similarly, the tensile energy used in Equation (38) is given by, The incompressibility condition (Equation 4) links J to the porosity field: Finally, K xx is derived in Appendix C and is given by 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.

III. TWO-PHASE PROBLEMS
Our full three-phase model captures several different physical effects, and ultimately consists of modelling fluid-fluid phase separation within the framework of large-deformation poroelasticity. To gain insight into the behaviour of our model, we begin by considering two limiting cases involving two-phase mixtures: a gas-liquid system (no solid) and a liquid-solid system (no gas). These two limiting cases are comparatively well understood, allowing us to benchmark the predictions of our model against previous results.

A. Gas-Liquid System (no solid)
In the no-solid limit, φ ≡ φ 0 ≡ 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 terms, becoming The gas chemical potential is now given by where we have defined Γ = Γ 1 + Γ 2 − Γ 3 = Γ g + Γ l as the characteristic interfacial coefficient between the fluid phases. Equations (46) and (47)  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 bub-bles, 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.

B. Liquid-Solid System (no gas)
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 φ g ≡ 0, Equation (36) for the evolution of the gas fraction degenerates and Equation (37) for the evolution of the porosity simplifies to ∂φ ∂t with and a dimensionless mobility given by M (φ) = φ. This equation is reminiscent of standard sharp-interface poroelasticity [e.g., 49], but with the addition of a potential, Ψ, 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 (φ ∼ 1) and a porous domain (0 < φ < 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 δ (t), with a fluid-filled domain for x < δ. According to the theory of large-deformation poroelasticity, the porosity in the domain x ∈ [δ (t) , 1] then evolves according to We anticipate that our phase-field model would reduce to this formulation in the sharpinterface limit. Conservation of mass requires that the influx of liquid at the left boundary, x = δ, 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 = −M ∂σ 0 ∂x = Q 0 . We also assume that the solid packing is relaxed at its left edge, so that the effective stress vanishes at x = δ. The boundary conditions for Equation (50) are then given by As discussed above, the steady-state porosity is piecewise: 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) With these boundary conditions, we solve Equations (38) and (48)   is approximated by a smooth interpolation in the phase-field model, the width of which is controlled by Γ 2 . As Γ 2 decreases, the phase-field simulations converge toward the sharpinterface solution (Figure 5, inset).

IV. THREE-PHASE RESULTS
To investigate the behaviour of our full three-phase system, we now consider the spontaneous formation of open gas cavities in a soft porous medium, starting from a nearly homogeneous initial distribution of gas within the pore space. This scenario mimics a natural system, such as a sea bed or lake bed, in which biological and/or chemical processes produce small gas bubbles throughout the otherwise undeformed and liquid-saturated pore space. As these bubbles grow, they may separate from the pore space by opening gas-rich (solid-free) cavities (Figure 3c). A similar situation could be achieved experimentally by saturating a porous packing with a fluid containing dissolved gas. Triggering gas exsolution would then lead to the formation and growth of gas bubbles. We study the evolution of this system through numerical simulations ( §IV A) and linear stability analysis ( §IV B).

A. Numerical Simulations
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 where S 0 is the initial gas saturation, η is a field of small, random fluctuations, and φ * 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 D (Figure 7). For a stiff porous medium (D 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 (D 1), gas is able to open macroscopic cavities within the packing as described above. We investigate the impact of other parameters (such as χ gs and S 0 ) in §IV B below.

B. Linear Stability Analysis
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 φ 0 = φ 0 − |φ * |, whose pore space is occupied by a homogeneous distribution of gas with saturation S 0 and volume fraction φ g0 = S 0 φ 0 . We assume that perturbations of size 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, where φ g and φ characterise the O( ) contributions to φ g and φ, respectively.
Substituting this ansatz (Equations 56) into our model (Equations 36-45) and linearising in leads to a set of equations that describe the evolution of small perturbations. To O( ), our model thus reduces to the following algebraic equations for φ g and φ, In the above, we have introduced the functions h αβ = ∂hα ∂φ β | φ 0 ,φ g0 , where h α is the homogeneous part of the chemical potential of phase α, and the parameters ν αβ , which are different combinations of interfacial coefficients. In general, h αβ and ν αβ depend on the dimensionless groups defined in §II F and the base states φ g0 and φ 0 , but they are constant with respect to k and s. Explicit expressions for h αβ and ν αβ are presented in Appendix F. Collecting like terms and finding the eigenvalues of Equations (57) leads to the dispersion relation where ζ 1 (k) = a 1 + a 2 k 2 and ζ 2 (k) = b 1 + b 2 k 2 + b 3 k 4 . Explicit forms of a i and b i are given in Appendix F.
As noted above, s is the growth rate of small perturbations to the homogeneous system described by Equations (56). The dispersion relation allows us to calculate s as a function of the wavenumber, k, of a particular perturbation. If s < 0 for all values of k, then perturbations will decay and the system is stable -gas will remain in the pore space of the solid skeleton. If s > 0 for any value of k, then perturbations with this particular wavenumber (or range of wavenumbers) will grow exponentially, and the system will be unstable -phase separation will lead to macroscopic gas cavities. The solution of Equation with small-k limit which shows that s > 0 at small k whenever either a 1 < 0 or b 1 < 0. In these cases, all perturbations with wavenumbers between zero and some cut-off value k c will be unstable and lead to phase separation. The cut-off value k c is defined by the non-zero roots of ζ 2 (k). In addition to the region of instability for k ∈ [0, k c ], our system has the unusual characteristic of forming an unstable band of wavenumbers k ∈ [k a , k b ], such that k a > 0 under certain circumstances. This characteristic occurs when a 1 > 0 and b 1 > 0, and ζ 2 (k) has two distinct non-zero roots, which requires that both b 2 < 0 and 4b 1 b 3 < b 2 2 . The condition for stability is thus that As a i , b i are complicated functions of many different parameters, we further explore the consequence of these stability conditions via a phase plane. Specifically, we plot the spinodal (or neutral stability) curves, where s = 0, which separate stable and unstable regions of the parameter space. Equation (59) also reveals the possibility of oscillatory modes of instability, for which s is complex, under certain parameter combinations. For the parameter space investigated below, these oscillatory modes occur so close to the spinodal curves that they are not observed in the simulations, and, as such, we do not explore them any further here.

Onset of Phase Separation
Three key parameters with regard to the formation of gas cavities are the deformability of the solid matrix, D, the strength of the interaction between the gas and solid phases, χ gs , and the initial gas saturation, S 0 . Recall that D measures the ratio of mixing energy to elastic energy, whereas χ gs measures the energetic cost of gas-solid interactions. Typically, χ gs is a function of physical properties such as the size and wetting characteristics of the solid particles. If the grains are smaller or if the gas phase is more strongly non-wetting to the solid, then the energetic cost of gas invading the pore space will be larger, resulting in in χ gs corresponds to an increase in capillary potential due to, for example, a smaller grain size, making it more energetically costly for gas to remain within the pore space and therefore promoting phase separation. All other parameters are the same as in Figure 6.
a larger value of χ gs . The initial gas saturation S 0 measures how much gas is present in the system. If there is insufficient gas present, then capillary effects will not be strong enough to open macroscopic cavities and phase separation will not occur.
The stability condition defined in Equation (61) allows us to identify the regions of the parameter space in which phase separation occurs and the regions in which it does not. In The deformability of the solid skeleton has a strong influence on the onset of phase separation. At low D, the solid matrix is too stiff to be deformed by capillary forces and the gas will remain within the pore space. For larger values of D, however, the solid skeleton is more easily deformed, and capillary forces may be able to overcome the elastic resistance and open macroscopic cavities, provided that S 0 is large enough (i.e., that enough gas is present). Note that, for a particular value of χ gs , there is a critical value of D below which phase separation cannot occur, regardless of S 0 .

C. Characteristic Cavity Size
In addition to the parameters that control the onset of phase separation, we also investigate the characteristic size of the gas cavities, ∆ g . The characteristic cavity size can be estimated from the linear stability analysis by finding the wavenumber of the most unstable mode at each point in the parameter space, k * , and estimating ∆ g ≈ π k * . Doing so for a range of D then allows us to predict cavity size as a function of deformability. We compare this prediction with the mean cavity size at steady state from numerical simulations ( Figure 9).
As D increases, ∆ g decreases. However, the total amount of gas within the cavities increases with D, so this decrease in ∆ g is more than offset by an increase in the total number of cavities. This behaviour can also be seen qualitatively in Figure 7: there are more, smaller cavities for D = 5 than for D = 2 and D = 2.5. Intuitively, we would expect larger cavities to form in a softer material since the solid skeleton can undergo larger deformations. Indeed, the total 'cavity volume' is larger for softer materials, but these cavities are smaller and more numerous. We attribute this observation to the fact that, in a softer material, a smaller amount of gas is needed locally to open a cavity. For a stiff porous material, a large amount of gas will need to accumulate in order to force open a cavity, and so when the cavity is formed, it will be correspondingly larger.
Our numerical simulations show remarkably good qualitative and quantitative agreement with the linear stability analysis for the value of ∆ g , which is surprising given the strongly nonlinear nature of this process. Figure 9 also shows that, as expected, there exists a certain critical deformability below which phase separation does not occur. The linear stability analysis and numerical simulations show this transition occurring at around the same value of D.

V. CONCLUSIONS
The key feature of two-phase flow in a deformable porous medium relative to a rigid porous medium is the ability of the non-wetting phase to form open cavities within the solid skeleton. Motivated by the formation of non-wetting gas cavities in sediments, we have derived a thermodynamically consistent phase-field model that describes the formation of cavities within a soft porous material as a result of gas-liquid-solid phase separation, and have explored our model in a one-dimensional setting. In the two limiting cases of no-solid and no-gas, our model reproduces the expected behaviour for these well-studied systems.
When no solid is present, we reproduce the fluid-fluid phase separation and coarsening behaviour characteristic of the Cahn-Hilliard equation; when no gas is present, our model matches well with analytic solutions to traditional sharp-interface Biot poroelasticity.
In the full three-phase system, phase separation is inhibited by elastic resistance from the solid, which imposes limits on the conditions under which cavities form: if the solid skeleton is too stiff, gas will remain within the pore space. We have shown how the onset of phase separation depends on different parameters, including the deformability of the solid matrix and the wetting characteristics of the two fluids, via both linear stability analysis and numerical simulation. We have also shown that, for a softer porous material, we expect smaller, but more numerous, cavities than for a stiffer porous material, owing to a smaller amount of gas needing to accumulate within the pore space in order to open a cavity in the former case.
Our work has been motivated by the venting of gas from non-cohesive granular sediments, in which the compressibility of the gas phase can be neglected when the capillary entry pressure of the granular skeleton is sufficiently small compared to the ambient liquid pressure.
Our model is also relevant for different fluid pairs in other contexts, such as phase-separating liquid droplets in polymer networks [16][17][18]. In polymer systems, the porous skeleton does not have a granular microstructure, and the energy required to fracture the polymer network to form cavities affects the final size of droplets [50][51][52]. Our model can capture fluid-fluid phase separation within such porous materials through an appropriate choice of the Griffith fracture energy.
We have solved our model assuming a neo-Hookean elastic response for the solid skeleton, but our theory allows for other elastic laws. We have assumed that the gas phase is incompressible, and that the gas and liquid phases are immiscible. We anticipate that relaxing these assumptions will lead to further interesting physical phenomena, the investigation of which will be the subject of future work. Our model makes use of a phenomenological free energy. Connecting the parameters used in our free energy to physical quantities remains an avenue for future study; considering the pore-scale thermodynamics of a gas-liquid-solid mixture could provide this link. Ongoing experimental work will allow us to compare our theory with experimental observations, and will focus on the dynamic transition between the phase-separated and homogeneous regimes via the application of external confining stress. As noted in §II D, our system displays distinct mechanical behaviours depending on whether the material is in tension or compression. As such, we must identify when each regime occurs. Following the approach of Tang et al. [41], we first rewrite the neo-Hookean strain-energy density in terms of the principal stretches, λ i , where we have used that J = λ 1 λ 2 λ 3 . The first part of this expression corresponds to stretching deformations and the second part to volumetric deformations. We then decompose this strain-energy function into tensile and compressive parts, given by W + = W el λ + i , J + and W − = W el λ − i , J − , respectively. We define λ ± i and J ± as piece-wise functions, with and As per Equations (22), the tensile part of the elastic energy, W + is degraded by an increase in the damage parameter.
The effective stress, σ , is found by taking the derivative of the elastic free energy (Equation 29), Using the chain rule, we have that (following [41]), where and This constitutive law reduces to a standard neo-Hookean behaviour in the case of an undamaged system (d = 0).
The solid mechanics of our system greatly simplify in the uniaxial case, as described in §II F. In 1D, u s = u s (x, t)ê x , and hence using Equation (2), we see that the deformation gradient tensor is diagonal, with Noting that J = det (F), we thus write that In this case, the principal stretches λ i are the eigenvalues of the deformation gradient, and so only one of these stretches λ = λ 1 = J is non-constant in this system. As such, the neo-Hookean strain-energy can be written solely in terms of J, and we can identify tension as being when J > 1. The above constitutive behaviour (Equation A4) thus simplifies to with the tensile and compressive components of the free energy defined as, and respectively, where W 0 el = 1 2 G (J 2 − 1 − 2 log J) + 1 2 K (log J) 2 . Undertaking the differentiation of the free energy with respect to J then gives, where σ 0 is the undamaged neo-Hookean stress, given by, To complete our description of the mechanics, we link J to the porosity, φ, using Equation (4).
where we have further defined the interaction potentials, Π α , as the component of the capillary potential resulting from the differentiation of ψ mix . These are given by: For incompressible, immiscible fluids, ρ α µ 0 α is constant. In our model, the capillary potential only appears as the argument of a gradient function, and so this constant bulk term has no effect on the overall state of the system. The spatial derivatives in the capillary potentials resist the formation of sharp gradients in volume fraction, and serve to regularise the system by preventing the formation of shocks in the gas fraction and porosity.

Appendix C: Korteweg Stress
In order to derive an explicit expression for the Korteweg stress, we first consider the free energy of a general interface between two unspecified phases, α and β, using the form given in Equation (24), ψ αβ = γ 2 (∇φ α · ∇φ β ) J.

Appendix E: Analytical Solution to Sharp-Interface Poroelasticity
To find a steady-state solution to the sharp-interface model described in §III B, we first note that the porosity is described by a piecewise function, we first define h α as the homogeneous part of the chemical potential of fluid α, given by and To derive h αβ , we then differentiate these expressions with respect to φ g and φ, and evaluate them at φ = φ 0 and φ g = φ g0 . This gives To find explicit forms for ν αβ , we collect the coefficients of fourth order derivatives to give The dispersion relation derived in §IV B (Equation 58) is written compactly in terms of the functions ζ 1 (k) = a 1 + a 2 k 2 and ζ 2 (k) = b 1 + b 2 k 2 + b 3 k 4 , where We identify M i , as linearised mobility coefficients.