1
Supplementary Information for
2
Crustal fingering facilitates free-gas methane migration through the hydrate stability zone
3
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
4
This PDF file includes:
5
Supplementary text
6
Figs. S1 to S7
7
Table S1
8
Captions for Movies S1 to S4
9
References for SI reference citations
10
Other supplementary materials for this manuscript include the following:
11
Movies S1 to S4
12
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
1 of
13
www.pnas.org/cgi/doi/10.1073/pnas.2011064117
Supporting Information Text
13
Laboratory experiments
14
Observation of solid crust formation.
During the depressurization phase, direct visualization of newly formed hydrate crust
15
becomes challenging, but is apparent in certain regions where the gas-liquid interface appears rough. In addition, the formation
16
of solid hydrate crust can be inferred indirectly. If solid hydrate formation does not take place, then the expansion of the
17
single gas bubble is expected to be stable and uniform, analogous to that of balloon inflation. However, during the controlled
18
depressurization, we do not observe stable expansion of gas into water (Fig. 2c). In fact, the gas volume remains unchanged for
19
a few minutes while pressure decreases, suggesting that the gas phase is initially trapped by the pre-existing hydrate crust.
20
The pressure difference between gas and liquid eventually builds up enough to rupture the hydrate crust (Fig. 2c,
t
=
2.5min),
21
allowing gas to displace into the ambient liquid. The consequent displacement pattern is finger-like, suggesting that a solid
22
layer of hydrate along the gas-liquid interface is spontaneously forming and modulating the direction of gas flow.
23
Validation of hydrate formation.
Here, we provide evidence that the observed solid that forms along gas fingers indeed represents
24
gas hydrates:
25
1.
Control experiment with water and air bubble
under the same condition (25
◦
C,
≤
7
.
5
MPa) does not show the
26
formation of a solid around the gas bubble;
27
2.
Direct visual observation of hydrate dissociation.
When the pressure is within the hydrate stability zone, the
28
roughness and jaggedness of the interface means that there is a solid forming (
3
,
4
) (see also SI video 1 and Figure 2a).
29
Once the pressure exits the stability zone, we observe that solid hydrate dissociation is two-staged. During the first stage,
30
the solid crusts along fingers gradually disappear, showing hydrate melting into the surrounding liquid (Figure S2a). The
31
dissociated hydrate leaves behind water that is supersaturated with xenon. During the second stage, the supersaturated
32
solution undergoes spinodal decomposition, leading to exsolution of gas bubbles (Figure S2b).
33
Effect of subcooling on hydrate growth along the fingers.
Subcooling (T
eq
-T) plays an important role on the rate of hydrate
34
growth and should control the pattern formation of crustal fingers. In particular, the rate of growth along the interface (
R
s
)
35
is positively linked to the subcooling with a power law behavior (
3
). In our experiments, T
eq
-T (the distance between the
36
red arrow and the blue phase boundary in Figure 1b) is not fixed as we depressurize. Rather, it decreases with time (and
37
pressure). Therefore, we should expect that crust that forms later in the experiment to grow slower and thus appears thinner.
38
We illustrate this in Figure S3, where we show experimental snapshots at two different subcoolings:
∆
T
1
≈
11
◦
C
(at
t
≈
280
s
)
39
and
∆
T
2
≈
5
◦
C
(at
t
≈
1000
s
). As expected, the thickness of the crust along the fingering front appears to be thinner at
∆
T
2
.
40
Phase-field modeling
41
We develop a continuum-scale phase-field model to study gas–liquid–hydrate systems far from thermodynamic equilibrium(
5
).
42
We denote by
φ
α
the volumetric fractions of phase
α
, where
α
=
g,l,s
refers to the gas, liquid and hydrate phase, respectively.
43
At any given point in the continuum domain they satisfy:
φ
g
+
φ
l
+
φ
s
≡
1
. The system is also characterized by the pointwise
44
mole fraction of CH
4
or Xe:
χ
=
N
CH
4
or Xe
/
(
N
CH
4
or Xe
+
N
H
2
O
)
.
45
Free energy design.
We start by designing a simplified version of the Gibbs free energy functional for the three phase as a
function of
χ
and temperature (
T
):
f
l
(
χ,T
) =
ω
mix
{
χ
log(
χ
)
−
(1
−
χ
) log(1
−
a
l
(
T
)
χ
)
−
χ
log(1
−
b
l
(1
−
χ
)) +
f
l
0
}
,
[1]
f
g
(
χ,T
) =
ω
mix
{
χ
log(
χ
)
−
(1
−
χ
) log(1
−
a
g
χ
)
−
χ
log(1
−
b
g
(
T
)(1
−
χ
)) +
f
g
0
}
,
[2]
f
s
(
χ,T
) =
ω
mix
{
a
s
(
T
)(
χ
−
χ
s
)
2
+
b
s
(
T
) +
f
s
0
}
,
[3]
where
ω
mix
[J/cm
3
] is a characteristic energy density. We account for nonlinear temperature dependence of
f
α
as suggested by
46
(
6
) for gas and liquid [Eqs. Eq. (1)–Eq. (2)], and as suggested by the solidification literature (
7
–
9
) for the solid phase [Eq. Eq. (3)]:
47
a
l
=
a
l
0
/
(
T/T
c
)
4
, b
g
=
b
g
0
/
(
T/T
c
)
2
, a
s
=
a
s
0
(
T/T
c
)
and
b
s
=
b
s
0
(
T/T
c
)
, where
T
c
= 1
K is the scaling temperature.
48
Under the phase-field framework, the
f
α
’s are incorporated into the total free energy
F
(
χ,
φ
,T
)
, which also considers the
energetic interactions between phases, and is composed of the bulk free energy
f
0
and the interfacial energy (gradient squared
terms):
F
=
∫
V
[
f
0
(
χ,
φ
,T
) +
2
c
(
T
)
|∇
χ
|
2
+
2
gl
(
T
)
∇
φ
g
·∇
φ
l
+
2
gs
(
T
)
∇
φ
g
·∇
φ
s
+
2
sl
(
T
)
∇
φ
s
·∇
φ
l
+
2
g
(
T
)
|∇
φ
g
|
2
+
2
l
(
T
)
|∇
φ
l
|
2
+
2
s
(
T
)
|∇
φ
s
|
2
]
dV.
[4]
2 of 13
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
A detailed description of
f
0
in
F
and its parameters for CH
4
and Xe hydrates can be found in (
2
,
5
). By calibrating the
49
parameters in the above energy (CH
4
parameters reported in Extended Table S1), we can recover the isobaric temperature–
50
composition phase diagram predicted by existing thermodynamic equilibrium predictions (Fig.S4).
51
Evolution equations.
The proposed free energy
F
is incorporated into a phase-field model to study the nonequilibrium
thermodynamics of the three-phase system. The evolution of the system variables (
χ
and
φ
α
’s) is driven by potentials
Ψ
, which
are variational derivatives of
F
:
Ψ
χ
=
∂F
∂χ
−∇·
∂F
∂
∇
χ
,
[5]
Ψ
φ
α
=
∂F
∂φ
α
−∇·
∂F
∂
∇
φ
α
, α
=
g,l,s.
[6]
To describe the evolution dynamics, we start by imposing mass conservation of the total mixture (methane plus water):
52
∂ρ
∂t
+
∇·
(
ρ
u
) = 0
.
[7]
53
Additionally, we prescribe the conservation of mass of methane using a Cahn-Hilliard-type equation for
χ
:
54
∂ρχ
∂t
+
∇·
(
ρχ
u
)
−
R
χ
∇·
(
D
(
{
φ
α
}
)
ρ
∇
Ψ
c
)
= 0
.
[8]
55
We complete the system with a non-conserved Allen-Cahn evolution equation for
φ
g
and
φ
s
in an advective form:
56
∂φ
α
∂t
+
u
·∇
φ
α
+
R
φ
Ψ
α
= 0
.
[9]
57
Here,
R
χ
is the effective rate of diffusion and
D
(
φ
) =
φ
g
D
g
+
φ
l
D
l
+
φ
s
D
s
is a dimensionless mixture diffusion coefficient
58
(where
D
g
,
D
l
and
D
s
are normalized by a characteristic gas-phase diffusion coefficient
D
gas
). We adopt
D
g
= 1
,
D
l
= 10
−
3
59
and
D
s
= 10
−
11
(whose relative magnitudes are consistent with experimental measurements (
10
,
11
) and emulate slow diffusion
60
in liquid and extremely slow diffusion within hydrate).
R
φ
α
in Eq. 9 is the rate of phase change for phase
α
. We assume
61
isothermal conditions in our model.
62
Multiphase Hele-Shaw flow.
We assume that flow occurs predominantly parallel to the glass plates at low Reynolds number,
63
which allows us to simplify the flow description using the Hele-Shaw approximation (
12
). Instead of identifying and tracking
64
velocities within each phase, we assume a mixture velocity (at each point in the domain) that obeys Darcy’s law:
65
u
(
x,y
) =
−
k
μ
(
φ
g
,φ
l
,φ
s
)
∇
p.
[10]
66
Gas compressibility is prescribed through a gas density that is linearly dependent on pressure. Here we simplify our
67
description of phase compressibility by imposing a linear pressure dependence in density for all three phases:
68
ρ
α
=
ρ
α
0
+
c
α
t
p.
[11]
69
Finally, it is convenient for modeling to define a mixture density, which is the mass per volume at a given point in the domain.
70
Here we describe the mixture density as a
φ
α
-weighted average:
71
ρ
=
ρ
l
φ
l
+
ρ
g
φ
g
+
ρ
s
φ
s
.
[12]
72
We model all three phases as fluid with a certain viscosity. Specifically, we assume that the hydrate phase is the most
73
viscous and gas the least viscous phase. We assume that the viscosities of all phases are independent of composition and
74
normalize the values by the liquid water viscosity
μ
water
so that:
75
μ
g
= 1
/M, μ
l
= 1
, μ
s
=
M.
76
We define a mixture viscosity by blending the three viscosities by a phase-weighted average:
77
μ
(
{
φ
α
}
) =
μ
water
(
μ
g
φ
g
+
μ
l
φ
l
+
μ
s
φ
s
)
,
[13]
78
similar to what is done to diffusion coefficients in Eq. Eq. (8).
79
The hydrate phase is described as an extremely viscous fluid, whose viscosity also hardens (increases) with its age:
dθ
dt
=
φ
s
(
r
−
|
u
|
θ
D
)
,
[14]
μ
s
=
θμ
s
0
.
[15]
We introduce some variability in the hydrate growth rate along the initial gas-liquid interface, which creates a thinner
80
segment along the crust that is prone to be broken through.
81
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
3 of 13
Numerical simulations
82
To understand how hydrate grows on a quiescent interface, we first perform a series of simulations in 1D. The left portion of
83
the domain is filled with gas and the other portion filled with gas-saturated liquid water (Fig. S5a). The model successfully
84
simulates the nucleation and growth of hydrate on the initially hydrate-free gas-liquid interface (Fig. S5b). By analyzing the
85
growth dynamics, we find that the parameter
R
s
controls the rate at which hydrate grows towards its finite thickness (Fig. S5b).
86
However, the nonequilibrium steady-state thickness of hydrate is independent of the rate of hydrate formation (
R
s
) (Fig. S5b).
87
A. Gas escaping from crusted bubble..
These simulations (Fig. 3a-b) are performed in a domain (rectangular or square)
88
initially filled with partially saturated liquid and a single gas bubble with pre-existing crust. The ambient liquid is withdrawn
89
from the upper left corner at a constant rate
Q
outlet
(Fig. S6a). After breaking away from the pre-existing crust, gas flows
90
within a single hydrate-crusted finger channel towards the fluid outlet (see Supplementary Video 3). The resistance to flow
91
provided by the crust allows for flow focusing within a thin channel. The meandering behavior of the gas finger is mainly due
92
to flow resistance by the crust. However, the imposed permeability heterogeneity (Fig. S6b) also plays a role in the randomness
93
of the pattern.
94
B. Field-scale simulations..
Assuming that the bottom simulating reflector (BSR)— the interface between an underlying
95
gas reservoir and its overlying sediments— is initially flat and a hydrate layer readily forms (Extended Fig .S7), we then
96
supply methane gas to the gas reservoir through periodic recharge events with an interval of
1
/f
and a maximum flow rate
97
of
Q
in
(Fig. S6a). After each recharge events, we assume an exponential decay of reservoir pressure with time, and thus an
98
exponentially decaying profile for the gas flux at the bottom boundary. The BSR is seeded with seven locations that are prone
99
to gas breakthrough. This is achieved by assuming
R
s
= 0
at these locations so that hydrate does not form locally. In practice,
100
these could be local weak spots formed by hydrofracturing of the hydrate-bearing permeability seal (
13
), roughness of the BSR
101
(
14
) or coupled fluid flow and solid deformation that leads to decompaction weakening (
15
). The simulations can be viewed in
102
Supplementary Video 4.
103
4 of 13
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
xenon port
water port
t
(I)
(II)
(III)
(IV)
P
(MPa)
T (C)
0
10
4
2
6
8
0
40
10
20
30
(a)
(b)
25 mm
25 mm
Fig. S1. Experimental setup. a,
Birds eye view of the microfluidic flow cell with dimensions. The gap thickness is 1 mm.
b,
P-T phase diagram for Xe-H
2
O systems, indicating
coexistence of hydrate with supercritical Xe (I), hydrate with liquid Xe (II), hydrate with gaseous Xe (III) and liquid water with gaseous Xe (IV). The red arrow indicates the
trajectory imposed during experiments in the microfluidic cell.
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
5 of 13
Fig. S2. Direct visual observation of gas hydrate dissociation in two stages of the experiment.
Yellow numbers correspond to time sequence. (a) First stage: hydrate
crust melts into the surrounding liquid, creating an aqueous solution supersaturated in xenon; (b) second stage: the supersaturated solution undergoes spinodal decomposition,
creating many xenon gas bubbles (blue circles mark regions of gas exsolution).
6 of 13
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
Fig. S3.
(Left) P-T phase diagram for Xe-H
2
O systems. The red arrow indicates the trajectory imposed during our experiments in the microfluidic cell. The yellow horizontal
bars measure the subcooling at two different times during the experiment. (Right) The experimental snapshots corresponding to the two subcoolings in the left figure. Green
circles mark the fingering front at the time the images are taken.
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
7 of 13
0
0.2
0.4
0.6
0.8
1
10
20
30
40
50
60
70
80
T
L
w
-V
V
H-V
L
w
-H
L
w
H
χ
T
H-V
L
w
-H
L
w
-V
V
L
w
H
χ
0
0.2
0.4
0.6
0.8
1
(
o
C)
Fig. S4. Thermodynamic phase diagram.
The isobaric phase diagram for CH
4
and H
2
O calculated by (left) CSMGem and (right) our model.
8 of 13
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
t=2
t=0
aqueous
solution
(aq)
gas
t=0
(a)
hydrate
gas
aqueous
solution
(aq)
〈
φ
s
〉
t
(b)
0
0.5
1
1.5
2
0
0.02
0.04
0.06
0.08
R
s
=1/1.
6
R
s
=1
R
s
=1/0.
8
R
s
=1/0.
4
Fig. S5. Hydrate growth on a quiescent interface. a,
Diagram illustrating the interfacial hydrate growth problem.
b,
The thickness of interfacial hydrate as a function of time
for different growth rates (
R
s
), calculated by 1D simulations of our model on the problem illustrated in
a
.
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.
9 of 13
(a)
(b)
Q
outlet
1
1.4
1.8
2.2
2.6
Fig. S6. Setup for pore-scale simulations. a,
In a square domain, the gas bubble is initially placed at the lower right corner and liquid withdrawn from the upper left corner.
b,
Imposed heterogenous permeability field that introduces randomness to the meandering pattern.
10 of 13
Xiaojing Fu, Joaquin Jimenez-Martinez, et al.