Methane Hydrate Formation in Ulleung Basin Under Conditions of Variable Salinity: Reduced Model and Experiments

Abstract

In this paper, we present a reduced model of methane hydrate formation in variable salinity conditions, with details on the equilibrium phase behavior adapted to a case study from Ulleung Basin. The model simplifies the comprehensive model considered by Liu and Flemings using common assumptions on hydrostatic pressure, geothermal gradient, and phase incompressibility, as well as a simplified phase equilibria model. The two-phase three-component model is very robust and efficient as well as amenable to various numerical analyses, yet is capable of simulating realistic cases. We compare various thermodynamic models for equilibria as well as attempt a quantitative explanation for anomalous spikes of salinity observed in Ulleung Basin.

1 Introduction

Gas hydrate is a frozen compound in which hydrocarbons are trapped in a water molecule lattice. Gas hydrates comprise a large and dynamic carbon reservoir; see Milkov et al. (2004) and Dickens (2003). In continental margin settings with high methane concentrations, gas hydrates occur naturally in hydrate stability zone, denoted by HSZ, at water depths H greater than 300–500 meters below see level (mbsl), wherever enough methane is present. Numerous laboratory and field studies at gas hydrate-bearing sites, including several drilling expeditions in the past decades, have provided critical background data on the conditions of gas hydrate stability, and have given an overall view of the composition and distribution of gas hydrates in nature. We refer to the recent review (Collett et al. 2014) and to the monograph (Sloan and Koh 2008) for an abundant list of references which illustrate the statements above.

Gas hydrate in these systems is known to occur in conditions of extreme variations in salinity. For example, gas hydrate in Ulleung Basin (offshore Korea) occurs in formations with salinities ranging from as low as 22 practical salinity units (psu) to brines with salinity values of 82.4 psu (Torres et al. 2011). Similar large range in salinity values has been reported in naturally occurring deposits along continental margins (Torres et al. 2004, 2011). Because of the need to understand methane hydrate evolution, there is growing interest in easy and robust mathematical and computational models which can be calibrated to experimental data and account for, e.g., the variable salinity. This paper is the first of two in which we present an approximate reduced model of methane hydrate evolution in subsea sediments under conditions of variable salinity. Our two-phase three-component physical model is a simplification of comprehensive models in Liu and Flemings (2008), Garg et al. (2008), and Daigle and Dugan (2011) and is simultaneously a significant generalization of the simpler models in Xu and Ruppel (1999), Nimblett and Ruppel (2003), and Torres et al. (2004), in which simplified kinetic or even simpler mechanisms for fluid equilibria were assumed. In contrast to Torres et al. (2004) and consistently with Liu and Flemings (2008), our model fits in the general framework of multiphase multicomponent models such as those in Lake (1989) and Class et al. (2002), and implements bona fide equilibrium phase constraints known from thermodynamics (Sloan and Koh 2008; Davie et al. 2004), albeit in an approximate manner. In the companion paper Peszynska et al. (2016) we present details of numerical discretization with a particular emphasis on the variants of the time-stepping, which are enabled by the approximations proposed here.

1.1 Model Construction

Our model accounts for both transport modes of methane and of salt: advective and diffusive, and it is derived from that in Liu and Flemings (2008) under the following simplifying assumptions.

• (I) The liquid and hydrate phases are incompressible.

• (II) The pressure is fixed and is close to hydrostatic.

• (III) The temperature gradient is fixed. In particular, the energy equation is not solved and the latent heat is not accounted for.

• (IV) The depth BHSZ of bottom of HSZ is fixed and is determined either from observations, or from phase equilibria using a fixed seawater salinity value. In addition, we consider NaCl as the only inhibitor and ignore the influence of other electrolytes.

After the simplifications, our model is still rich enough to allow the study of complex dynamics of hydrate formation over thousands of years (kyr) under the conditions of variable salinity and yet is robust and very efficient compared to the published comprehensive approaches. In particular, it solves a system of two mass conservation equations for three variables, of which one is eliminated via an approximate phase equilibrium relationship. This relationship is fixed for the entire simulation, but it allows the two-way coupling between the (transport of) salt and equilibria, and therefore, the model can predict the occurrence of salinity anomalies. In contrast, the comprehensive models available to date solve four equations (mass conservation plus pressure and energy equations) for five variables and must reevaluate the phase equilibria at every grid point, time step, and at every iteration of the nonlinear solver. We acknowledge that due to the simplification following from (I–IV), the model presented here cannot be used when significant pressure or temperature changes occur. Thus, in particular, it is inadequate for simulations of gas production from hydrate.

The crux of our model rests on how the equilibrium phase behavior is formulated. The common approach in fully implicit comprehensive models is to use multivariate lookup tables for the thermodynamics constraints, and to apply variable switching (Class et al. 2002; Liu and Flemings 2008). However, the complexity and sparsity of the phase equilibrium data published in the literature makes the simulation of even simple case studies quite delicate, as we have seen in Peszyńska et al. (2010). Therefore, we use the assumptions (I–IV) and approximate the precise thermodynamics data to formulate a robust, reduced, physically consistent, phase equilibria model. We use the software CSMGem (Sloan and Koh 2008, http://hydrates.mines.edu/CHR/Software.html) and compare its results to several empirical and semiempirical algebraic approaches. These comparisons show general consistency but also differences.

Furthermore, we follow our recent work (Gibson et al. 2014; Peszynska et al. 2015) in which the methane–salinity phase behavior is realized as an (inequality) nonlinear complementarity constraint; we will refer to this elegant explicit construction as NCC-MS. NCC-MS allows to implement easily a range of models from fully comprehensive to the simpler approximate time-stepping variants in which one or more variables are assumed known. With the reduced approximate phase equilibria in the NCC-MS formulation, each part of our model can be carefully analyzed, specialized, tested, and validated, while such an endeavor is nearly impossible in the comprehensive models. In fact, rigorous analysis of the diffusive transport model of methane was first given in Gibson et al. (2014), followed by more general analysis in Peszynska et al. (2015) for the advective/diffusive transport. The NCC-MS approach enables various variants of numerical discretization and of time-stepping discussed in the second paper Peszynska et al. (2016).

Table 1 Notation and definitions

1.2 Model Application

To demonstrate the application of our model, we choose an extreme example of a site from Ulleung Basin where methane gas is known to migrate through the gas hydrate stability field and gas hydrate is present in near-seafloor sediments characterized by the presence of brine (Torres et al. 2011). We compare the model results with the data from 2010 UBGH2 expedition in which salinity spikes were observed close to the ocean floor (Kim et al. 2013). We use our model to hypothesize on what could have been the dynamics of hydrate formation that can explain these spikes. In accordance with Torres et al. (2011), we argue that large fluxes of dissolved methane cannot explain these anomalies, and the Ulleung Basin data argue against the presence of a high-salinity front as postulated by Liu and Flemings (2006).

The outline of the paper is as follows. We present the model in Sect. 2 and describe how it is calibrated using CSMGem in Sect. 3. In Sect. 4, we describe the setup of simulations and in Sect. 5 compare their results to the experimental data from 2010 UBGH2 expedition, and discuss the limitations of the current models to explain the salinity spikes. We close in Sect. 6 with conclusions. The “Appendix” provides details on some of the calculations which relate our model to that in Liu and Flemings (2008).

2 Reduced Model of Hydrate and Salinity Transport with Methane Hydrate Formation

We now describe our methane–salt transport model. The notation is summarized in Table 1. The transport takes place in the sediment reservoir $\varOmega$ under the ocean bottom; $\varOmega \subset {\mathbb {R}}^d,d=1,2,3$. Each point $x=(x_1,x_2,x_3)\in \varOmega$ is at some depth D(x) below the sea surface. In this paper, we assume that $x_3$ points in the direction of gravity upwards and that the origin $x=0$ is somewhere in, or beneath the hydrate reservoir. In the general case of a 3D reservoir, the bathymetry is variable; thus, D(x) is measured relative to the sea surface rather than to the seafloor. In 1D case, $x=x_3$, and it is customary to consider a fixed reference depth $D_\mathrm{ref}=H$ equal to the water depth H at seafloor, i.e., at the top of the reservoir.

In this paper, we assume that the conditions in $\varOmega$ are favorable for hydrate presence: i.e., the pressure is high enough and the temperature is low enough in $\varOmega$, and that there is a sufficient methane supply to the system. The latter may result from upward advection of methane gas originating at depth (Torres et al. 2011); methane may also be generated in situ via microbial methanogenesis (Hong et al. 2014). The high-pressure and low-temperature conditions are possible at large depths H, or in Arctic regions. At higher temperatures, such as those occurring at depth within the sediment, methane exist in the gas (“vapor”) phase. Upward methane transport in the gas phase has been documented, but transport in such conditions is not considered in this paper. We refer to the gas phase only when discussing phase equilibria.

The liquid and hydrate phases have respective densities $\rho _\mathrm{l},\rho _\mathrm{h}$ which, in general, are mildly dependent on the pressures and temperature, but in our model we assume (I),

$$\begin{aligned} \rho _\mathrm{l}\approx {\hbox {const}},\; \rho _\mathrm{h} \approx {\hbox {const}}. \end{aligned}$$

(1a)

Similar incompressibility assumptions are commonly made in two-phase water–oil reservoir models (Peszyńska et al. 2000; Lu et al. 2002), and (1a) is entirely reasonable over the timescale considered here.

Per assumptions (II) and (III), the pressure P(x) is usually assumed to be close to the hydrostatic pressure, and the temperature usually follows the geothermal gradient

$$\begin{aligned} P(x)\approx & {} P_\mathrm{ref} + G_\mathrm{H}(D(x)-D_\mathrm{ref}). \end{aligned}$$

(1b)

$$\begin{aligned} T(x)= & {} T_\mathrm{ref}+G_\mathrm{T}(D(x)-D_\mathrm{ref}). \end{aligned}$$

(1c)

The use of (1c) is common (Davie et al. 2004; Rempel 2012); in Peszyńska et al. (2010), we showed little influence of a particular energy model for variable T(x) on methane fluxes over long time period. Instead of (1b), one can find P(x) from the pressure equation defined in the “Appendix.”

The presence of the liquid and hydrate phase is accounted for by their void fractions, $S_\mathrm{l},S_\mathrm{h}$, respectively, also called saturations (Lake 1989; Class et al. 2002). Since $S_\mathrm{l}+S_\mathrm{h} \equiv 1$, only one of these phase saturations is an independent variable.

The liquid phase (also called “aqueous phase”) consists of water, salt, and dissolved methane components, and their corresponding mass fractions in the liquid phase are denoted by $\chi _\mathrm{lW}, \chi _\mathrm{lS}, \chi _\mathrm{lM}$, respectively. In the hydrate literature, the mass fractions $\chi _\mathrm{lM}, \chi _\mathrm{lS}$ are also called the “solubilities.” The hydrate phase is made of molecules of water and of methane, with the mass fractions denoted by $\chi _\mathrm{hW}, \chi _\mathrm{hM}$. Because of the physical nature of hydrate crystals built from a fixed proportion of methane and water molecules, it is common to assume the last two are constants, while $\chi _\mathrm{lW}, \chi _\mathrm{lS}, \chi _\mathrm{lM}$ are variables. Since for mass fractions in the same phase we have $\chi _\mathrm{lW}+\chi _\mathrm{lS} + \chi _\mathrm{lM} \equiv 1$ (Lake 1989, 2.2.8a), therefore only two of the variables $\chi _\mathrm{lW},\chi _\mathrm{lM}, \chi _\mathrm{lS}$ can be independent. In what follows, we choose the salt mass fraction $\chi _\mathrm{lS}$ and one of methane-related variables as the independent variables.

The porosity $\phi _0$ and permeability $K_0$ of sediment typically decrease with overburden pressure, i.e., with increasing D(x). If hydrate is present, then the actual porosity $\phi (x)$ available to the liquid phase is $\phi (x,t)=\phi _0(x)S_\mathrm{l}(x,t)$. The actual permeability K(x) in the presence of hydrate is an important property; however, it is only required when the pressure equation is solved.

2.1 Mass Conservation

In region $\varOmega$, we write the mass conservation equations for methane and salt components as in Liu and Flemings (2008). Each equation includes a sum of mass fractions over all phases in which a given component is present. These equations can be derived from first principles as a simplification of the comprehensive model from Liu and Flemings (2008).

$$\begin{aligned} \frac{\partial \phi _0 N_\mathrm{M}}{\partial t} -\nabla \cdot {D_\mathrm{M} \nabla {\chi _\mathrm{lM}}} + \nabla \cdot (q \chi _\mathrm{lM})= & {} f_\mathrm{M}, \end{aligned}$$

(2a)

$$\begin{aligned} \frac{\partial \phi _0 N_\mathrm{S}}{\partial t} -\nabla \cdot {D_\mathrm{S} \nabla {\chi _\mathrm{lS}}} + \nabla \cdot (q \chi _\mathrm{lS})= & {} 0. \end{aligned}$$

(2b)

Here, we have denoted by $N_\mathrm{S}$ and $N_\mathrm{M}$ the (nondimensional) concentrations of methane and salt relative to water density

$$\begin{aligned} N_\mathrm{M}= & {} S_\mathrm{l} \chi _\mathrm{lM} + R (1-S_\mathrm{l}), \end{aligned}$$

(2c)

$$\begin{aligned} N_\mathrm{S}= & {} \chi _\mathrm{lS} S_\mathrm{l}. \end{aligned}$$

(2d)

where R is a positive constant made precise below.

The flux q is the volumetric Darcy flux of the liquid phase assumed known, and the source term $f_\mathrm{M}$ is given. The diffusivities $D_\mathrm{M},D_\mathrm{S}$ are functions of $S_\mathrm{l}$

$$\begin{aligned} D_\mathrm{C}=D_\mathrm{C}^0 \phi =D_\mathrm{C}^0 \phi _0 S_\mathrm{l}, \end{aligned}$$

(2e)

where $D_\mathrm{C}^0$ is the (molecular) diffusivity of the component C in bulk brine, and $\phi _0S_\mathrm{l}$ accounts for the decrease in solubility due to the presence of porous medium (Lake 1989, 2.2–20). For components with (small) molecules of similar size, $D_\mathrm{C}^0 \approx D^0=10^{-9}{\mathrm {m}}^2$/s. We note that more complicated formulas for $D_\mathrm{C}$ involving, e.g., tortuosity, and Archie’s exponent, can be found, e.g., in Bear and Cheng (2010), Sect. 7.1C and Dullien (1979), Sect. 6.2.4.

In (2), we have four equations and five unknowns: $N_\mathrm{M},N_\mathrm{S},\chi _\mathrm{lS},\chi _\mathrm{lM}$ and $S_\mathrm{l}$. After we eliminate $N_\mathrm{M},N_\mathrm{S}$ using (2c) and (2d), we have the two mass conservation equations (2a) and (2b) with three unknowns. The additional relationship which closes the system is the phase constraint.

The quantity $\chi _\mathrm{lM}^\mathrm{max}$ determines how the methane $N_\mathrm{M}$ is partitioned between the liquid and hydrate phases. If $N_\mathrm{M}(x,t) < \chi _\mathrm{lM}^\mathrm{max}$, then only the liquid phase is present, i.e., $S_\mathrm{l}(x,t)=1$, $N_\mathrm{M} =\chi _\mathrm{lM}$, and $\chi _\mathrm{lM}$ is the independent variable which describes how much methane is dissolved in the liquid. On the other hand, when the amount present reaches the maximum amount that can be dissolved, i.e., $N_\mathrm{M} \ge \chi _\mathrm{lM}^\mathrm{max}$, the excess forms the hydrate phase with $S_\mathrm{h}=1-S_\mathrm{l}>0$. In this case, $S_\mathrm{l}$ becomes the independent variable while $\chi _\mathrm{lM}= \chi _\mathrm{lM}^\mathrm{max}$ fixed.

These constraints can be written concisely as a nonlinear complementarity constraint referred to as NCC-MS

$$\begin{aligned} \left\{ \begin{array}{ll} {\chi _\mathrm{lM}} \le \chi _\mathrm{lM}^\mathrm{max}, \; &{} {S_\mathrm{l}} =1,\\ {\chi _\mathrm{lM}} = \chi _\mathrm{lM}^\mathrm{max}, &{} \; {S_\mathrm{l}} \le 1,\\ (\chi _\mathrm{lM}^\mathrm{max}-{\chi _\mathrm{lM}} ) (1-{S_\mathrm{l}}) &{}=0 \, \end{array} \right. \end{aligned}$$

(2f)

The companion paper Peszynska et al. (2016) gives details on how (2f) is implemented in the numerical solver. Below, we discuss the data for $\chi _\mathrm{LM}^\mathrm{max}$.

The model (2) must be supplemented with boundary and initial conditions appropriate to a given case study.

2.2 Phase Behavior: Solubility Constraints

From the hydrate literature Liu and Flemings (2008) and Sloan and Koh (2008), it is known that maximum solubility constraint $\chi _\mathrm{lM}^\mathrm{max}$ depends on $P,T,\chi _\mathrm{lS}$

$$\begin{aligned} \chi _\mathrm{lM}^\mathrm{max}=\chi _\mathrm{lM}^\mathrm{max}(P,T,\chi _\mathrm{lS}), \end{aligned}$$

(3)

and there are tabulated data, or complex thermodynamics models, for $\chi _\mathrm{lM}^\mathrm{max}$. Conversely, the variables $P,T,\chi _\mathrm{lS}$ determine the circumstances in which $S_\mathrm{l} <1$ and $S_\mathrm{h}>0$, i.e., when the hydrate phase can be present. The dependence of $\chi _\mathrm{lM}^\mathrm{max}$ on the type of sediment from Daigle and Dugan (2011) will not be discussed here.

Per assumption (IV), we consider a particular approximation to (3)

$$\begin{aligned} \chi _\mathrm{lM}^\mathrm{max} \approx \chi _\mathrm{lM}^\mathrm{max}(x,\chi _\mathrm{lS}) \approx \chi _\mathrm{lM}^{\mathrm{max},0}(x)+\alpha (x){\chi _\mathrm{lS}}, \end{aligned}$$

(4)

calibrated for the case study in Ulleung Basin. To find $\chi _\mathrm{lM}^{\mathrm{max},0}(x)$ and $\alpha (x)$, we use thermodynamics models and data from the literature.

2.3 Numerical Model

The numerical model corresponding to (2) is based on a nonuniform structured grid in 1D and 2D/3D. Discretization is cell-centered finite differences (FD) with harmonic averaging and upwinding. We use operator splitting and treat advection explicitly and diffusion/equilibria implicitly, in several variants of time-stepping applied to the coupled methane–salt system. Details and sensitivity studies are provided in the companion paper Peszynska et al. (2016).

3 Model Calibration

In order to apply the model (2) to realistic cases, we need data, in particular, for $\chi _\mathrm{lM}^\mathrm{max}$ in (2f). In comprehensive models such as Liu and Flemings (2008), the data for $\chi _\mathrm{lM}^\mathrm{max}$ are provided via multivariate lookup tables based on sparse data. The sparsity contributes to the roughness of the multivariate sampling, which in turn creates difficulties for a numerical solver. These difficulties can be exacerbated by switching of the primary unknowns as in Liu and Flemings (2008), and by the use of numerical derivatives calculated from multivariate approximations, which can lead to further complications, even if the underlying case study is fairly simple.

In this section, we derive an approximate reduced model (4) for $\chi _\mathrm{lM}^\mathrm{max}$ which simplifies the phase behavior solver substantially but which honors the well-known qualitative properties of $\chi _\mathrm{lM}^\mathrm{max}$. In particular, it is known that the values of $\chi _\mathrm{lM}^\mathrm{max}$ in HSZ are most strongly controlled by the temperature (Rempel 2012; Davie et al. 2004), with only a mild dependence on salinity, and with negligible dependence on the pressure.

We also compare various theoretical and experimental approaches to provide the context for our approximation. As one of the approaches, we consider the tabulated results of CSMGem. The code CSMGem was developed by Sloan and Koh (2008), http://hydrates.mines.edu/CHR/Software.html and calculates $\chi _\mathrm{lM}^\mathrm{max}$, also called methane hydrate saturation, based on the statistical thermodynamics models proposed in Barrer and Stuart (1957), Platteeuw and Waals (1959), and Ballard (2002). CSMGem is an extension of CSMHYD which is publicly available http://hydrates.mines.edu/CHR/Software.html. Since this model is most detailed, and up to date, we select it for our numerical simulations in Sect. 4. We provide comparisons with the model by Tishchenko et al. (2005) which uses a semiempirical approach based on the theoretical work from Pitzer (1991) to derive $\chi _\mathrm{lM}^\mathrm{max}$ in conditions for $\chi _\mathrm{lS}=0$ (freshwater) to $\chi _\mathrm{lS}=2 \chi _\mathrm{lS}^\mathrm{sw}$ (twice of seawater salinity). We also consider available experimental data. Some models for $\chi _\mathrm{lM}^\mathrm{max}$ require the knowledge of methane hydrate stability pressure $P_\mathrm{eq}$. We note that in the literature, $\chi _\mathrm{lM}^\mathrm{max}$ is frequently called MHSAT, and $P_\mathrm{eq}$ is called MHEQ; we use these symbols in figures.

In practice, to get a model for $\chi _\mathrm{lM}^\mathrm{max}$, we first determine the HSZ where hydrate can coexist with liquid phase. Our main simplifying assumption (IV) is that the salinity at large depths is close to the seawater value as suggested in Davie et al. (2004). With this, we calculate the pressure $P_\mathrm{eq}$ at the three-phase equilibrium (aqueous–hydrate–vapor). The knowledge of $P_\mathrm{eq}$ fixes the depth BHSZ of the bottom of HSZ. Alternatively, as was done for Ulleung Basin, we determine BHSZ from seismic-inferred observations.

Next, above BHSZ, we only consider the two-phase aqueous–hydrate equilibria, and for this, we prepare (offline) the tabulated data on $\chi _\mathrm{lM}^\mathrm{max}$ depending on $(T,P,\chi _\mathrm{lS})$ within the range realistic for Ulleung Basin. We recognize that in some settings within the Ulleung Basin and elsewhere, there is evidence for methane transport in the gas phase within the HSZ. In this paper, however, we do not consider the gas transport. The presence of gas phase within the HSZ is the exception, and in most systems, there is no gas within the HSZ.

In general, the data for $\chi _\mathrm{lM}^\mathrm{max}$ $(T,P,\chi _\mathrm{lS})$ are trivariate. However, we can simplify further, since for a given position x within HSZ, we recall that T(x), P(x) are known. In the end, our reduced model is a fit to (4) of the tabulated data against $\chi _\mathrm{lS}$.

In this paper, we consider the stability and saturation of only structure I (sI) hydrate, with methane as the only guest component in the clathrate structure. Also, as included in assumption (IV), we consider NaCl as the only thermodynamic inhibitor. More generally, other electrolytes such as KCl or CaCl$_2$ also serve as inhibitors (Sloan and Koh 2008; Dholabhai et al. 1991); however, their effect is by an order of magnitude smaller than that of NaCl and will be neglected.

3.1 Calculation of $P_\mathrm{eq}$

The equilibrium pressure $P_\mathrm{eq}$ is the pressure at which the three phases: liquid, hydrate, and vapor, can coexist. In general, $P_\mathrm{eq}$ increases with the temperature T and decreases with the salinity $\chi _\mathrm{lS}$.

Various estimates of the dependence of $P_\mathrm{eq}$ on T and $\chi _\mathrm{lS}$ are shown in Fig. 1 including those from CSMGem, Maekawa et al. (1995), and Tishchenko et al. (2005). The model for $P_\mathrm{eq}(T,\chi _\mathrm{lS})$ from CSMGem is obtained by running CSMGem for tabulated values of $T,\chi _\mathrm{lS}$.

Fig. 1

Methane hydrate stability $P_\mathrm{eq}$ denoted by MHEQ for different salinity, pressure, and temperature estimated by various models. Available experimental data were shown for comparison. For salinity values below that of seawater, all models agree well with each other and the experimental data. The stability field estimated by Tishchenko et al. (2005) strays away from the theoretical estimation by CSMGem and from the estimation by Maekawa et al. (1995) based on the interpolation of experimental data

The algebraic model for $P_\mathrm{eq}$ from Maekawa et al. (1995) is obtained by fitting the laboratory measurements of $P_\mathrm{eq}$ with the following relationship

$$\begin{aligned} \ln \left( \frac{P_\mathrm{eq}}{P_0}\right)= & {} -926.815+\frac{31979.3}{T}+144.909 \ln (T)\nonumber \\&+\,5847.92\chi _\mathrm{lS}^\mathrm{m} + 322.026 (\chi _\mathrm{lS}^\mathrm{m}) ^2 + 5840.5 \ln (1-\chi _\mathrm{lS}^\mathrm{m}). \end{aligned}$$

(5)

Here, $P_\mathrm{eq}$ (MPa), T (K), and $P_0=0.101$ MPa are the atmospheric pressures, and $\chi _\mathrm{lS}^\mathrm{m}$ (mol/mol) is the mole fraction of NaCl in the aqueous phase. The relationship (5) is valid in conditions with salinity up to $\sim$ 8.5 times higher than seawater value and is in good agreement with laboratory data obtained under high-salinity conditions (Roo et al. 1983; Kobayashi et al. 1951).

As shown in Fig. 1, CSMGem values are close to those given by (5) and to those given by the semiempirical model from Tishchenko et al. (2005). However, for fluids with high salinity, the $P_\mathrm{eq}$ estimated in Tishchenko et al. (2005) is greater than that estimated by CSMGem and the empirical relationship derived in Maekawa et al. (1995).

However useful and accurate, the model from Maekawa et al. (1995) is not accompanied by a $\chi _\mathrm{lM}^\mathrm{max}$ model. Thus, in what follows, we use CSMGem as the model for $P_\mathrm{eq}$ with largest validity range providing both $\chi _\mathrm{lM}^\mathrm{max}$ and $P_\mathrm{eq}$.

3.2 Three-Phase Equilibrium Point(s) and the Depth $D_\mathrm{eq}$ of BHSZ

The knowledge of $D_\mathrm{eq}$ and $P_\mathrm{eq}$ and $T_\mathrm{eq}=T(x_\mathrm{eq}(t))$ is needed in the estimates of $\chi _\mathrm{lM}^\mathrm{max}$.

From (5), since $T=T(x,t)$ and $\chi _\mathrm{lS}=\chi _{lS}(x,t)$, we see that $P_\mathrm{eq}=P_\mathrm{eq}(x,t)$. If $P=P(x,t)$, then at a given time t there may be a point or points $x=x_\mathrm{eq}(t)$ at some depth $D_\mathrm{eq}=D(x_\mathrm{eq}(t))$ at which

$$\begin{aligned} x: \;\; P(x)=P_\mathrm{eq}(T(x),\chi _\mathrm{lS}(x,t)). \end{aligned}$$

(6)

In general, this means that $\chi _\mathrm{lM}^\mathrm{max}$ can vary in time t; this is allowed in the comprehensive models in Liu and Flemings (2006), Liu and Flemings (2008), Peszyńska et al. (2010), and Daigle and Dugan (2011). Further, the depth of points $x_\mathrm{eq}$ needs not be unique. These considerations must be taken into account when modeling nonhydrostatic pressure, dynamically changing temperature, and in particular when modeling the production of gas from hydrates. Unfortunately, these general considerations also make the numerical model very complex, since a recalculation of $P_\mathrm{eq}$ and $\chi _\mathrm{lM}^\mathrm{max}$ must be done at every point, at every time step, and/or even within every iteration of an iterative solver. Furthermore, if $\chi _\mathrm{lM}^\mathrm{max}$ varies in time, the model is not amenable to even the general mathematical analysis of well-posedness in Peszynska et al. (2015).

However, in basin modeling, it is reasonable to make some approximations. Following the main assumptions (II, III) we adopted, with hydrostatic pressure and a linear temperature profile as in (1b) and (1c), we see that P, T are monotone in x. If, in addition, the salinity $\chi _\mathrm{lS}\approx const$, there is at most one such depth $D_\mathrm{eq}$ where (6) holds; this is the base of HSZ. For depths above $D_\mathrm{eq}$ (or temperatures lower than $T_\mathrm{eq}$), liquid in $\varOmega$ can coexist with hydrate phase.

If the salinity within HSZ is nonconstant, the conundrum is that we do not know $\chi _\mathrm{lS}(x)$ when calculating $D_\mathrm{eq}$ from (6). However, we can assume, as suggested in Davie et al. (2004), that the salinity at the depths close to $D_\mathrm{eq}$ equals that of $\chi _\mathrm{lS}^\mathrm{sw}$. This means that the base $D_\mathrm{eq}$ of HSZ is calculated only once and is fixed; we identify BHSZ as the set of points $x_\mathrm{eq}$ for which

$$\begin{aligned} x_\mathrm{eq}: \;\; P(x_\mathrm{eq})=P_\mathrm{eq}(T(x_\mathrm{eq}),\chi _\mathrm{lS}^\mathrm{sw}). \end{aligned}$$

(7a)

This approximation is clearly reasonable given the fact that it only determines BHSZ.

Alternatively, one may have additional information about $D_\mathrm{eq}$ from the seismic-inferred depth of hydrate stability zone. Such was the case of Ulleung Basin where we know the depth of BHSZ (Table 2).

Table 2 Range of validity of $P_\mathrm{eq}$ and $\chi _\mathrm{lM}^\mathrm{max}$ models in Sect. 3.3.1

3.3 Model for $\chi _\mathrm{lM}^\mathrm{max}$

Once we know $D_\mathrm{eq}$, the values $P_\mathrm{eq}$ and $T_\mathrm{eq}$ are fixed. With these, one calculates the maximum methane mass fraction at the three-phase equilibrium, which is used in turn to get $\chi _\mathrm{lM}^\mathrm{max}(T(x),\chi _\mathrm{lS}(x,t))$ at a given x, t.

We recall first the parametric model from Davie et al. (2004) which provides a linear fit to data generated by the theoretical thermodynamics calculations from Zatsepina and Buffett (1997); see also Table 1 in Davie et al. (2004). The model

$$\begin{aligned} C_3(T,P)= & {} C_3(T_0,P_0)\nonumber \\&+\, \partial _\mathrm{T} C_3(T_0,P_0) (T-T_0) +\partial _P C_3(T_0,P_0) (P-P_0), \end{aligned}$$

(7b)

provides the solubility of methane at the three-phase equilibrium point Davie et al. (2004) based on an estimate of $C_3, \partial _\mathrm{T} C_3, \partial _P C_3$ at some given $(T_0,P_0)$. We provide these for completeness in Table 3.

In particular, knowing $D_\mathrm{eq},P_\mathrm{eq},T_\mathrm{eq}$, we can calculate from (7b) the solubility $C_3(T_\mathrm{eq},P_\mathrm{eq})$ at the base of HSZ. To correct for the influence of salinity, and to find $\chi _\mathrm{lM}^\mathrm{max}$ at a given depth D(x) within HSZ, we follow Davie et al. (2004) and use

$$\begin{aligned} C_\mathrm{eq}(T(x),\chi _\mathrm{lS}) = C_3(T_\mathrm{eq},P_\mathrm{eq}) exp\left( \frac{T(x)-T_\mathrm{eq}}{a}\right) (1-\beta \chi _\mathrm{lS}^\mathrm{M}). \end{aligned}$$

(7c)

Here, $a=14.4\;{\mathrm {K}}$, $\beta =0.1 \;{\mathrm {mol}}^{-1}$ are the parameters determined from the theoretical calculation of (Zatsepina and Buffett 1997), Eq. (7). The variable $\chi _\mathrm{lS}^{rm M}$ is the salinity in the unit of molality. See also Rempel (2012, Eq. 11), for (7c) calculated for pure water in heterogeneous sediments. Finally, we obtain $\chi _\mathrm{lS}^\mathrm{max}$ via the conversion factor

$$\begin{aligned} \chi _\mathrm{lM}^\mathrm{max}({\mathrm {kg/kg}}) = C_\mathrm{eq} ({\mathrm {mM}}) 10^{-3} \frac{16.04}{1030}. \end{aligned}$$

(7d)

Here, we have used molecular weight of methane equal 16.04 g/mol, the seawater density $1030\ {\mathrm {g/L}}$, and recalled that $1 {\mathrm {mM}}=10^{-3}\ {\mathrm {mol/L}}$.

Combining (7d) with (7c), we see that the dependence of $\chi _\mathrm{lM}^\mathrm{max}$ on $\chi _\mathrm{lS}$ is linear, which is consistent with the model postulated in (4).

We compare the model (7d) and various other parametizations and experiments of $\chi _\mathrm{lM}^\mathrm{max}$ including CSMGem, (Tishchenko et al. 2005; Sloan and Koh 2008; Davie et al. 2004; Kim et al. 2008) in Fig. 2. Estimates using freshwater and low pressure in Tishchenko et al. (2005) and Sloan and Koh (2008) agree well with each other and with experimental results. As salinity increases, the estimates from both Tishchenko et al. (2005) and Davie et al. (2004) suggest a reduction in $\chi _\mathrm{lM}^\mathrm{max}$ (i.e., the reduction in the maximum methane mass fraction in equilibrium with hydrate), in agreement with the laboratory results from Kim et al. (2008). CSMGem, however, suggests an increase in $\chi _\mathrm{lM}^\mathrm{max}$, consistent with the theoretical calculation of Zatsepina and Buffett (1998), which also suggest an increase in $\chi _\mathrm{lM}^\mathrm{max}$ at salinities higher than about 0.1 mol/kg of water, or 7 g/kg. Finally, since only few experimental data for high salinity are available (Kim et al. 2008), the evaluation of accuracy of theoretical analyses for high salinity is difficult.

We remark that, if the position $D_\mathrm{eq}$ of BHSZ changes, one should recalculate $C_3$ and $C_\mathrm{eq}$ in (7b) and (7c). This is done in comprehensive models, but has not been included in our model.

Table 3 Parameters required in Eq. (7b) to calculate methane hydrate stability and saturation following Davie et al. (2004)

Fig. 2

Methane hydrate saturation $\chi _\mathrm{lM}^\mathrm{max}$ for different salinity, pressure, and temperature estimated by various models. Note that only few experimental data for pure water and $\chi _\mathrm{lS}\approx 2\chi _\mathrm{lS}^\mathrm{sw}$ are available. The value $\chi _\mathrm{lM}^\mathrm{max}$ estimated by CSMGem is always higher than the one estimated by Davie et al. Davie et al. (2004), while $\chi _\mathrm{lM}^\mathrm{max}$ estimated by Tishchenko et al. Tishchenko et al. (2005) overlaps with one or the other approaches. The experimental data are from Lu et al. (2008) and Kim et al. (2003). More figures are available as supplementary material

3.3.1 Use of CSMGem to Get $\chi _\mathrm{lM}^\mathrm{max}$

First, for a given T(x), we calculate $P_\mathrm{eq}$. Then, we use P(x) to find the depth $D_\mathrm{eq}$ of BHSZ assuming seawater salinity at BHSZ. Next, we use CSMGem to estimate $\chi _\mathrm{lM}^\mathrm{max}$. We first construct a lookup table in which the input values of pressure $P_i$, temperature $T_j$, and salinity $\chi _{\mathrm{lS},k}$ cover the range of interest. For the pressures, we consider the range between the seafloor pressure and that at BHSZ. Since pressure has relatively small effect on $\chi _\mathrm{lM}^\mathrm{max}$, we only use these two values $P_1=P_\mathrm{ref}$ and $P_2=P_\mathrm{BSR}$ as the grid points. The temperature dependence is very significant, and we consider the interval $T_j \in (273\,{\mathrm {K}}, 291\,{\mathrm {K}})$, with $\varDelta T = 2$ K. We also consider salinity values $\chi _{\mathrm{lS},k} \in (0, 0.125) {\mathrm {kg/kg}}$, where the right end point is four times the seawater salinity $\chi _\mathrm{lS}^\mathrm{sw}$, with $\varDelta \chi _\mathrm{lS} = 0.0156$ for the total of nine grid points.

Next, we use CSMGem to estimate $\chi _\mathrm{lM}^\mathrm{max}$ for each of the grid points $(P_i,T_j,\chi _{\mathrm{lS},k})$. This is done by trial and error: We provide CSMGem with some guess of $\chi _\mathrm{lM}$, and CSMGem predicts the phase conditions for $(P_i,T_j,\chi _{\mathrm{lS},k}, \chi _\mathrm{lM})$. We try different values of $\chi _\mathrm{lM}$ until we locate the maximum methane mass fraction $\chi _\mathrm{lM}^\mathrm{max}\vert _{(P_i,T_j,\chi _{\mathrm{lS},k})}$ for which methane is only in two phases, i.e., as dissolved methane and methane hydrate. This process gives a table of values

$$\begin{aligned} (P_i,T_j,\chi _{\mathrm{lS},k},\chi _\mathrm{lM}^\mathrm{max}\vert _{(P_i,T_j,\chi _{\mathrm{lS},k})}) \end{aligned}$$

with $1\le i\le 2, 1\le j\le 20,1\le k \le 9$.

Next, for each grid point $(P_i,T_j)$, we estimate the regression between $\chi _{\mathrm{lS},k}$ and $\chi _\mathrm{lM}^\mathrm{max}\vert _{(P_i,T_j,\chi _{\mathrm{lS},k})}$. The regression provides us, for each $(P_i,T_j)$ in the gridded table, with the coefficients $A_{ij}$ and $B_{ij}$ of the linear model so that

$$\begin{aligned} \chi _\mathrm{lM}^\mathrm{max}\vert _{(P_i,T_j,\chi _{\mathrm{lS},k})} = A_{ij} + B_{ij} \chi _{\mathrm{lS},k}. \end{aligned}$$

As shown in Fig. 3, the values $A_{ij},B_{ij}$ are not very sensitive to the pressure; thus, we approximate further

$$\begin{aligned} \chi _\mathrm{lM}^\mathrm{max}\vert _{(P_i,T_j,\chi _{\mathrm{lS},k})} = \bar{A}_{j} + \bar{B}_{j} \chi _{\mathrm{lS},k}. \end{aligned}$$

(8)

where $\bar{A}_j=A_{1j}$, $\bar{B}_j=B_{1j}$.

Fig. 3

Dependence of $\chi _\mathrm{lM}^\mathrm{max}$ on the salinity $\chi _\mathrm{lS}$ estimated from CSMGem as in (9). Positive value of the slope $\alpha$ indicates that the methane hydrate is more difficult to form at higher salinity. Higher temperature elevates $\chi _\mathrm{lM}^\mathrm{max}$ and makes methane hydrate more difficult to form. $\chi _\mathrm{lM}^\mathrm{max}$ decreases only slightly when increasing pressure at the same temperature and salinity

For the cases where $P(x), x \in \varOmega$ changes by more than 1-2 MPa, one may consider a more accurate multivariate model than (8).

In the last step, we connect (8) to the model (4). With constant geothermal gradient according to (1c), each $T_j$ corresponds to a unique depth $D_j$; thus, we set up a lookup table extending (8) to

$$\begin{aligned} \chi _\mathrm{lM}^\mathrm{max}(x,t) = \bar{A}(D(x)) + \bar{B}(D(x)) \chi _\mathrm{lS}(x,t) \end{aligned}$$

(9)

where $\bar{A}(x), \bar{B}(x)$ are the appropriate piecewise linear functions built from $(D_j,A_j)$ and $(D_j,B_j)$, respectively.

We observe that there is qualitative agreement between the approximation (4), the parametric model (7c), and the regression (9). In particular, $\chi _\mathrm{lM}^{\mathrm{max},0}(x):=\bar{A}(D(x))$ is the temperature dependent coefficient, and $\alpha (x) :=\bar{B}(D(x))$ is the salinity-dependent coefficient.

However, we note that $\bar{B}$ found from tabulated data can have any sign. In fact, we find that it is positive, in contrast to the model (7c). In the companion paper Peszynska et al. (2016) we discuss the sensitivity of simulations to the particular value and to the sign of $\bar{B}$.

4 Application to the Ulleung Basin Case

In this section, we describe how the model (2) was calibrated using data from Ulleung Basin. The second drilling expedition to the Ulleung Basin (UBGH2) (Bahk and Kim 2013) offshore South Korea (Fig. 4) drilled four sites that targeted the acoustic blanking chimneys in the seismic reflection data. These acoustic features extend from below the HSZ to near the seafloor, where they are usually accompanied by the presence of pockmarks or mounds on the seafloor bathymetry (Horozal et al. 2009). The seismic blanking zones have been interpreted to image conduits for gas migration, because of the low impedance of seismic waves as they travel through gas. Gas hydrates with different modes of occurrence were recovered from all four sites. From three of the sites (UBGH2-3, UBGH2-7, UBGH2-11), massive gas hydrates related to fracture-filling (or grain displacing) morphology were observed at depths shallower than 6 mbsf (Bahk and Kim 2013). Disseminated gas hydrates related to either fracture-filling or pore-filling modes were recovered from UBGH2-2_1 (Bahk and Kim 2013). Finally, the porosity values were found to be

$$\begin{aligned} \phi _0 \in (0.6, 0.87), \end{aligned}$$

(10)

with a few local anomalies down to 0.4. In the simulations, we use the actual nonconstant porosity data for this site.

Fig. 4

a Bathymetric map of the Ulleung Basin, offshore Korea, showing the location of the four drilled sites in this study. b Seismic profiles across the drilled sites. The rectangles are the location and the depth covered by drilling. These locations were chosen inside the seismic blanking zone in the chimneys. (From KIGAM)

4.1 Salinity Data

Salinity, pressure, and temperature conditions are fundamental in constraining the stability of gas hydrate. For the Ulleung Basin, the salinity data obtained shipboard are of less precision than dissolved chloride (Scientists 2010). We therefore use the chloride data and convert it to salinity (see Fig. 5) using the empirical relationship obtained by fitting all data from UBGH2 sites with

$$\begin{aligned} S=61.6 {\hbox {Cl}}_\mathrm{M} + 1.4301, \end{aligned}$$

(11)

where ${\hbox {Cl}}_\mathrm{M}$ (M) is chloride concentration in mol/L (M).

Fig. 5

Left profiles of Cl in pore water for the four study sites. Right empirical relationship between salinity and chlorinity for UBGH2 data. The arrow denotes seawater values. The freshening observed at depth precludes inferences for the presence of a high-salinity front as the mechanism to support methane transport in the gas phase as postulated by Liu and Flemings (2008)

Pore water chloride profiles from these four sites reflect gas hydrate kinetics and fluid origins; see Fig. 5. Chloride concentrations at the bottom of the recovered sections are always lower than seawater, which have been interpreted as reflecting input of freshwater from clay mineral dehydration reactions at depth (Kim et al. 2013). The shallower sediment sections show different degrees of chloride enrichment at each site. At UBGH2-3, we have the most prominent chloride peak, with concentrations almost three times the seawater value. At UBGH2-11 and UBGH2-7, the enrichments range from a few millimolar to $\approx$180 M above seawater concentration, respectively. The site UBGH2-2_1 shows the strongest signal of deep-sourced freshwater input, but has no enrichments in chloride. It is worth noticing that these enrichments in chloride concentration are minimum values, since they may be affected by gas hydrate dissociation during core recovery (Scientists 2010). We also refer to Malinverno et al. (2008) for more discussion.

Finally, the salinities we infer from shipboard measurements were not measured in situ. Indeed, the “real” Cl and salinity are likely higher than what we measured, because of gas hydrate dissociation during recovery lowers the pore fluid salt and chloride concentrations, but none of our conclusions drawn in Sect. 5 will be different if using the “real” salinity.

4.2 Temperature and Pressure Data

The data from Ulleung Basin include temperature $T_\mathrm{ref}$ at the seafloor and downhole temperature measurements from which we estimate $G_\mathrm{T}$, see Table 4.

Further, with known hydrostatic gradient $G_\mathrm{H}$, the pressure at the seafloor, the pressure at the first gas hydrate appearance, and the pressure at the base of the HSZ are listed in Table 4. In a typical reservoir of thickness of 100 to 200 m, the pressure difference in the hydrostatic distribution is about $\varDelta P_\mathrm{H} \le 2$ MPa, and it significantly exceeds the contributions to pressure difference that may occur due to advective fluxes that have been observed. Thus, it makes sense to assume hydrostatic relationship (1b).

Table 4 Basin parameters of the four study sites in Ulleung Basin

5 Model Results and Discussion

In this section, we apply our model to the case from UBGH2-7 in an effort to illustrate the applicability of the model to a natural system and to explain the coupled methane and salinity dynamics resulting in salinity spikes accompanying hydrate deposits. We provide background with motivation, details on the setup of the cases, and we discuss the results.

Background. Based on purely thermodynamic considerations, water and gas hydrate will coexist in the sediment section that lies within the HSZ. As the temperature in the sediment increases according to the attendant geothermal gradient, a depth is reached where gas hydrate becomes unstable. Below this depth, water and free gas coexist, but as long as there is water available in the formation, free gas should not be present within the HSZ. There is however ample evidence of methane migration through the HSZ at gas hydrate provinces worldwide. Observations of methane discharge at the seafloor, pressure core sampling imaging, and analyses of methane concentrations at in situ pressures, acoustic blanking in seismic data, and logging data all support the vertical migration of gas through the HSZ, which in most cases result in the formation of massive gas hydrate deposits at or near the seafloor (Torres et al. 2011).

The report of the presence of near-surface brines associated with massive gas hydrate deposits on Hydrate Ridge (Oregon) led to the development of hypotheses to explain this observation. Torres et al. (2004) used a one-dimensional transient model to simulate the observed chloride enrichment and show that in order to reach the observed high chloride values, methane must be transported in the gas phase from the depth of the BSR to the seafloor. Methane transport exclusively in the dissolved phase is not enough to form methane hydrate at the rates needed to generate the observed chloride enrichment. As shown by Trehu et al. (2004), when enough free gas accumulates below the HSZ, the excess (nonhydrostatic) pressure at the top of the gas layer may be sufficient to fracture the sediments and drive gas toward the seafloor. Alternatively, Liu and Fleming argue in Liu and Flemings (2006) that as gas migrates from below the HSZ, gas hydrate formation depletes water and elevates salinity enough to shift the local three-phase equilibrium to the point where the aqueous water, hydrate, and vapor (free gas) coexist, thus allowing vertical migration of free gas through the HSZ. The role of salinity in the thermodynamics of hydrate is important here, since there is a $1.1\,^{\circ }{\mathrm {C}}$ offset in dissociation temperature of methane hydrate in 33%NaCl, relative to that for pure water. Rapid increase in salinity due to recent gas hydrate formation poses negative feedback on hydrate crystallization by shifting the phase boundary.

There have been additional observations of pore fluids highly enriched in dissolved chloride at sites of massive gas hydrate occurrence in northern Cascadia accretionary margin (Canada), the Krishna–Godavari Basin (India), and the Ulleung Basin (Korea). The sites drilled on seismic acoustic chimneys indicative of free gas transport in the Ulleung Basin all show chloride enrichments of up to 1440 mM from near-seafloor to depths of  100 meters below seafloor (mbsf). Below the depth of chloride maxima, however, chloride values approach concentrations that are lower or equal to seawater values, with minor negative chloride anomalies superimposed on baseline that reflect discrete gas hydrate-bearing horizons (Torres et al. 2011). None of these sites, however, show any evidence for the elevated salinity values beneath the shallow lens of massive hydrate formation (Torres et al. 2011). Extreme high-salinity values (of up to 3 times seawater values (Liu and Flemings 2008) have been postulated by current models, as these high values are needed to create a shift in the gas hydrate thermodynamic equilibrium and sustain gas transport from the base of the gas hydrate stability front to the seafloor.

Below we apply our model in an effort to explain the observed salinity anomalies. It turns out that we are only partially successful.

Model setup. We use fully implicit numerical solver implementation of (2) with d$x=1$m, d$t=1$ year; see details in Peszynska et al. (2016).

The data from UBGH2-7 are along the vertical transect, and thus, the case is essentially 1D, and we set up $\varOmega =(0,L)$ where $L=124$ m is the reservoir thickness. The bottom of the reservoir is at $x=0$ and is at BHSZ. We assume T and P as in Sect. 4.2. We use relatively small advective flux q, and thus, solving pressure equation is not necessary.

We set up the following boundary and initial conditions. The boundary conditions for methane and salt components are needed at $x=0$ and $x=L$. For the top of reservoir $x=L$, i.e., sea bottom, we use seawater salinity and zero methane concentrations

$$\begin{aligned} \chi _\mathrm{lM}(L,t)=0, \;\; \chi _\mathrm{lS}(L,t) = \chi _\mathrm{lS}^\mathrm{sw}. \end{aligned}$$

(12)

At $x=0$, we assume conditions above BHSZ and set up boundary condition for methane to be given by $\chi _\mathrm{lM}^\mathrm{max}$ at the corresponding depth. For salinity at $x=0$, we use the observed salinity values $\chi _\mathrm{lS}^0=0.0273$ kg/kg shown in Fig. 5 following (Kim et al. 2013)

$$\begin{aligned} \chi _\mathrm{lM}(0,t)=\chi _\mathrm{lM}^\mathrm{max}(0,\chi _\mathrm{lS}^{0}), \;\; \chi _\mathrm{lS}(0,t) = \chi _\mathrm{lS}^{0}. \end{aligned}$$

(13)

The initial conditions are

$$\begin{aligned} \chi _\mathrm{lM}(x,0)=0, \;\; \chi _\mathrm{lS}(x,0) = \chi _\mathrm{lS}^I(x), \end{aligned}$$

(14)

where $\chi _\mathrm{lS}^I(x)$ is a linear function between $\chi _\mathrm{lS}(0,0)$ and $\chi _\mathrm{lS}(L,0)$.

We use reservoir parameters listed in Table 4 and set up five different scenarios to investigate how the profiles of dissolved methane concentration, salinity, and gas hydrate saturation respond to different modes of aqueous fluid transport. The cases are summarized in Table 5.

Table 5 Parameters of the five simulation cases

5.1 Scenarios with Different Advection Rates and Sources

Cases 1, 2, and 3 compare simulation scenarios with different Peclet numbers as in Fig. 6, 7, and 8. Advection transports the fluids with abundant methane from sources below HSZ, which facilitates the formation of hydrate, see Fig. 8. With a strong advective flux (Case 3), gas hydrate saturation reaches more than 30 % after 100 kyr of simulation. This is in contrast to Cases 1 and 2 with Peclet numbers smaller or equal to 1.

Fig. 6

Model results of Case 1 for the fluid system with a small Peclet number (4E-6). In this case, diffusion alone is not sufficient to deliver enough methane to form gas hydrate

Fig. 7

Model results of Case 2 for the fluid system with Peclet number close to 1. Even though the advection component is stronger in this case compared to Case 1, still not enough methane delivered for gas hydrate formation within the simulation time

Fig. 8

Model results of Case 3 for the fluid system with larger Peclet number. Methane is rapidly delivered by advection to form the gas hydrate in the entire sediment column. Salinity, however, decreases due to the effective delivery of fresh fluid from the bottom. This salinity trend is different from the observations

However, even with very strong advection in Case 3, no brine is formed at any depth in the sediments. On the contrary, due to the strong fluid advection prescribed in this scenario, the whole sediment column is flushed with the freshwater. Such result contradicts the observations from our study sites, where shallow brine coexists with the abundant gas hydrate in the sediments in the upper 100 mbsf as in Fig. 5. A similar case study applied in Torres et al. (2004) to Hydrate Ridge led the authors to conclude that the methane transport exclusively by advection is not sufficient to sustain the hydrate formation rate required to produce the observed salinity enrichment. A different source of methane other than aqueous transport from depth was postulated in Torres et al. (2004) to be required.

In Case 4, we postulate therefore the existence of a source of methane $f_\mathrm{M} \ne 0$ in the sediment section where abundant gas hydrate was observed (17 mbsf at UBGH2-7). In this simulation, we use minimum advective flux (Peclet number $\ll 1$ as in Table 5) and show that in response to the strong methane input, gas hydrate saturation exceeds the highest saturation obtained in Case 3 within 5 kyr. Because of the rapid formation of methane hydrate, dissolved ions accumulate in the pore fluids faster than are lost by diffusion to the overlying bottom water, leading to a brine patch above 50 mbsf. After running the model for 10 kyr, the hydrate saturation exceeds 60 % and the salinity is 1.5 times higher than $\chi _\mathrm{lS}^\mathrm{sw}$ in bottom seawater, a value that is similar to what we observed in the pore water profiles in Fig. 9.

Fig. 9

Model results of Case 4 for the fluid system with a small Peclet number and a methane source term at 25 mbsf. The observed salinity enrichment is similar to that for experimental data by adding an arbitrary source of methane. The source term contributes large quantity of methane in a short time sufficient for rapid hydrate formation which in turn creates the salinity spike. Due to the insignificant advection component assigned in this case, diffusion is not strong enough to erase such salinity spike

In Case 5 shown in Fig. 10, we include both large advective flux q and an arbitrary methane source $f_\mathrm{M} \ne 0$. Similarly as in Case 4, gas hydrate saturation increases rapidly around the depths where methane source is present. However, the salinity enrichment in Case 5 is different than that observed in Case 4: The highest value is smaller, and the profile is nonsymmetric because some of the salt is transported toward the seafloor by strong fluid advection.

Fig. 10

Model results of Case 5 and the fluid system with a large Peclet number and a methane source term. Similar to the results from Case 4, large quantity of gas hydrate forms in less than 2000 years. The salinity enrichment is smaller compared to Case 4, and its profile is nonsymmetric due to the more effective fluid transport by larger advection

We note that in Cases 4 and 5, one might argue that pressure Eq. (21) should be solved to account for the local value of $\nabla \cdot q = f=f_\mathrm{M}$ instead of assuming $\nabla \cdot q=0$. However, the methanogenesis represented by $f_\mathrm{M}$ turns carbon from solid phase (organic matter) to dissolved phase (dissolved methane in pore fluid) and does not introduce new carbon into the overall system, thus $f=0$.

5.2 Discussion

The model (2) appears to reproduce the two-way coupled dynamics, and the hydrate and salinity profiles, in a manner consistent with the intuition. Furthermore, Case 4 gives results which are close to the profiles recorded in experiments. However, the presence of large source of methane $f_\mathrm{M}$ is needed to create the shallow brine patches, and the magnitude of $f_\mathrm{M}$ is not fully explained.

5.2.1 Limitations of the Model in Its Ability to Explain the Experimental Data

As shown in Hong et al. (2014), microbial methane production through organic matter degradation initiates at the depth where sulfate in the pore water is depleted and methane concentration starts to increase, i.e., in sulfate–methane transition zone (SMTZ). The depth of microbial methane production may correspond to the location of the brine patches observed in Ulleung Basin. Therefore, in Case 4, we tested whether in situ methanogenesis could provide the methane required to sustain the rapid hydrate formation. Methanogenesis rates in Ulleung Basin, estimated from one chimney and one nonchimney site using a kinetic model constrained by pore water data, range from a few to $\approx$25 mmol/m$^3$/year Hong et al. (2014). Using the unit conversion (22), we see that the rate $f_\mathrm{M}$ assumed in Case 4 is significantly higher than the realistic rate of methanogenesis estimated in Hong et al. (2014). In other words, the rate $f_\mathrm{M}$ proposed in Hong et al. (2014) is not large enough to account simultaneously for rapid gas hydrate formation and the associated shallow brine observed in Ulleung Basin.

As another possible explanation, one might argue that there might be a lateral advective transport of gas which might provide the source of methane. However, the seismic and chemistry analyses presented in Kim et al. (2013), Hong et al. (2014) suggest that most of the methane is generated below the SMTZ, or even deeper, and move upward as imaged in seismics, with no lateral advection.

Thus, while the simulation gives results consistent with the data, further hypotheses are needed to explain the observations.

5.2.2 Inclusion of Gas Phase

Similarly to the reasoning used in Torres et al. (2004) for the Hydrate Ridge case, we are led to conclude that the methane in the Ulleung Basin sites discussed here must be advecting in the gas phase from below the model domain. The methane solubility is too low for fluid advection to supply enough methane, with advection rate slow enough not to erase the positive salinity lense. Most likely, there is a source of gas below the HSZ, as imaged in seismic data, but free gas cannot travel through HSZ in the model (2) nor in the comprehensive models (Liu and Flemings 2008) since these assume that water is abundant. Liu and Flemings in Liu and Flemings (2006) hypothesized that the positive salinity anomaly that results from rapid hydrate formation at the base of the HSZ sustains a local three-phase equilibrium that allows methane gas to migrate upward and extends the saline tongue throughout the HSZ. Such extended positive salinity anomaly is, however, not observed in Ulleung Basin. Rather, the observed profiles as in Fig. 5 show that the brine is confined to shallow depths less than 50 mbsf, and to salinities lower than seawater salinities at depths greater than that.

5.3 Salinity Dependence

Furthermore, according to the $P_\mathrm{eq}$ calculations for UBGH2-7 in the pressure range (21.65–22.90) and the temperature range (273.55–294.8), at the depth of the salinity spikes between 20–30 mbsf, we cannot have free gas phase, even if salinity equals double the seawater value. Therefore, our conclusions from simulation results are not affected by the particular approximations made to obtain the reduced phase equilibria model. In addition, the difference in salinity data that can be attributed to the measurements shipboard versus in situ does not change our conclusions.

Further extensions of the model, and in particular the inclusion of methane transport in the gas phase, are therefore needed to explain the particular salinity spikes and are outside the present scope.

6 Conclusions

In this paper, we presented a reduced model of transport of methane and salt dissolved in liquid phase, with accompanying methane hydrate formation. The model was obtained from the comprehensive model in Liu and Flemings (2008) after several simplifying assumptions were made. These assumptions are easily justified for basin modeling and make our reduced model very compact, efficient, and easily amenable to the various analyses. The model is easily calibrated using phase behavior described in the literature, and we described in detail good agreement between various empirical and algebraic models. Thus, our paper provides a bridge between the practical and useful models and the rigorous mathematical model and computational analyses, and thus represents a useful tool for modeling the dynamic gas hydrate evolution in marine systems. In addition, it opens the door to various new computational simulations while it can be calibrated with the experimental data.

We were able to obtain good quantitative agreement between the model results and the data from Ulleung Basin by providing an additional methane source within the modeled domain, but note that in situ methanogenesis is not sufficient to generate the needed methane. In addition, the presence of fresh fluids at depth in Ulleung Basin sites that host near-seafloor brine patches argues against the development of a large positive salinity anomaly rising from the base of the HSZ to the seafloor, which could support methane transport through the HSZ as proposed by others, e.g., Liu and Flemings (2006). Our results are consistent with previous work by Torres et al. (2004) and Torres et al. (2011). However, we stress that since neither in situ methanogenesis nor transport within a salinity front are consistent with Ulleung Basin data, there must be a separate process supplying enough methane, so that the salinity spikes that accompany near-surface gas hydrate patches can be sustained.

Appendix

For completeness, we recall here the details on modeling as well as certain auxiliary conversion factors.

1.1 Derivation of Reduced Model

The conservation of mass for the methane component in hydrate zone (Liu and Flemings 2008) takes the form

$$\begin{aligned} \frac{\partial }{\partial t} \left( \phi _0 (S_\mathrm{l} \rho _\mathrm{l} \chi _\mathrm{lM} + S_\mathrm{h} \rho _\mathrm{h} \chi _\mathrm{hM})\right) + \nabla \cdot (q\rho _\mathrm{l} \chi _\mathrm{lM} ) - \nabla \cdot (D_{lM} \rho _\mathrm{l} \nabla \chi _\mathrm{lM} )= & {} \bar{f}_\mathrm{M}. \end{aligned}$$

(15)

In this equation, $\bar{f}_\mathrm{M}$ is an external source of methane, e.g., due to bacteria-induced methanogenesis.

The accumulation part (the term under the time derivative) can be rewritten with $N_\mathrm{M}= \frac{\bar{N}_\mathrm{M}}{\rho _\mathrm{l}}$ as

$$\begin{aligned} N_\mathrm{M} \rho _\mathrm{l} = \bar{N}_\mathrm{M} = S_\mathrm{l} \rho _\mathrm{l} \chi _\mathrm{lM} + S_\mathrm{h} \rho _\mathrm{h} \chi _\mathrm{hM} =S_\mathrm{l} \chi _\mathrm{lM}+(1-S_\mathrm{l})R. \end{aligned}$$

where the (dimensionless) quantity R is

$$\begin{aligned} R :=\frac{\rho _\mathrm{h} \chi _\mathrm{hM}}{\rho _\mathrm{l}} \approx {\hbox {const}}, \end{aligned}$$

(16)

Note that $\bar{N}_\mathrm{M}$ is the total mass of methane per unit volume which accounts for the methane present both in the liquid and hydrate phases, and $N_\mathrm{M}$ is its dimensionless counterpart, relative to brine density.

Finally, it is useful to see that $S_\mathrm{l}(x,t;N_\mathrm{M})$ is a function

$$\begin{aligned} S_\mathrm{l}= \frac{N_\mathrm{M} - R}{\chi _\mathrm{lM} - R}=\left\{ \begin{array}{cc} 1,&{}\; N_\mathrm{M} \le \chi _\mathrm{lM}^\mathrm{max}(x,t),\\ \frac{N_\mathrm{M} - R}{ \chi _\mathrm{lM}^\mathrm{max} (x,t)-R},&{}\; N_\mathrm{M} > \chi _\mathrm{lM}^\mathrm{max}(x), \end{array} \right. \end{aligned}$$

(17)

Upon ${f_\mathrm{M}} :=\frac{\bar{f}_\mathrm{M}}{\rho _\mathrm{l}},$ and rescaling, we rewrite (15) in the form (2a).

Next, mass conservation for salt has the form

$$\begin{aligned} \frac{\partial }{\partial t} (\phi _0 (S_\mathrm{l} \rho _\mathrm{l} \chi _\mathrm{lS} )) + \nabla \cdot (q \rho _\mathrm{l} \chi _\mathrm{lS} ) - \nabla \cdot (D_{lS} \rho _\mathrm{l} \nabla \chi _\mathrm{lS} ) = \bar{f}_\mathrm{S}. \end{aligned}$$

(18)

and that for water

$$\begin{aligned} \frac{\partial }{\partial t} (\phi _0 (S_\mathrm{l} \rho _\mathrm{l} \chi _\mathrm{lW} + S_\mathrm{h} \rho _\mathrm{h} \chi _\mathrm{hW})) + \nabla \cdot (q \rho _\mathrm{l} \chi _\mathrm{lW} ) - \nabla \cdot (D_{lW} \rho _\mathrm{l} \nabla \chi _\mathrm{lW} ) = \bar{f}_\mathrm{W}. \end{aligned}$$

(19)

The structure of these equations is similar to that of (15) except that the salinity component is not present in the hydrate phase. Based on (1a), we can divide (18) by $\rho _\mathrm{l}$; renaming $\chi _\mathrm{lS}=\chi _\mathrm{lS}$, we obtain the salinity part of (2).

1.2 Pressure Equation

To derive the pressure equation, we add (15), (18), (19). Collecting terms and taking into account volume constraints $S_\mathrm{l}+S_\mathrm{h}=1$ as well as $\chi _\mathrm{lM}+\chi _\mathrm{lS}+\chi _\mathrm{lW}=1$ and $\chi _\mathrm{kM}+\chi _\mathrm{hW}=1$, we see that the accumulation term becomes $\frac{\partial }{\partial t} \phi _0 (S_\mathrm{l} \rho _\mathrm{l} + S_\mathrm{h} \rho _\mathrm{h})$. The advection term becomes $\nabla \cdot (q\rho _\mathrm{l})$, and the diffusion term

$$\begin{aligned} R_D:= & {} - \nabla \cdot ( S_\mathrm{l} \rho _\mathrm{l} \phi _0 ( D_\mathrm{M}^0 \nabla \chi _\mathrm{lM} + D_\mathrm{W}^0 \nabla \chi _\mathrm{lW} + D_\mathrm{S}^0 \nabla \chi _\mathrm{lS} )\nonumber \\= & {} - \nabla \cdot ( S_\mathrm{l} \rho _\mathrm{l} \phi _0 ( D_\mathrm{M}^0 \nabla \chi _\mathrm{lM} + D_\mathrm{W}^0 \nabla (1-\chi _\mathrm{lM}-\chi _\mathrm{lS}) + D_\mathrm{S}^0 \nabla \chi _\mathrm{lS} )\nonumber \\= & {} - \nabla \cdot ( S_\mathrm{l} \rho _\mathrm{l} \phi _0 ( (D_\mathrm{M}^0 -D_\mathrm{W}^0)\nabla \chi _\mathrm{lM} + (D_\mathrm{S}^0-D_\mathrm{W}^0) \nabla \chi _\mathrm{lS} ) . \end{aligned}$$

(20)

Assuming that all diffusivities are equal, $\nabla \sum _\mathrm{C=M,W,S} \chi _\mathrm{lC}=\nabla (1)=0$ and the diffusion term $R_D$ vanishes. After further simplifications based on (1a), we obtain

$$\begin{aligned} \frac{\partial }{\partial t} \bigg (\phi _0 \big (S_\mathrm{l} (1-\frac{\rho _\mathrm{h}}{\rho _\mathrm{l}}) + \frac{\rho _\mathrm{h}}{\rho _\mathrm{l}}\big )\bigg ) + \nabla \cdot q = f, \end{aligned}$$

where $f = f_\mathrm{M}+f_\mathrm{S}+f_\mathrm{W}$. If furthermore $\phi _0$ is assumed constant in time, and $\rho _\mathrm{h} \approx \rho _\mathrm{l}$, then, after some algebra, we obtain the steady-state pressure equation

$$\begin{aligned} \nabla \cdot q =f. \end{aligned}$$

(21a)

The Eq. (21a) is coupled with Darcy’s law

$$\begin{aligned} q = - \frac{K}{\mu } (\nabla P_\mathrm{l} - \rho _\mathrm{l} G \nabla D(x)) \end{aligned}$$

(21b)

In the absence of sources and $f\equiv 0$, q is divergence free. In the 1D case, $q = const$ and is equal to the flux across the boundary $\partial \varOmega$. In fact, due to low fluxes q, the pressure is usually close to the hydrostatic $P_\mathrm{l}^0(x)$ pressure defined by (1b).

1.3 Conversion Factors

The conversion factor between $\chi _\mathrm{lS}$ and $\chi _\mathrm{lS}^\mathrm{m}$ is computed as follows. Assume we have 1L seawater with weight 1.03 kg=1030 g. Let the salinity be $\chi _\mathrm{lS} = 0.035\, {\mathrm {kg/kg}}$. In the volume of 1 L, this corresponds to $0.035 \times 1030=36.05 {\mathrm {g}}$. Using molecular mass of $58.44 {\mathrm {g/mol}}$ of NaCl, we see that this gives $36.05/58.44=0.62$ moles of NaCl. The same volume 1 L of seawater has $N=1030/18.02=57.2$ moles, since $18.02\,{\mathrm {g/mol}}$ is the molecular weight of water. Thus, the mole fraction corresponding to $\chi _\mathrm{lS}=0.035 {\mathrm {kg/kg}}$ and $S=35 {\mathrm {g/kg}}$ is therefore $\chi _\mathrm{lS}^\mathrm{m}=0.62/(0.62+57.2) =0.01$.

Further conversion factors are needed. In particular, for $f_\mathrm{M}$, we have

$$\begin{aligned}&1 \frac{\hbox {mmol}}{{\mathrm {m}}^3 {\mathrm {year}}}\nonumber \\&\quad = 10^{-3} \frac{{\mathrm {mol}}}{{\mathrm {m}}^3 {\mathrm {year}}} \times 16 \frac{g}{mol} \times 10^{-3} \frac{{\mathrm {kg}}(\hbox {CH}_4)}{\hbox {g}} \times \frac{1}{1030} \frac{\hbox {m}^3}{\hbox {kg} (\hbox {SW})} \times 1000 \frac{\hbox {year}}{\hbox {kyr}}\nonumber \\&\quad = 1.55 \times 10^{-5} \frac{\hbox {kg } (\hbox {CH}_4)}{\hbox { kg(SW)}\hbox { kyr}} \end{aligned}$$

(22)