Numerical Simulation of Hydrate Formation in the LArge-scale Reservoir Simulator (LARS)

: The LArge-scale Reservoir Simulator (LARS) has been previously developed to study hydrate dissociation in hydrate-bearing systems under in-situ conditions. In the present study, a numerical framework of equations of state describing hydrate formation at equilibrium conditions has been elaborated and integrated with a numerical ﬂow and transport simulator to investigate a multi-stage hydrate formation experiment undertaken in LARS. A veriﬁcation of the implemented modeling framework has been carried out by benchmarking it against another established numerical code. Three-dimensional (3D) model calibration has been performed based on laboratory data available from temperature sensors, ﬂuid sampling, and electrical resistivity tomography. The simulation results demonstrate that temperature proﬁles, spatial hydrate distribution, and bulk hydrate saturation are consistent with the observations. Furthermore, our numerical framework can be applied to calibrate geophysical measurements, optimize post-processing workﬂows for monitoring data, improve the design of hydrate formation experiments, and investigate the temporal evolution of sub-permafrost methane hydrate reservoirs.


Introduction
Gas hydrates are ice-like crystalline compounds made of lattices of hydrogen-bond water molecules, in which hydrate-forming gas molecules are embedded [1,2]. Various hydrate-forming gases have been identified thus far, including some typical smaller hydrocarbons (e.g., CH 4 , C 2 H 6 ) and inorganic compounds (e.g., H 2 S, CO 2 ), with methane (CH 4 ) being the most common one [3]. Naturally occurring gas hydrates are stable at elevated pressures and low temperatures, commonly present in marine environments and permafrost regions, where these conditions are fulfilled [1]. Generally, approximately 97% of natural gas hydrates (NGH) are reported to be concentrated in marine environments, whereas the rest are accumulated below permanently frozen strata [4], such as the Mallik NGH-bearing site in the Mackenzie Delta of Canada [5].
Gas hydrates are formed when gas and water molecules are in contact at high-pressure and low-temperature conditions [6]. Due to the distinct cage-like hydrate structures, a high amount of gas can be embedded in the three-dimensional network of water cages. According to conservative predictions, the total amount of CH 4 (ca. 500-2500 Gt) [7] stored in worldwide NGH reservoirs is estimated to exceed the proven CH 4 inventory of conventional gas reservoirs by about one order of magnitude [8]. Hence, NGH are considered as an alternative fossil energy source, and in-depth research is still required to study the NGH formation [9] and dissociation kinetics [10] as well as their efficient utilization by sustainable extraction technologies [11]. Although CH 4 extraction from NGH reservoirs has been investigated for almost four decades, it is still far from commercial production, and various knowledge gaps need to be addressed by scientific studies.
The "excess-gas" method, also termed as "gas injection" by Fitzgerald et al. [16] and originally presented by Handa and Stupin [12] in 1992, is the most widely used approach for hydrate formation in the laboratory. CH 4 hydrate forms in the space where injected water accumulations are trapped around grain contacts, exhibiting a cementation habit by employing the "excess-gas" method. Generally, the hydrate growth rate in a gas-rich environment generated by the "excess-gas" method proves to be orders of magnitude higher than the "excess-water" method achieves in water-dominated systems. In comparison, Priest et al. [14] reported that the outer layers of injected gaseous CH 4 bubbles are surrounded by CH 4 hydrate shells, showing a matrix-supporting habit by employing the "excess-water" method. CH 4 hydrate formation via the "excess-water" method requires a simpler laboratory setup and shorter experimental time periods in comparison to the "dissolved-gas" method. For mimicking the natural conditions present in hydrate-rich sediments, the "dissolved-gas" method was initially proposed by Spangenberg et al. [52], using a sample cell filled with glass beads outfitted for a micro-scale experimental setup. Here, it was demonstrated that hydrate saturation reached approximately 95% until the termination of the experiment by the decrease in permeability. Although the occurrence of pore-filling hydrate habit was reported for the "dissolved-gas" method [52,53], pore-filling hydrate naturally turns into matrix-supporting hydrate when the local hydrate saturation reaches 25-40% [36,54]. Moreover, Waite et al. [42] and Choi et al. [44] achieved a good balance between hydrate growth rate and high bulk hydrate saturation by using the "excess-gas" method to initiate hydrate nucleation before the continuation of hydrate growth by circulating dissolved CH 4 -rich fluid.
In addition, Stern et al. [55] suggested a special CH 4 hydrate formation technique called the "ice-seeding" method, which was originally designed to generate core-scale CH 4 hydrate samples from ice seed (pure H 2 O), for further mechanical testing and the investigation of hydrate dissociation patterns. Although this method is rarely used to form CH 4 hydrates, Spangenberg et al. [56] employed the "ice-seeding" method in combination with partial freezing to form CH 4 hydrates in sand samples.
According to our knowledge, most of these laboratory experiments involving CH 4 hydrate formation in synthetic sediments (Table 1) have not yet been reproduced by numerical simulations, except for the work of Yin et al. [3,45,46], who conducted several numerical investigations by means of the TOUGH+HYDRATE simulator [57]. The authors aimed to explore different kinetic and equilibrium CH 4 hydrate formation models to improve the description of the "excess-water" method and establish a sensitivity analysis on the hydrate saturation distribution concerning different multi-stage cooling schemes within a core-scale cylindrical sample chamber [3,17,45,46]. Although many hydrochemical models and numerical codes have been developed and implemented to study CH 4 hydrate production as summarized by White et al. [58], only a few of these are capable of reproducing hydrate formation by the "excess-gas" and "excess-water" methods. Among those numerical implementations, HydrateResSim (HRS) [59,60] is the only available open-source and open-access code, describing both equilibrium and kinetic models of hydrate formation, and only a few studies [61,62] have made use of it recently.
The LArge scale Reservoir Simulator (LARS) [15] has been established to study intermediate processes during hydrate formation via dissolved CH 4 at reservoir conditions and various hydrate dissociation strategies. Laboratory tests previously undertaken in LARS offer data for calibrating numerical models to further improve process understanding as well as experimental strategies and workflows. To the authors' knowledge, there is currently no numerical modeling study published that represents the observed hydrate formation or dissociation processes in LARS. Furthermore, CH 4 hydrate formation using the "dissolved-gas" method has not been simulated at laboratory scale. Consequently, the present study aims at developing a suitable numerical framework, verifying it against an established numerical simulator, and calibrating and validating it using a hydrate formation experiment undertaken in LARS.
For that purpose, a framework of equations of state (EOS) to simulate the physical properties of water with dissolved NaCl as well as CH 4 and equilibrium CH 4 hydrate formation has been developed, as demonstrated in Appendix A.2. The EOS was then implemented and integrated with the TRANsport Simulation Environment (TRANSE) [63] to investigate time-dependent and spatial CH 4 hydrate formation in a porous medium at pressure and temperature (p-T) conditions representative for the Mallik site [27,64]. The resulting simulation tool is referred to as T plus H (TRANSE+Hydrate) in this study. Our simulation results demonstrate that the numerical model implementation is capable of reproducing the main processes of hydrate formation in LARS, so that it can substantially support the further development of the experimental design and investigation of hydrate formation in water-dominated hydrate-rich sediments at the field scale.

Experimental Data
So far, eleven laboratory experiments have been successfully conducted in LARS, illustrated in Figure 1, including five different investigations into hydrate dissociation induced by thermal stimulation [15,39]. Additionally, three other experiments focus on CH 4 hydrate formation along with dissociation triggered by depressurization [26,27,40,41], while the rest of the tests focused on CH 4 -CO 2 or CH 4 -CO 2 -N 2 exchange processes [21,25]. Those previously undertaken and well-analyzed laboratory tests provide numerous resources and supplement materials for the calibration and validation of numerical models. However, the major processes of those studies in LARS have neither been reproduced by numerical simulations nor addressed the fundamental hydro-chemical characteristics of the mechanisms of CH 4 hydrate formation from supersaturated dissolved CH 4 in saline fluids.
Priegnitz et al. performed two experiments [27] to replicate the in-situ natural settings at the Mallik site [64] and form hydrate using the "dissolved-gas" method, before the depressurization-induced CH 4 hydrate dissociation was studied in LARS [27,40]. In this course, key parameters such as the time-dependent spatial temperature distribution and bulk pressure within the sediment sample were continuously recorded based on an automatic protocol executed during the hydrate formation processes. Additionally, other crucial variables, for instance bulk hydrate saturation (S h,bulk ) and fluid flow rate, were measured manually at regular intervals. Moreover, the workflow for quantification of the spatial hydrate saturation (S h ) distribution relied on achieving a thermal equilibrium before undertaking electrical resistivity tomography (ERT) [26,27] measurements to map the spatial S h distribution. Consequently, for the determination of S h from ERT, the CH 4 -loaded brine circulation was stopped to initiate temperature equilibration throughout the sandy sediment sample before performing the ERT measurement.
After a careful analysis of the experimental datasets, the early stages (ca. 15 days) of one experiment conducted by Priegnitz et al. [27] have been selected to serve as benchmark for model calibration and validation in the present study, as shown in Figure 2. Particularly, our research focuses on the conformance of simulated hydrate saturations with those derived from fluid sampling, temperature distribution, and ERT data collected during the hydrate formation experiment.  Table 2 for their coordinates) during the hydrate formation experiment in LARS (left) with their relative location (right, not to scale). Observations at T0-T12 are modified from a hydrate formation experiment conducted by Priegnitz et al. [27]. The accuracy of the RTD locations is in a range from 0.01 to 0.05 m, and the measurement error of the applied Pt100 RTD is ±(0.3 + 0.005T)°C.

Hydrate Formation Experimental Schedule
At the onset of any hydrate formation experiments in LARS, a plastic mesh plate with fourteen mounted Pt100 temperature sensors (resistance temperature detector (RTD)) was first installed into the cylindrical sample chamber that was isolated by a neoprene rubber jacket. Subsequently, the sample chamber was filled with quartz sand and sealed by the lid of the pressure vessel from its top. Thereafter, the vessel lid combined with the sample chamber was inserted into a cylindrical autoclave, where the sample is separated from the cooling liquid, circulating in between the neoprene rubber jacket of the sample chamber and autoclave wall. Finally, the nuts and bolts of the autoclave were secured to complete the installation of the sample in LARS, and hydraulic integrity along with the availability of all types of sensors installed were verified before initiating the experiment.
Before the start of the experiment, a CH 4 -free saline solution (3.68 g NaCl·L −1 ) [27] sourced from the pore fluid container was circulated through a stainless porous filter plate from the bottom of the sample chamber to drive the air out. After the saturation procedure, a confining pressure (ca. 14-15 MPa) was applied to the sandy specimen by pressurizing the coolants circulating through the cooling chamber. The amount of water injected in the saturation procedure and expelled during the build-up of confining pressure was determined as a prerequisite for the estimation of the intrinsic porosity and S h,bulk . In the next step, the brine was loaded with CH 4 and pressurized (ca. 11 MPa) in the dissolved gas charging vessel ( Figure 1). It was then pumped into the sample chamber from the top through the hydrate-bearing sand at constant pressure and temperature. The prescribed temperature was slightly above hydrate equilibrium temperature at the given pressure (ca. 13.6°C at 11 MPa) to avoid the porous filter plate at the inlet from clogging by forming hydrates.
Inside the sample chamber, the temperature of the flowing pore fluid was reduced by the cooling circulation system to approximately 3.5°C to achieve hydrate stability conditions. A considerable temperature drop of the inflowing warm CH 4 -rich brine occurs close to the neoprene rubber jacket, and additional hydrate is thus formed when CH 4 solubility is decreased in the presence of hydrate nucleation [53]. Consequently, the CH 4 -supersaturated pore fluid continuously releases CH 4 to form hydrate until the CH 4 concentration is reduced to maximum CH 4 solubility. The outflowing brine is then heated and reloaded with additional CH 4 in the gas charging vessel before re-entering the sample chamber for the next flow-through cycle.
Three major pore fluid circulation stages marked by Roman numbers (I to III) were considered in the numerical simulations, whereby Stage I has been divided into two Substages (I-1 and I-2) due to an unintentional interruption of the warm inflowing fluid flow for around 9 h, as indicated in Substage I-1-1 ( Figure 2). Excluding this interruption, other intentional interruptions (Substages I-1-2, II-2, and III-2) of the inflow of warm CH 4 -charged water resulted in a temperature "equilibration" of the sandy sediment sample close to the temperature of the circulating confining pressure fluid and a decline of CH 4 available for hydrate formation.
The ERT measurements taken at the end of Substages I-1-2, II-2, and III-2 ( Figure 2) produced the most reliable resistivity distributions because the stationary conditions reduced the fluctuations in temperature and also improved the quality of the collected data. As CH 4 hydrate is an electrical insulator, ERT measurements allow for the determination of the spatial S h distribution in the sample chamber [27]. Hereby, ions from the dissolved salt accumulate in the pore fluid, as only CH 4 and pure water are consumed during hydrate formation. As a result, the electrical conductivity of the pore fluid increases. Furthermore, the mass of accumulated CH 4 hydrate can be determined by using the S h -dependent electrical conductivity approach presented by Waite et al. [42] and Spangenberg et al. [52]. Based on that approach, the spatial S h can be derived from spatial variation of the electrical conductivity in the hydrate-bearing sand.
According to Waite and Spangenberg [53], the amount of CH 4 available for hydrate formation at about 5°C amounts to approximately 42% of the initial CH 4 solubility in brine at 20°C. By circulating the warmer CH 4 -rich pore fluid through the sandy sediment sample, the accumulation rate of S h,bulk increases by 2 to 4.5% per day, filling almost 31% of the sample's pore space after Stage III.
Out of the fourteen installed RTDs, twelve operated and two malfunctioned (T10 and T13), with the latter excluded from previous studies [27]. The calibration of RTDs was conducted before their installation inside LARS during the preparation of the first hydrate formation test in 2011. In addition, the original measurement deviations of the RTDs T4 and T8 were both 4.2°C, falling out of the average measurement deviation of the other RTDs (3.3°C). According to the correlation between the distance of these RTDs to the fluid inlet and their temperature correction listed in Table 2, temperatures for RTDs T4 and T8 have been revised. Therefore, a re-correction was made in the present study using the proven temperatures of 3.2°C and 3.3°C for the originally calibrated temperatures for RTDs T4 and T8, respectively.

Mathematical Model
An equation of state (EOS) module for hydrate formation has been developed and coupled with a flow and transport simulator [63] in the scope of the present study to investigate the coupled hydro-thermo-chemical processes in LARS as discussed in Section 2.1.1 [26,27]. Kowalsky and Moridis [65] demonstrated that an equilibrium reaction model is a feasible alternative to a kinetic approach for simulating gas hydrate behavior at the reservoir scale.
However, temperature measurements in the LARS experiments were made every few seconds by the RTDs, and the sample volume of LARS is approximately 210 L. Therefore, we were not able to state until now that hydrate formation processes can be described by an equilibrium reaction approach given the aforementioned conditions. For the representation of short-term and core-scale (typically around 0.1 to 10 L) hydrate formation processes, the kinetic model is accurate and able to capture the transitional results of intermediate states [60,65], but its requirements in terms of computational power and numerical model convergence are substantially higher. Therefore, one objective of the present study was to investigate whether an equilibrium reaction approach is capable of representing hydrate formation via dissolved CH 4 in LARS using multi-stage cooling.

Modeling Assumptions
The developed equilibrium model utilizes the temperature and pressure-dependent relation proposed by Moridis [66] at the hydrate-aqueous equilibrium. According to Kashchiev and Firoozabadi [67], the aqueous solution has to be supersaturated with the hydrate-forming gas at the given pressure and temperature conditions; hydrate crystallization can then occur as the supersaturated gas is encased by the hydrate structure. Consequently, CH 4 hydrate formation or precipitation can be defined by the following reaction [2]: where the hydration number, n, commonly equals 5.9 in series of experimental studies undertaken in LARS [27,39,41], with CH 4 hydrates of cubic structure I (S I ) [6] being formed. The applied numerical framework allows for conducting quantitative descriptions of the involved coupled thermal, hydraulic, and chemical processes in hydrate-bearing sand, which are presented in Appendix A. Hereby, fluid migration is governed by density-driven flow in porous media (Darcy's Law), considering advective and diffusive transport of dissolved CH 4 and NaCl in the pore fluid. Moreover, heat transport and thermal energy exchange occur via conduction and convection [63], complemented by the equilibriumbased CH 4 hydrate formation reaction.
In order to maintain the accuracy of the numerical solution of the non-linear system of partial differential equations, the underlying simplifications were considered to maintain computational efficiency and numerical convergence requirements: The porous medium is completely filled by pore fluid and/or CH 4 hydrate, with single-phase flow considered in the entire modeling domain; 2.
Deformation of the porous medium (hydrate-bearing sand) is assumed to be negligible due to the applied confining pressure of 14 to 15 MPa, with the porous medium matrix being evenly compacted and homogeneous; 3.
Thermophysical properties of the aqueous solution do not consider the effects of the dissolved CH 4 , as these are negligible for the present study. The dissolved inhibitor (NaCl) influences neither the molecular structure of the formed CH 4 hydrate nor the rate of hydrate formation, but fluid density, viscosity, heat conductivity and capacity as well as CH 4 solubility, only; 4. CH 4 from the supersaturated aqueous phase is directly consumed by equilibrium hydrate formation without any intermediate phase changes and side reactions; 5.
Mobile components contain the aqueous phase with dissolved CH 4 and NaCl. All water-soluble species and liquids are non-volatile at the applied temperature range (0-25°C) and pressure conditions (ca. 11 MPa).
The simplification of the inhibition effect of NaCl and other salts on hydrate formation is attributed to the fact that the salt ions bind water molecules, which are then no longer available for dissolving CH 4 molecules in the aqueous phase. Thus, the amount of CH 4 available for hydrate formation is reduced in saline aqueous solutions compared to deionized water, as presented by Malagar et al. [68] and literature cited within it. In LARS, the saline fluid is almost fully saturated with CH 4 when it leaves the dissolved CH 4 charging vessel ( Figure 1) and enters the sample chamber. Due to the continuous flow, CH 4 is continuously supplied as a hydrate former, so that the salting-out effect described above is reduced or even completely eliminated. Therefore, an instantaneous formation of CH 4 hydrate within LARS is observed when the salinity-corrected p-T condition is met under the equilibrium CH 4 hydrate formation approach.

Numerical Model Verification
The objective of the benchmark study discussed in the following was to verify the coupling between the CH 4 hydrate formation EOS implemented in the present study with the fluid flow and transport simulator presented by Kempka [63]. For that purpose, the well established numerical simulator HydrateResSim [59,60] (HRS) has been used as reference. Figure 3a shows the 1D modeling domain, where the first left element acts as a cooling boundary at a constant temperature of 4°C under the assumption of the presence of a negligible amount of hydrate nucleation. The pore space of all other elements is filled with CH 4 -saturated water at an initial temperature of 16°C. At the impermeable cooling boundary, heat exchange is allowed between it and its neighboring element. Figure   The supersaturated dissolved CH 4 is instantly consumed by hydrate formation as the water temperature drops from 16°C (blue dot) to 4°C (green dot) during 50 days of simulation. Additionally, Figure 3c shows the temperature distribution in the model at simulation times of 0.1, 1, 3, 5, 10, and 50 days. Figure 3d presents the temporal evolution of the dissolved CH 4 concentration (C CH4 ) and CH 4 hydrate saturation (S h ) at the right boundary for 50 days of simulation time as computed by our model (T plus H) and HRS.
Overall, the maximum relative deviation between the main simulation results, i.e., temperature, S h and C CH4 , calculated by T plus H and HRS is < 0.5%. The main reasons for the deviations are attributed to the application of different equations of state, as well as the distinct realization for the same initial and boundary conditions in both simulators (i.e., the cooling boundary is implemented as an element with infinite volume in HRS, whereas it is a finite volume element in T plus H). Further, additional error sources for these deviations may be attributed to the application of different temporal and spatial discretization schemes. Following the results of this benchmark, the T plus H simulation results show a similar high accuracy, so we conclude that our code is capable of addressing the main objective of the study: the simulation of the hydrate formation experiment undertaken in LARS as discussed in the following section.

Model Implementation to Reproduce the LARS Experiments
Following the successful model verification, T plus H is calibrated using the experimental LARS data. For that purpose, temperature profiles determined by the installed RTDs and spatial hydrate saturations derived from geophysical monitoring as well as fluid sampling were reproduced numerically.

Initial and Boundary Conditions
The initially homogeneous thermophysical properties of the porous medium in LARS change with the increase in S h , i.e., effective porosity and permeability as well as effective heat conductivity of the immobile components decrease. The sampled S h,bulk exceeded 89% and the local S h observed by ERT reached ca. 94.2% at the end of the hydrate formation experiment, while the minimum local effective permeability calculated by the Carman-Kozeny relation (Equation (A13)) was 28.8 mDarcy [27]. In contrast to the ERT observations, the measured effective permeability was 2 mDarcy [31] at the final stage, maintained by a local S h of 97.5% and residual pore fluid saturation of 2.5%. Applicable thermophysical properties for the porous medium in LARS were determined using an iterative matching approach and are summarized in Table 3. Only limited information is provided in the description of the hydrate formation experiment in [27], comprising the estimated temperature ranges of the fluid at the inlet and the surrounding coolants as well as the average S h,bulk accumulation rates (ca. 2% per day). Other data required for model parametrization, such as the average inlet fluid rate (ca. 80 L per day), were derived from experimental records of an identical experiment. As the actual fluid temperature after passing through the porous filter plate was not determined in the laboratory study, it was estimated from the provided temperature thresholds (13.6-16°C). Hereby, the lower limit was chosen to ensure that the inflowing brine temperature remains above hydrate stability conditions (ca. 13.6°C) at the given pressures. The evaluated upper temperature limit (ca. 16°C) was determined based on a measurement in the dissolved CH 4 charging vessel undertaken at the start of the experiment.
Moreover, the pore fluid control system (Figure 1) was implemented by means of a Neumann boundary condition in the numerical model (Figure 4a). The coolant circulation system and the confining chamber were introduced as Dirichlet boundary conditions with impermeable hydraulic properties, fixed pore pressure, and constant temperature gradient linearly increasing from the fluid inlet (ca. 3.5°C) to the outlet (ca. 4.0°C). The coolant temperature was measured once at the cooling chamber inlet, and then the coolant was assumed to be heated gradually by the thermal energy transmitted through the neoprene jacket from the sediment chamber. The temperature of the recycled coolant was determined once before entering the heat exchanger ( Figure 1) at less than 4.0°C. The main variables used to determine thermo-physically reasonable initial and boundary conditions, implemented for model calibration by means of an iterative history-matching procedure, are summarized in Table 4.

Results and Discussion
One multistage CH 4 hydrate formation experiment was chosen for the following model validation as indicated in Figure 2. It is also referenced as "LARS RUN2" [27] in the research on spatial S h characterization by ERT as well as "experiment A" [40] in the investigation of gas production triggered by a multistage depressurization scheme. The temporal evolution of temperature profiles recorded at twelve RTDs is compared with the simulation results for model calibration in Sections 3.1.1 and 3.2.1. Hereby, the location of the temperature sensors was accordingly calibrated, and thus revised as discussed in Section 3.1.2. After calibration, the model was validated by comparison of the temporal and spatial evolution of the simulated and observed hydrate saturations presented in Section 3.2.2. An overview of the experimental temperature evolution and numerical predictions at an early period of the hydrate formation experiment is shown in Figure 5.  The experiment and simulation started with the circulation of the CH 4 -saturated brine sourced from the gas-charging vessel (Figure 1), defined as hour zero of the experimental time in Figure 5. Before hour zero, the hydrate-bearing sand is assumed to be filled with a negligible amount of hydrate crystals, and the porous filter plate is partially clogged by a mixture of hydrate and sand. Stage I is regarded as the most representative period of the hydrate formation experiment, considering that each period is influenced by different effects discussed in the following.
As outlined in Figure 2, Stage I has been divided into two Substages (I-1 and I-2) due to the occurrence of an unexpected discontinuity in the provision of the warm CH 4 -rich fluid during the time period 33.4 to 47.5 h (Table 6). After the fluid flow interruption, the warm CH 4 -rich fluid re-entered LARS from hour 47.5 on, until the pumping system was shut down 12.5 h later for ERT measurements to be taken at hour 90. Considering the uncertainties related to the manual temperature and rate control of the injected fluid, data in Figures 2 and 5 suggest that the inlet fluid parameters were occasionally not fully maintained according to the experimental plan.
The earliest rapid temperature increase was captured at RTD T0, with the shortest distance to the fluid inlet, and the peak temperature (slightly below 14.5°C) was recorded by RTD T1 during Stage I-1 (Figure 5a). The simulated warm and CH 4 -rich fluid reached T1, T2, and T3 at almost the same time, whereas the measured temperature front arrival delay between RTD T0 and T1 was about 2 h. In contrast to the respective experimental results, the numerically predicted arrival time of the elevated temperature front between T0 and T1 is over 0.8 h, and those between T1 and T2 as well as T3 are approximately 1.2 h, as the distances of T2 and T3 to the fluid inlet are identical.
RTD T0 was expected to record the highest temperatures during the entire experiment duration, because the inflowing warm fluid should reach RTD T0 first under the assumption of a homogeneous porous medium. However, the temperature curve at T0 barely reached 12.5°C at the beginning, and then it gradually declined to about 10.8°C after 4 h. The highest temperature reading at RTD T1 was more than 3°C higher compared to the corresponding simulated one at T0 in Stage I-1. This anomalous behavior did not occur during the remaining experiment, where the numerical predictions at the T1 location were in line with the observed temperatures. Despite the possibility of instrumental failures and spatial displacement of RTDs during their installation, it is reasonable to assume that the region near RTD T0 had a higher CH 4 hydrate saturation than that obtained from the simulations. The region near RTD T1 was assumed to have a lower CH 4 hydrate saturation than the simulated one, and thus more inflowing warm water was redistributed from the area near sensor T0 to the location of T1. This may be explained by a considerable amount of hydrate being present at the top of the sample chamber before hour zero of the experiment, resulting in RTD T0 being coated or surrounded by hydrate much earlier than in our simulation.
The simulated temperatures obtained for the RTD positions T2-T6, T8, T11, and T12 match very well with their corresponding observations ( Figure 5). Additionally, the general tendency of the simulated temperatures at RTDs T3, T7, and T9 is consistent with the observations, even though it shows maximum deviations of 8% during the time period of 10 to 33 h.

Calibration of RTD Locations
RTD locations in the model were adjusted to calibrate the simulated temperatures by the observed ones. For the simulation results presented in Figure 5, the obtained numerical predictions were not extracted from the exact coordinates plotted in Table 2, because (1) the RTDs' actual spatial detection range as well as pressure and flow rate sensitivity regarding its measurement accuracy are unknown; (2) it is very likely that unquantifiable deformation has been introduced during the installation of the sample chamber when it was hoisted for mounting into the pressure vessel; (3) further immeasurable deformation may occur during compaction of the sediment sample when the confining pressure is initially applied; (4) inevitable position deviation may emerge during the manual installation of the temperature sensors onto the reserved holes of the plastic frame.
To improve the match between the simulated and experimental data, the spatial RTD positions have been adopted within a range of 0.01 m (region close to sample chamber top) to 0.04 m (region close to sample chamber bottom), except for RTD T12 (Table 5). This noticeable deviation of the revised location of RTD T12 may be attributed to the heterogeneity introduced by local compactions of the hydrate-bearing sand during the installation of the counter-current heat-exchange reactor [15,39] (Figure 5b) in the experimental setup. In summary, data in Figure 5 show that the majority of the simulated temperatures are in very good agreement with the observed ones, excluding a few RTDs positioned close to the sample chamber boundaries. Consequently, the applied numerical model has been successfully calibrated using the temperatures recorded at the RTDs.

Model Calibration and Validation
By implementing the previously introduced initial and boundary conditions within the bandwidths listed in Table 4, the best-fit combinations with minimum deviations from the observations were obtained from the model calibration based on the Stage I results in the previous section. Subsequently, further model calibration based on Stages II and III was achieved via an iterative history-matching procedure. As a result, the model parametrization was revised accordingly ( Table 6). As shown in Figure 2 and Table 6, the CH 4 -loaded brine inflow periods are represented by Substages I-1-1, I-2-1, II-1, and III-1, with the rest of the experimental duration determined by brine inflow suspension periods. Table 6. Main initial and boundary conditions employed in the simulation study on the CH 4 hydrate formation experiment performed in LARS (cf. Figure 2 for the division of Substages).

Model Calibration via Comparison of the Temporal Evolution of Simulated and Observed Temperature Profiles
During the LARS experiment and numerical simulation, a contribution to the temperature increment within the sample chamber is made by the combined effect of inflowing warmer CH 4 -loaded fluid and the latent heat released from hydrate formation. Generally, the heat release of hydrate formation is ca. 54.4 kJ (mol CH 4 ) −1 [71]. As reported by Waite and Spangenberg [53], given a small dissolved CH 4 consumption rate, the temperature increment is limited to < 0.5°C even under the assumption that the components present in each representative elementary volume instantaneously absorb all released heat. However, the forming hydrate would gradually fill the available pore space in at least 200 h so that the contribution to the temperature increment from the inflowing warmer CH 4 -rich fluid is an order of magnitude higher than that from hydrate formation.
From 50 h on (Figure 6a), the measured and simulated peak temperatures are always observed at RTD T2. This shows that the temperature increment contributed by warm fluid in the vicinity of RTDs T0, T1, and T3 is less than that in Stage I-1 due to the local permeability decrease induced by hydrate accumulation. In contrast, warmer fluid flowed along T2 via the center of the sample, where permeability was much higher than at the nearby model boundary. This phenomenon indicates that a high permeable flow channel existed along the central axis of the model geometry near T2. Consequently, CH 4 hydrate is primarily formed at the model boundary close to T0, T1, and T3, where temperature is determined by external cooling in contrast to the vicinity of T2.  Figure 5b and Table 5). LARS T0-T12 data were adapted from the CH 4 hydrate formation experiment conducted by Priegnitz et al. [27].
Moreover, the vertical distance sequence of the RTDs in the upper sediment was T0 < T1 < T3 < T4 (sorted from top to bottom), whereas the observed temperature sequence was T1 > T0 > T3 > T4 (sorted from high to low temperatures) in Stages I-1, I-2, and II-1. The observed temperature sequence then changed to T1 > T3 > T0 > T4 in the early period of Stage III-1 and ended with T3 > T0 > T1 > T4 afterwards (Figure 6a). This phenomenon confirms the previous explanation that the warm inflowing fluid was redirected to RTD T1 at the start of the experiment due to pre-existing hydrate in the region of RTD T0. Moreover, hydrate was constantly amassed around RTD T1 until almost full hydrate saturation was achieved at this location before Stage III. The redistributed warmer and CH 4 -rich inflowing fluid then advanced to the region around RTD T3, and thus the observed temperature at RTD T1 became lower than that at RTD T0, due to its larger distance to the fluid inlet at the top.
However, the numerically predicted temperature sequence for the revised RTD positions maintains T0 > T1 > T3 > T4 (sorted from high to low temperatures) until the end of Stage II-1, complying with the aforementioned distance relation. In Stage III-1, the temperature sequence changes to T0 > T1 > T4 > T3, showing that the warmer inflowing fluid was redirected to T4 when the nearby region of T3 was occupied by hydrate. Hereby, the nearby regions of T0 and T1 were not saturated with hydrate at the end of Stage III, as illustrated in Figure 7a. These findings further support the hypothesis that hydrate must have been present near RTDs T0 and T1 at the start of the experiment. In Stage II-1, similar temperature changes are observed at T0, T1, and T3, whose simulated temperatures gradually drop by 0.6 to 1.0°C. However, a reverse temperature change is observed during the same period at T4, whose simulated temperature steadily increases by approximately 1.1°C (Figure 6a). This phenomenon indicates that the inflowing warm CH 4 -loaded fluid was redirected from the nearby regions of lower permeability (T0, T1, and T3) to those around T4 with higher permeability. In addition, the simulated temperature development at T4 reflects that the hydrate formation process successively generated heat in the region around T4. This causes the simulated temperature at T4 to deviate from the corresponding observations by up to 13% during the experimental time period of 140 to 160 h.
During Stage III, the most noticeable difference between the simulated and observed temperatures is found at T3 with a maximum deviation of 40% during the time period from 220 to 310 h (Figure 6a). In this period, the region around T3 exhibits a decreasing permeability with the continuous accumulation of hydrate, as illustrated in Figure 7c,d. Considering these deviations at the RTD at the sample top near the neoprene jacket, the influence of a buffer layer [34,[72][73][74] is expected. Although digital rock modeling [72] of the hydrate-bearing sample [34,73] is beyond the scope of the present study, the buffer layer of the unconsolidated sample built by gravity-driven sedimentation [74] is relevant to this study.
The latter indicates that the hydro-physical properties of the interface between the hydrate-bearing sand and the neoprene jacket as well as surrounding metal structures (the porous filter plates and the counter-current heat-exchange reactor) were probably influenced by a buffer layer. As a result, remarkably high porosities and one order of magnitude higher local permeabilities were observed at these locations [74]. Consequently, a certain amount of warmer inflowing fluid migrated downwards along this high-permeable outer layer to reach RTD T3, supplying extra heat and increasing the observed temperatures at RTD T3 beyond the corresponding numerical predictions.
A substantial difference between the temperature evolution at the RTDs in the upper sample (T0-T4) and those in the lower one (T5-T9, T11, and T12) becomes obvious from Figure 6. The simulated and observed temperatures in the upper sample are higher than those in the lower one by more than 2°C before Stage III. This shows that the temperature of the downward-flowing warm CH 4 -loaded fluid is substantially reduced after passing RTDs T0-T4, whereby some amount of the initially dissolved CH 4 is consumed by hydrate formation, accompanied by the release of latent heat. Consequently, the highest hydrate accumulations occur in the upper sample according to the simulations and ERT measurements plotted in Figure 7a and b, respectively.
Although the distance from the fluid inlet to RTD T12 is farther than that to RTD T11, the lowest monitored and simulated temperatures are always observed at T11 rather than at T12. For the boundary region around T11, the heat transmitted from the inflowing fluid is superimposed by the external heat sink. In comparison, T12 is located in the aforementioned warm fluid flow channel, where the thermal conditions are contrary to those observed in the region near T11. Moreover, similar simulated and observed temperature change characteristics did not only exist at RTDs T8 and T9 within the warm fluid flow path, but also at RTDs T5 and T6 near the buffer layer. The same temperature change characteristics are also present for the observed and simulated results at RTDs T7 and T12. Despite their maximum deviations extending to up to 9% for the time period of 140 to 160 h, simulated temperatures in the lower sample (T5-T9, T11, and T12) matched well with the corresponding observations after Stage I (cf. Figure 6).

Model Validation by Comparison of the Temporal and Spatial Evolution of Simulated and Observed Hydrate Saturations
The evolution of spatial hydrate saturation and hydraulic permeability distributions can be described by the numerical simulations and ERT observations, separately. At hour zero of the experiment shown in Figures 7b and 8b, S h,bulk > 7.2% was determined by ERT [27]. At the end of Stage I (90 h), it is indicated by the simulation results that the majority of the formed hydrate was accumulated under the bottom of the top porous filter plate and distributed along the warm-cool fluid interface present at 1.3-0.65 m of the specimen height. Moreover, the simulated S h,bulk reaches 9.0% (Figure 7a), showing good agreement with the pore fluid sampling data (10%), disregarding the deviation to S h,bulk of 20.4% determined from ERT measurements (90 h in Figure 7b).

Figure 8.
Observed and simulated bulk permeability alongside with the bulk hydrate saturation (S h,bulk ) evolution during the hydrate formation experiment: (a) comparison of the simulated volumeaveraged permeability with the numerically inversed ERT results [27], the Carman-Kozeny relation as a function of S h,bulk , and the bulk permeability determined by hydraulic testing (blue dot) [27] at the end of experiment [27]; (b) comparison of the simulated S h,bulk against the sampled and ERT-measured ones, modified from Priegnitz et al. [27].
In Stage II, hydrate formation mainly occurred at 1.3-0.95 m of the specimen height, as shown for 190 h in Figure 7a,b. Furthermore, the maximum simulated S h increased to 97% and the simulated S h,bulk to approximately 18.9% simultaneously, which closely matches with the fluid sampling result (21%) at 190 h. Subsequently, the front of accumulated hydrate advanced to 0.95-0.65 m of specimen height, indicating a similar hydrate distribution pattern to the ERT-measured findings (cf. 360 h in Figure 7a,b). Finally, the simulated S h,bulk increased to about 30.1%, and is thus almost identical with the corresponding S h,bulk of 31% determined by fluid sampling at the end of Stage III. Despite the relative error of the ERT-measured S h,bulk , it has been confirmed that hydrate initially formed at the top of the specimen, and then the front of hydrate formation advanced to the specimen center along the neoprene jacket.
In Figure 7a,b, the S h observed by ERT and simulated are virtually zero at the bottom and top of LARS. At the top region near the inlet of LARS, the p-T conditions of the inflowing CH 4 -loaded fluid are outside of the CH 4 hydrate stability range. These conditions were chosen to avoid undesired blockages by hydrate formation at the inlet. As the fluid flows through the sandy sediment sample, it cools and eventually meets the CH 4 hydrate stability conditions, whereby excess CH 4 is consumed. When the fluid flows to the bottom of the sandy sample, it is undersaturated with CH 4 , so that CH 4 hydrate can no longer form despite the favoring p-T conditions.
The Carman-Kozeny Equation (A13) shows that the bulk permeability is a function of S h,bulk , as the predefined initial permeability is constant (Table 3). Accordingly, the change in local permeability reflects the variation of S h . At hour zero of the experiment (Figure 7d), a slightly heterogeneously distributed permeability was observed in a range of 300 to 500 Darcy. Subsequently, accumulating hydrate reduced the hydraulic permeability in the upper part of specimen. Some warmer CH 4 -loaded brine advanced along the neoprene jacket, and a small amount migrated along the aforementioned buffer layer. Moreover, the electrical resistivity of the relevant region was appreciably increased [41]. On account of the constantly increasing hydrate accumulation, the permeability of the affected region was thus declining until 360 h. Finally, the minimum simulated local permeability is 2.1 mDarcy (Figure 7c) and identical to the hydraulic test results (2 mDarcy) represented by the blue dot in Figure 8a.
Overall, the simulation results confirm that ERT monitored the spatial hydrate distribution within the specimen mainly qualitatively in the early stages (before 190 h in Figure 7b) of the experiment, and the volume-averaged S h,bulk was estimated at high precision by pore fluid sampling in previous laboratory studies [27].

Uncertainties of Critical Parameters in the Experimental Study
Initially, several essential details were collected to determine practical ranges of the input parameters (Table 4) and establish an equivalent geometrical model for calibration (Figure 4). The optimized input parameters were then determined by a history-matching procedure during model calibration, and the documented initial and boundary conditions were revised (Tables 5 and 6). Additionally, the exact locations of RTDs could not be derived by any measurement method once the sample chamber was filled with quartz sand. Thus, these had to be revised iteratively in the present study (Table 5).
There are further uncertainties similar to those noted above. For instance, the initial intrinsic permeability of the specimen was too high and out of the measurement range of the experimental setup, thus it could exclusively be estimated empirically from porosity and grain size distribution. However, the final bulk permeability could be measured at the end of the hydrate formation experiment by a hydraulic test [27], whose result is illustrated in Figure 8a. In addition to this, the boundary condition for the inflowing fluid became timedependent, caused by the hydrate likely forming inside the porous filter plate. Furthermore, some critical variables required for the simulation study were exclusively measured at the onset of the hydrate formation study, e.g., the external cooling and inflowing fluid temperatures as well as the dissolved CH 4 concentration in the inflowing fluid. Despite these uncertainties, most of them can be prevented or eliminated in future experimental designs by optimizing the execution of the experimental planning.
With regard to the uncertainties in the estimation of the initial permeability and S h,bulk , the initial interpretation of Figure 8 was given in a previous study [27]. However, the remarkable differences between hydrate saturations derived from ERT measurements, conductivity measurements, and simulation are still worth noting and require further analysis. Figure 8b demonstrates that the simulation curve increases in the form of multiple steps, which is inconsistent with the simple linear growth pattern determined by pore fluid sampling and ERT measurements. It can be explained by the artificial acquisition intervals of experimental data being much larger than those in the numerical simulations. For example, the pore fluid sampling for pore fluid conductivity measurements was undertaken on a daily basis, whereas ERT measurements were only made at the end of each formation stage. Hence, these cannot characterize the detailed transition during hydrate accumulation.
Theoretically, the remarkable contrasts in electrical resistivities between the coexisting pore-filling components (pore fluid and hydrate) granted hydrate to be distinguished from the pore fluid and localized even at low S h [27], as hydrate is electrically non-conductive. However, it should be noted that the most significant deviation of the ERT-measured S h,bulk (ca. 7.2%) from the corresponding simulated and sampled S h,bulk (0%) was observed at the onset of the experiment (hour zero) in Figure 8b. This is attributed to the initial assumption and relevant post-processing routine for the conversion of electrical resistivity into S h [27]. Hereafter, the difference in hydrate content from ERT and pore fluid electrical resistivity becomes less notable as the hydrate saturation increases, particularly at S h,bulk > 65%, when the sampled S h,bulk finally agrees with the ERT-observed one. Eventually, the ERT-observed S h,bulk resulted in a notable deviation of more than 10% from the corresponding sampled S h,bulk (ca. 89%) in Figure 8b.
The ERT measurements provided useful information regarding the location where hydrate started to form and how its distribution generally changed with time, which helped to adjust the way the experiment was conducted. However, the upper and lower 0.15 m of the sample chamber are not covered with electrodes, and ERT is not able to capture effects close to the top and bottom steel closures with the fluid in-and outlets. In addition, the resolution of the ERT method is limited, and small-scale accumulations of hydrate with only little effect on electrical resistivity might not be recognized. Furthermore, the inversion process converting the resistance measurements on the sample surface into a 3D electrical resistivity distribution is not unique (as for all potential methods). The transformation of electrical resistivity into S h relies on an empirical rock-physical model, without any specific calibration for hydrates [26,27]. This further explains the observed discrepancies to the numerical modeling.

Summary and Conclusions
The present study numerically reproduced a previously conducted multi-stage CH 4 hydrate formation experiment. The effectiveness and accuracy of the developed coupled numerical framework have been evaluated and demonstrated by a benchmark and comparisons to the experimental observations. Consequently, the key parameters and an optimum combination of initial and boundary conditions were determined. Our findings allow for the following conclusions: 1.
The general consistency of the experimental observations with the simulation results proves that the employed equilibrium CH 4 hydrate formation model can represent the main processes of hydrate formation in LARS. The equilibrium reaction model is a practicable alternative to kinetic approaches at macro-scale (vessel volume > 0.2 m 3 ) given the application of the "dissolved-gas" method. In contrast, kinetic reaction approaches tend to be irreplaceable for modeling hydrate formation by other methods, because their CH 4 hydrate growth rates are orders of magnitude faster than that of the "dissolved-gas" method.

2.
The deviations among the experimental observations (i.e., continuously recorded temperature profiles, periodically gathered S h,bulk , and ERT-tomography derived spatial S h distributions) and the corresponding numerical predictions were minimized through an iterative optimization procedure. It has been indicated that the combination of the thermal properties of inflowing CH 4 -loaded fluid and the hydrate-bearing sand determine the spatial distribution of hydrate accumulations.

3.
The presented spatial S h distribution illustrates a heterogeneous accumulation within the hydrate-bearing sand at an early experimental period when S h,bulk < 30%, with the feature becoming less prominent until S h,bulk > 80%. 4.
In the LARS hydrate formation experiment, a relatively large temperature gradient (ca. 10°C/0.23 m) is generated between the inflow of warm brine and its surrounding coolants, leading to a heterogeneous hydrate distribution. In contrast, the subpermafrost and sub-seafloor geothermal gradients in natural settings are substantially lower (3°C/100 m) [75] and steady for long time periods, causing a lower and almost constant dissolved CH 4 concentration gradient in the saline fluid. Therefore, relatively uniformly distributed S h were found within the NGH intervals with ignorable lateral variations at the Mallik site. These NGH accumulation intervals could be simplified as CH 4 hydrate layers formed via the continuous supply of dissolved CH 4 , migrating through the up-dip natural faults in the Canadian Beaufort-Mackenzie Basin region.
The proposed numerical framework can be utilized to improve experimental designs and optimize post-processing workflows of monitoring data. Thereby it could contribute to calibrate the advanced geophysical identification techniques and investigate dynamic hydrate accumulation processes in water-dominated geological settings at the field scale. Acknowledgments: The first author would like to thank Isaac K. Gamwo for providing the simulator HydrateResSim (HRS) used for numerical model verification in the present study. The authors thank the four reviewers for their constructive comments that supported improving the quality of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The following nomenclature is used in this manuscript: The continuity equation for the mobile components in the fluid flow and transport simulator [63] is represented by In Equation (A1), φ is the effective porosity of the porous medium (hydrate-bearing sediment); P stands for the pore fluid pressure; and α and β are the compressibilities of the porous medium and fluid, respectively. Source/sink terms are represented by W.
The velocity of the mobile components, − → v , is defined to obey the single-phase Darcy's Law where the effective permeability of hydrate-bearing sediment is described by k; µ f is the dynamic fluid viscosity and − → g is the gravity vector. The density of the mobile components, ρ f , is expressed as where X w , X m , and X i are the mass fractions of water, dissolved CH 4 and inhibitor (NaCl here), respectively. Accordingly, ρ w , ρ m , and ρ i are the densities of water, dissolved CH 4 and inhibitor, respectively. In particular, the fluid density change due to dissolved CH 4 can be ignored (ρ m = 0). Species transport by diffusion and advection is described by the mass balance equation [63]: In Equation (A4), the concentration of each mobile component is stored in the matrix of C; D represents the hydrodynamic dispersion tensor, and the source/sink term is given by Q.
Conductive and convective heat transport is taken into account by the energy balance equation [63], written as where c p f is the specific heat capacity of the mobile components, and H is the source/sink term. The average thermal conductivity of immobile and mobile components, λ a (W·m −1 ·K −1 ), is defined as where λ f is the thermal conductivity of the mobile components, and the thermal conductivity of the immobile components, λ r , is expressed as In Equation (A7), λ s and λ h are thermal conductivities of the matrix of the hydratebearing sediment and CH 4 hydrate, respectively; S h is the CH 4 hydrate saturation of the pore space in the hydrate-bearing sediment.
The specific heat capacity of the immobile components, c pr (J·kg −1 ·K −1 ), is where c ps and c ph are the specific heat capacities of the matrix of the hydrate-bearing sediment and CH 4 hydrate, respectively. The density of the immobile components, ρ r (kg·m −3 ), is where ρ s and ρ h are densities of the matrix of hydrate-bearing sediment and CH 4 hydrate, respectively.
To solve the aforementioned governing equations, additional equations restricting the behaviour of the related components are required. The conservation relation of mass fractions of the mobile components is and the saturation summation of all components in the pore space is equal to 1: The effective porosity, φ, of the hydrate-bearing sediment is proposed by Spangenberg [76] as where ϕ is the intrinsic porosity of the hydrate-bearing sediment. Based on the assumption of a pore filling hydrate formation mechanism [41], the effective permeability, k (m 2 ), is assumed to obey the modified Carman-Kozeny relation [77] and is defined as a function of hydrate saturation; that is In Equation (A13), κ is the intrinsic permeability of the hydrate-bearing sediment; n is the linear relation with respect to the hydrate saturation by n = 0.7S h + 0.3 [27,76,78].