A Monte Carlo Model of Gas-Liquid-Hydrate Three-phase Coexistence Constrained by Pore Geometry in Marine Sediments

Gas hydrates form at relatively high pressures in near-surface, organic-rich marine sediments, with the base of the hydrate stability field and the onset of partial gas saturation determined by temperature increases with depth. Because of pore-scale curvature and wetting effects, the transition between gas hydrate and free gas occurrence need not take place at a distinct depth or temperature boundary, but instead can be characterized by a zone of finite thickness in which methane gas bubbles and hydrate crystals coexist with the same aqueous solution. Previous treatments have idealized pores as spheres or cylinders, but real pores between sediment grains have irregular, largely convex walls that enable the highly curved surfaces of gas bubbles and/or hydrate crystals within a given pore to change with varying conditions. In partially hydrate-saturated sediments, for example, the gas–liquid surface energy perturbs the onset of gas–liquid equilibrium by an amount proportional to bubble-surface curvature, causing a commensurate change to the equilibrium methane solubility in the liquid phase. This solubility is also constrained by the curvature of coexisting hydrate crystals and hence the volume occupied by the hydrate phase. As a result, the thickness of the three-phase zone depends not only on the pore space geometry, but also on the saturation levels of the hydrate and gaseous phases. We evaluate local geometrical constraints in a synthetic 3D packing of spherical particles resembling real granular sediments, relate the changes in the relative proportions of the phases to the three-phase equilibrium conditions, and demonstrate how the boundaries of the three-phase zone at the base of the hydrate stability field are displaced as a function of pore size, while varying with saturation level. The predicted thickness of the three-phase zone varies from tens to hundreds of meters, is inversely dependent on host sediment grain size, and increases dramatically when pores near complete saturation with hydrate and gas, requiring that interfacial curvatures become large.


INTRODUCTION
Natural gas hydrates are ice-like compounds that commonly form in permafrost and marine sediments from mixtures of methane and water (Sloan and Koh, 2007). As a promising source for future energy, methane hydrate has attracted much attention from the oil and gas industry, with further motivation for their study coming from the need to quantify methane migration in sediments (e.g., Nole et al., 2016), assess submarine landslide risk (e.g., Sultan et al., 2004;Handwerger et al., 2017), and understand the material cycle in benthic ecology (e.g., Suess et al., 1999). Seismic data and drilling logs from natural hydrate reservoirs have identified anomalies of high saturation level (i.e., hydrate pore volume fraction) within layers of comparatively coarse sediments, suggesting heterogeneous hydrate accumulation rates that depend not only on temperature and pressure but also on sediment properties (e.g., Borowski, 2004;Malinverno, 2010;Wang et al., 2011;Bahk et al., 2013). Experimental studies also demonstrate that pore sizes play an important role in controlling the spatial and temporal distribution of hydrate deposits (e.g., Chong et al., 2015;Liu et al., 2015).
The formation of gas hydrate in permafrost and marine sediments is often approximated using the constraint of local bulk equilibrium between a combination of up to three methanebearing phases: free methane gas (G), methane hydrate (H), and dissolved methane in aqueous solutions (L), in which the methane solubility is a unique function of temperature, pressure and salinity (e.g., Sloan and Koh, 2007). A more precise understanding of these systems must account for perturbations to this bulk phase behavior imposed by the surface properties and geometry of sediment particles, with the gas and hydrate acting most commonly as non-wetting phases, whereas the aqueous solution wets particle surfaces. As hydrate forms or dissociates, hydrate crystals or gas bubbles approach the pore walls, and the phase behavior is affected by the surface energy of the curved L-H or L-G interface, with high curvature causing elevated local dissolved methane concentrations. By constraining allowable interface curvatures, heterogeneously distributed sediment pores introduce deviations in the equilibrium methane concentration at in situ temperature and pressure conditions Henry et al., 1999;Daigle and Dugan, 2011;Rempel, 2011;Dai et al., 2012;Cook and Malinverno, 2013;VanderBeek and Rempel, 2018), thereby affecting both the growth of hydrate deposits and their decomposition. Existing works that approximate the role of pore geometry mostly focus on the average pore size, often simplifying the pores as circular cylinders (e.g., Millington and Quirk, 1961;Wilder et al., 2001;Denoyel and Pellenq, 2002) or spheres connected by cylindrical throats (e.g., Jang and Santamarina, 2011;Liu and Flemings, 2011). These simple pore models provide useful insight into how hydrate forms and dissociates in sediments, but they fail to capture variations in curvature as phase boundaries evolve. Rempel (2011) avoided this limitation by considering triangular pores, and a subsequent two-dimensional treatment (e.g., Rempel, 2012) examined the crevice spaces between random close-packed spheres. By treating granular porous media as packed three-dimensional spherical grains, Chen et al., (2020) used Monte Carlo sampling to effectively approximate the constraints of pore geometry on phase boundary curvatures in a two-component system within randomly packed, poly-dispersed sediments. In this work, focused on three-phase coexistence, we first outline the basic phase behavior expected within 2D triangular pores, and then extend the treatment using an averaging method to approximate the behavior in pores between spherical grains, before examining the fully 3D problem with a Monte Carlo method.

EQUILIBRIUM METHANE CONCENTRATION GRADIENT IN SEDIMENTS
In marine sediments that are sufficiently coarse-grained for porescale curvature effects to be negligible, bulk three-phase equilibrium at the base of the hydrate stability zone (BHSZ) occurs at a distinct depth that is uniquely determined by the pressure, temperature and salinity. Above the bulk BHSZ, the equilibrium methane solubility of the binary L-H system increases with depth, whereas below the BHSZ, the equilibrium is between liquid and free gas, and the methane solubility decreases with depth, driven by increases in the ambient temperature. In typical circumstances with heterogeneously distributed micron-scale pores, however, the hydrate phase, gas phase and aqueous methane solution may coexist in a zone of finite thickness where the upper and lower boundaries are shifted according to the solubility perturbations associated with confining the hydrate and gas phases in tight, and variable effective pore sizes.
The shift in methane solubility from bulk conditions in L-H and L-G two-phase equilibrium can be approximated as follows. For the L-G equilibrium, the thermodynamic relations for the methane solubility in molar fraction require where Δ gl sol H m is the molar heat of solution of methane gas (negative for this exothermic reaction), V m is molar volume of methane gas, and V m is the partial molar volume of methane in water, which is negligible compared with V m . Partial pressure from water vapor is also negligible because in the temperature range of interest, the saturation vapor pressure is less than 1 kPa, which is much smaller than the hydrostatic pressure. Similarly for the L-H equilibrium, Adopting a coordinate axis with the z-direction pointing vertically downwards, the bulk BHSZ is at depth z 3 below the seafloor, corresponding to a three-phase equilibrium condition where T 0 and P 0 are the temperature and pressure at the seafloor, and G T and G P are the temperature gradient and pressure gradient, respectively, in the sediment. The respective solubilities vary with depth near z 3 according to where the pressure dependence of L-H solubility is in fact negligible because the volume change is relatively small without the presence of a free gas phase. For illustration, we consider perturbations around the three-phase equilibrium T 3 ≈ 295 K and P 3 ≈ 30 MPa, and use nominal values for G T and G P at Blake Ridge (Table 1) so that g gl −0.309 km − 1 , g hl 2.09 km − 1 .
Because g hl > g sl (i.e., |d ln x hl z| > d ln x gl z ), the gradient of the gas solubility is much gentler than that of the hydrate solubility, as depicted in Figure 1. The bulk solubilities at z z 3 + Δz are approximately Curved surfaces of gas bubbles and methane hydrates within the confined pore space elevate the chemical potential of the nonwetting gas and hydrate phases. At a depth z where three phases coexist, setting the radius of the methane bubble to r g , and the radius of the hydrate crystal to r h , the shifted solubilities are x ′ gl r g x gl 1 + c gl P 3 + G P Δz 2 r g , Equilibrium between the phases requires which is expanded to Equation 12 describes the chemical equilibrium when three phases coexist. With r h and r g constrained from the pore distribution, the offset Δz gives the thickness of the threephase zone. In simple porous sediment model with a single pore size so that r h r g , the three-phase zone shrinks to one unique depth. With heterogeneously distributed effective pore sizes, however, we expect the BHSZ to be characterized by a zone of three-phase coexistence bounded by depths corresponding to equilibrium conditions for which the hydrate crystals and gas bubbles each have different interfacial curvatures ( Figure 1) that are nevertheless related by the constraint that each of these nonwetting phases must also be in equilibrium with a wetting aqueous solution containing the same concentration of dissolved methane. Combined with the geometric constraints derived below, Eq. 12 enables us to determine the thickness of the three-phase zone.

GEOMETRIC CONSTRAINTS WITHIN PORES
Models consisting of regular pores with concave interior walls, such as spheres or cylinders, permit very little variation to phase boundary curvature, as non-wetting phases fill pore centers and the wetting phase occupies thin films that coat pore walls. In natural irregular pores with predominantly convex interior walls, as the non-wetting phase grows, the phase boundary intrudes further into crevices between solid grains where the wetting phase persists in ever-shrinking convexly bounded pockets. One simplified pore model with features resembling such diminishing crevices is a 2D triangular pore; a more realistic 3D model can be constructed using a conglomerate of packed  (Anderson, 2004;Gupta et al., 2008).

Model parameters Value
Methane molar dissolution heat (Duan and Mao, 2006) Δ gl sol H m [kJ mol −1 ] −12.59 Hydrate molar dissolution heat (Lu et al., 2008) Δ hl sol H m [kJ mol −1 ] 41.96 Molar volume of water (Wagner and Pruss, 1993) V w [cm 3 mol −1 ] 17.93 Partial molar volume of methane in water (Duan and Mao, 2006) V m [cm 3 mol −1 ] 38.87 Molar volume of hydrate (Sun and Duan, 2005) V h [cm 3 mol −1 ] 135.4 Hydration number n ∼ 6 Geothermal gradient (Ruppel, 1997) G T [K m −1 ] 3 .69 × 10 − 2 Hydrostatic pressure gradient (Ruppel, 1997) G (Hardy, 1977) c particles, idealized here as spherical grains. With changing temperatures and/or pressures, for example with increased depth below the seafloor, the hydrate phase (H) is expected to dissociate and a new non-wetting phase (G) will emerge from the wetting phase (L) so that a two-phase equilibrium (L-H) configuration gives way to a new three-phase equilibrium (G-L-H). The zone of three-phase coexistence may have a finite thickness with varying saturation levels for the non-wetting phases, before reverting to a different two-phase equilibrium (L-G) at still greater depths. Importantly, at the onset of threephase coexistence, surface energy considerations imply that the emergent phase (in this scenario, G at the top, H at the bottom) is bounded by the largest surface of constant curvature that can fit within the pore space (i.e., a sphere). This simplifies the geometry of the emergent phase considerably and facilitates determination of the three-phase zone thickness while avoiding the need to consider the regions of variable curvature adjacent to the extended wetting films that coat both non-wetting phases elsewhere. A second useful constraint is that the curvature of the residual phase (in this scenario, H at the top, G at the bottom) must remain continuous across the three-phase boundary.
With the two constraints, we can describe the evolution of saturation levels of the non-wetting phases when crossing the boundary from regions of two-phase equilibrium into the zone of three-phase equilibrium. For example, in the L-H region immediately above the three-phase zone, hydrate is the only non-wetting phase in pores, separated from pore walls by films of liquid phase. The hydrate crystals have a radius r h controlled by the surface energy. At low hydrate saturations, the crystals may take a spherical form, whereas at high saturation levels (with abundant methane), small spheres may coalesce, and occupy the largest pore with bumps growing into nearby crevices with the radius r h . Across the three-phase boundary, a portion of hydrate dissociates, and spherical methane gas bubbles emerge with radius r g > r h tangent to the walls of the largest pores and hydrate crystals, whereas remaining hydrate resides in smaller pores and extends into crevices with the same interfacial radius r h as the crystals in the L-H region right above. At the base of the three-phase zone, where hydrate is the emergent phase and free gas the residual phase, a parallel set of constraints applies with continuous gas radius r g and spherical hydrate crystals characterized by r h > r g . For the 2D triangular pore model, an analytical description of these geometrical constraints is available; for the 3D model, we developed an averaging method to FIGURE 1 | The three-phase coexisting zone near the bulk BHSZ. Above the BHSZ, no methane gas is present, and the methane solubility is determined by L-H equilibrium, increasing with depth (green curves). Below the BHSZ, hydrate dissociates so that dissolved methane is instead constrained by equilibrium with free methane gas, and the solubility decreases with depth due to increasing temperature (red curves). Bulk solubility curves correspond to the scenario where pore-scale effects can be neglected. In smaller pores, however, two-phase solubility curves shift toward higher values. The hydrate and methane gas phases can first coexist with the same aqueous solution when the emergent free gas phase at the upper boundary of the zone of three phase coexistence has the smallest possible curvature (i.e., largest radius), while the curvature that characterizes crystals of the residual hydrate phase must remain continuous with the value set by the hydrate saturation level in the two-phase L-H zone above; a parallel set of restrictions pertains at the lower boundary with the roles of gas and hydrate reversed. The dark L-H-G line labels the methane solubility such that even the smallest pores are filled with one non-wetting phase.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 600733 represent a mono-dispersed scenario. We describe the geometrical constraints for these two idealized cases next, before outlining our Monte Carlo approach to addressing a more realistic synthetic sediment consisting of randomly packed spherical particles in Section 3.2.

Simplified 2D Triangular Pores
In a 2D equilateral triangular pore with sides of length W, the radius of the residual non-wetting phase I near the vertices is R 1 (Figure 2A). The total pore area is In the case where its boundaries are idealized as spherical, the total area of non-wetting phase I in three vertices is where the inequality means not all vertices are necessarily hosting phase I. When a new non-wetting phase II emerges, it must have lower surface energy (i.e., larger radius) so it locates near the center of the pore, with R 2 ≈ 3R 1 , and where h W/2 3 √ , and d is the film thickness far from vertices along pore walls. The saturations of each phase are In the limit that R 2 ≫ d, the saturation of the emergent phase is Here, the pore geometry requires R 2 ≤ h W/(2 3 √ ) so phase II has minimum saturation S 2 ≈ 0.6, and at the onset of three-phase coexistence S 1 may have a range of values below 0.2. However, it remains possible for the two non-wetting phases to occupy separate nearby pores as long as R 2 ≈ 3R 1 , and S 2 remains a valid description of the saturation level of the emergent phase in equilibrium with L alone at the onset of three-phase coexistence. These geometric constraints apply only at the boundaries of the three-phase zone, and not further within the zone itself. Instead, the gas and hydrate interfacial curvatures within the interior of the three-phase zone depend on the amount of methane present, since it must be partitioned between L, H, and G. It is possible that inside the zone, pores are under-filled, i.e., neither residual phase I and emergent phase II are above a saturation level of 0.6, which we will discuss later.

Mono-Dispersed 3D Pores
A similar approach combined with an averaging method can be used to obtain saturation estimates in 3D pores. In a monodispersed 3D sediment with particle radii R, the entire volume of each pore in a virtual triclinic cell bounded by eight grains with internal angles α, β, and c is Using the hyper-volume formula (Mackay, 1974), the radius of the largest inscribed sphere is 3 − 2 cos α − 2 cos β − 2 cos c 1 − cos 2 α − cos 2 β − cos 2 c + 2 cos α cos β cos c + 2 cos α cos β + 2 cos β cos c + 2 cos c cos α 1 − cos 2 α − cos 2 β − cos 2 c + 2 cos α cos β cos c − cos 2 α + cos 2 β + cos 2 c 1 − cos 2 α − cos 2 β − cos 2 c + 2 cos α cos β cos c Phase II attains equilibrium as the new non-wetting phase with radius R 2 ≤ r(α, β, c) given by Eq. 19. Similar to the 2D case, the emergent phase II may appear as spheres with radius R 2 ≤ r(α, β, c) adjacent to phase I, or as a contorted body intruding into all possible interstitial sites, with a surface characterized by small bumps of R 2 tangent to the bounding sediment particles to form crevices. The residual non-wetting phase I may stay inside one or more crevices, as shown in Figure 2B. Neglecting the small volumes contributed by liquid films, as before, we have FIGURE 2 | Schematic of three-phase coexistence in (A) one 2D equilateral triangular pore and (B) residual liquid reservoir inside one crevice in 3D spherical grains. As the non-wetting phase II emerges, the area or volume occupied by non-wetting phase I shrinks.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 600733 where δ arcsin[R/(R + R 2 )]. The wetting phase volume is wellapproximated as filling a trio of minor crevices (two of which are bounded on one side by phase II while the other sides approach particle surfaces), so that the wetting volume is Here, the total volume V 1 occupied by phase I is bounded between the volume of balls of radius R 1 and that of three tori while the volume occupied by phase II is Finally, the saturation levels can be written as For the collective values of R 1 , R 2 , and S 2 over numerous pores, these values are averaged over angles α, β, and c (see Appendix 1 for details). At the top and bottom of the three-phase coexisting zone, we recognize the emergent phase II as the gas and hydrate phases, respectively. The analyses for the 2D and 3D scenarios suggest that the saturation of the residual phase S 1 can vary within a range, while the saturation of the emergent phase S 2 is better constrained.

Monte Carlo Simulation of 3D Pores
In natural, randomly packed sediments, clearly the virtual cell of Section 3.1.2 may be heavily distorted, and the distributions of angles are affected by grain radii, so the averaging method may not work properly. We develop a Monte Carlo scheme to simulate the growth of the emergent phase as constrained by the pore geometry, as well as the requirement imposed by continuity of solubility. We test the method here using the mono-dispersed random close pack of Finney (1970), and sample the cross-section of the pack with N 2000 random test points. For each test point located in the pore space, we find the largest inscribed sphere containing the test point, which is recognized as emergent phase II with a radius R 2 , and in the crevices formed between the sphere and two tangent particles, we calculate a tangent coplanar sphere as residual phase I. There may be more than one possible crevice in each pore because the phase II sphere may touch as many as four particles, and we choose the largest residual phase I sphere, with radius R 1 . We record all pairs of phase I and II sphere radii (R 1 , R 2 ), and sort them according to the values of R 2 . This sorting procedure enables us to approximate the saturation level of phase II with radius R 2 as the proportion of sampled points that are encompassed by phase II (see Chen et al., 2020, for a more detailed discussion). After scaling with particle radius, we solve for the corresponding Δz using Eq. (12). The Monte Carlo sampling procedure results in different estimates of Δz for estimates of the emergent-phase saturation S 2 , and we recognize the envelope of extremal values as approximating the depth range of three-phase coexistence.

MODEL RESULTS
We seek the upper and lower depth limits that define the zone where three-phase equilibrium may occur. At the top of this zone, free gas is the emergent phase II and hydrate is the residual phase I, whereas at the bottom these roles are reversed, with the gas constituting the residual phase I and hydrate the emergent phase II. By applying the geometric constraints derived in the 2D and 3D scenarios just described to Eq. 12, we can determine the dependence of Δz on S 2 . For uniform 2D triangular pores ( Figure 3A), the pore geometry requires R 2 ≤ W/(2 3 √ ), so that the minimum S 2 ≈ 0.6. For 3D pores in mono-dispersed grains ( Figure 3B), the thickness of the three-phase zone is the average over all pores, and the minimum S 2 is around 0.75. Because the z-direction points downwards, the figures are plotted with flipped Δz so that shallower locations (negative Δz) are above deeper locations (positive Δz).
The two scenarios behave similarly. With larger pores, the zone of three-phase coexistence is thin, but as pore size decreases (represented by the different lines in Figure 3, with sizes noted in the legends), the upper and lower boundaries of the three-phase zone deviate further from zero, corresponding to a thicker zone of three-phase coexistence. In the smaller pores of finer sediments, the hydrate phase begins to dissociate and the gas phase emerges at a depth much shallower than the bulk BHSZ, but the hydrate phase may also persist to a depth far below the bulk BHSZ. The thickness of the zone of three-phase coexistence is constrained as well by the requirement that solubility remain continuous across the boundaries with adjacent two-phase zones, leading to the dependence on S 2the saturation of emergent phase. Figure 4 compares the 3D average result from Section 3.1.2 with the Monte Carlo simulation result from Section 3.2. The upper and lower bounds of the Monte Carlo results match well with the average curves at the beginning, but deviate further at high saturation levels. In finer sediments, boundaries marked by the Monte Carlo results begin to deviate further from the averaging result. We attribute this discrepancy to errors in the averaging procedure produced by distortions to the virtual cell.
We emphasize that it is not necessarily the case that any particular pore in the zone of three-phase coexistence can hold all three phases, and in fact such a configuration is only possible at very high methane input. Nevertheless, since the volume occupied by phase I determines its interfacial curvature and hence the methane solubility in the adjacent two-phase zone, together with the pore size constraint on the geometry available for phase II, this implies that the three-phase thickness (i.e., bounded by the first appearance of a secondary nonwetting phase) must depend on both the saturation of the primary (residual) non-wetting phase and the pore size distribution.

Growth of the New Phase
The geometric constraints applied in the 2D and 3D scenarios treat the emergent phase as volumetrically dominant, limited in Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 600733 extent by the pore walls and residual phase. It must be noted that the radius of the emergent phase may be restricted by the presence of the residual phase; when the new phase nucleates, it is assumed to form near an interface with the residual phase, essentially replacing much of the pre-existing phase I to reach the minimum free-energy configuration while maintaining the same phase I curvature as that which pertains outside the three-phase coexisting zone.
Alternatively, the new phase could grow in the largest pores, either without being adjacent to the residual phase, or by completely replacing the residual phase that would have occupied those pores under the slightly perturbed conditions in the adjacent two-phase zone. In this situation, to satisfy the continuity of phase I curvature with that outside the three-phase zone, phase I can persist either in the form of small residual inclusions within pore crevices, or as a body filling almost all of FIGURE 3 | The shifted three-phase boundary Δz b and Δz t at the first appearance of the emergent phase II as a function of the saturation level S 2 , shown for the different values of (A) pore size W and (B) grain size R noted in the legend. Note that since z is pointing downwards, we have the z-axis flipped so that shallower location stays above deeper locations.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 600733 7 smaller pores with bumps of small radii R 1 . In the first case, the small radius characterizing the residual phase in the two-phase region suggests high methane solubility, which is unstable because the chemical potential can be further minimized by increasing R 1 and reducing methane solubility. The second case leads to unrealistically high saturation levels for the hydrate crystals or methane bubbles.

Under-Filled Region Within the Zone
In our models, both non-wetting phases are mobile in the pores as long as not limited by the pore walls and the other non-wetting phase. At the boundaries, the emergent phase spans the pore center while the residual phase stays in the crevices. If the residual phase occupies a significant fraction of the pore, emergent phase may not be able to touch the pore walls, resulting in a lower saturation (< 0.6 for the 2D pores and < 0.75 for 3D pores). This under-filling scenario can be intuitively investigated for the 2D model, where the residual phase has a radius R 1 > W/(6 3 √ ), or S 1 > 0.2 ( Figure 5). The value of R 2 is smaller accordingly, and so is the value of S 2 . Ignore d and let R 1 αW, and we can find the corresponding R 2 where 2 3 √ ≤ α −1 ≤ 6 3 √ . In this configuration, R 1 and R 2 are symmetric, and it is easy to calculate that the offsets is smaller than the configuration in Figure 2A. We postulate that for the 3D pores, the under-filled configuration also gives smaller offsets. The three-phase zone is bounded by the maximum offset possible with given pore structure, pressure, and temperature, and under-filling cases are located in between the maxima. In pores located near the middle of the three-phase zone, both non-wetting phases may be under-filling, allowing the transition from L-H above the three-phase zone to L-G beneath the three-phase zone. Therefore, when seeking the boundaries of the three-phase zone, we need only consider configurations in which the emergent phase fills pore centers.

Shifted BHSZ and BSR
The bottom simulating reflector (BSR) is commonly interpreted as marking the BHSZ, which is the boundary separating the hydrate phase above from the free gas phase below. However, due to perturbations in salinity and the pore scale effects described here, hydrates can still be present at equilibrium below the bulk FIGURE 4 | Monte Carlo simulation of shifted three-phase boundary Δz b and Δz t for residual phase I before the appearance of emergent phase II as a function of grain size R, and comparison with 3D averaging results (black and cyan solid lines). The upper and lower bounds begin to deviate when the saturation is high and R is significantly reduced.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 600733 BHSZ while free gas bubbles can persist above the bulk BHSZ. Our calculations show that the resulting zone of three-phase coexistence can vary in thickness from only a few meters to many tens of meters. This may cause the temperature and pressure at the BSR to deviate from three-phase equilibrium conditions, and the observed depth of BSR may differ from the bulk BHSZ. For example, the Ocean Drilling Project Leg 164 at Blake Ridge found that the temperatures at the BSR are 0.5°-2.9°C lower than the theoretical equilibrium temperature; this corresponds to an upwards shift of 30-100 m above the bulk BHSZ, assuming a geothermal gradient of ∼ 30°C/km (Ruppel, 1997). Liu and Flemings (2011) showed that it is also possible to have the top of the three-phase zone below the bulk equilibrium depth, resulting in a deeper BHSZ. Such discrepancies have been attributed to shifted hydrate equilibria in porous media. In Figure 4 the simulated top positions in some cases are lower than the bulk BHSZ. Our model predicts that a zone of three-phase coexistence can be expected to span the bulk BHSZ, with its upper boundary shifted toward colder temperatures by an amount controlled both by the hydrate saturation in the L-H region above, and the largest pore sizes available to emergent gas bubbles, giving a quantitative explanation for the observed discrepancy.

Implications for Poly-Dispersed and Non-Spherical Granular Sediments
Our model deals with simplified pores in mono-dispersed spherical grains, and our averaging method is strictly valid only for pores bounded by grains in direct contact. Theoretically, crevices can occur between separated grains, but as the distance between the two grains increases, it is much more difficult for liquid connecting the grains to form a concave meniscus with positive mean curvature. Hence, at low liquid saturations, most liquid stays in crevices between contacting grains. In real sediments, grains are poly-dispersed and irregular. If the grains are silt-sized and assumed spherical and contacting, we can model a random packing using the drop-androll method (Chen et al., 2020) with the particle sizes following a specified distribution, and apply the Monte Carlo method similarly. When the particle sizes follow a log-normal distribution ln N (μ, σ 2 ), our simulation results suggest that the shifted three-phase zone will remain mostly the same as in the mono-dispersed situation, except that the relevant grain size R should be comparable to the median radius R m exp(μ) (defined by particle count rather than weight).
Constructing realistic synthetic packings that incorporate highly non-spherical grains, as expected of sediments with significant clay contents, is a more challenging numerical problem. Chen et al., (2020) pursued a simplified strategy in which two-phase saturation predictions were performed on a mono-dispersed packing with particle radii chosen so that the specific surface area matched the measured value for a silt loam (73 m 2 /g: 33% sand, 49% silt, 18% clay by weight; Or and Tuller, 1999). Comparisons with partial saturation measurements showed excellent agreement when the non-wetting phase occupied more than 90% of the pore space, and the agreement remained acceptable down to about 60% non-wetting phase saturation. Further exploration of these results suggests that most of the residual wetting phase under these conditions is found in the increased numbers of small crevice-like regions in the vicinity of particle contacts. In sediment with significant amounts of non-spherical grains such as clay minerals, the crevice-like regions may be smaller, which limits the size of residual phase, and possibly will cause a thicker three-phase zone.

Sensitivity of Emergent Phase Stability to Residual Phase Saturation
Because residual phase saturation has a strong influence on the thickness of the zone of three-phase coexistence, it is possible to estimate the saturation of the hydrate or gas phase from BSR observations if the median particle size is known or can be estimated. However, since both the saturation level of the residual phase and the particle size of the host sediment can vary over short distances, we may expect that in some patches there may be interleaving of two-phase and three-phase zones. This further complicates the interpretation of BSR observations, and may be responsible for discontinuities in BSR location.

CONCLUSION
By approximating porous sediments as consisting of pores with diminishing crevices, we have demonstrated that near the base of the gas hydrate stability field, the upper (cold) boundary of a three-phase region is set by the gas-liquid surface energy of the first spherical bubbles that can form in the partially hydratesaturated sediments, while the lower (warm) boundary is controlled by the surface energy of the first hydrate crystals that can form in the partially gas-saturated sediments. Of more fundamental importance, our analysis shows that the thickness of the three-phase zone depends not only on the grain-size distribution, increasing dramatically from tens of meters in porous sediments with a median grain size of 1 µm to hundreds of meters when the median grain size is 0.1 µm, but that in a given sediment the thickness is also sensitive to the saturation levels of hydrate and gas at the boundaries with the two-phase zones above and below.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JC and AR conceived the original idea of how the residual phase must affect the thickness three-phase coexisting zone, and how it might be calculated. SM helped JC with the numerical implementation of these ideas in the Monte Carlo calculations. All authors contributed toward the writing, interpretations, and editing.