Numerical Simulation of Coastal Sub-Permafrost Gas Hydrate Formation in the Mackenzie Delta, Canadian Arctic

: The Mackenzie Delta (MD) is a permafrost-bearing region along the coasts of the Canadian Arctic which exhibits high sub-permafrost gas hydrate (GH) reserves. The GH occurring at the Mallik site in the MD is dominated by thermogenic methane (CH 4 ), which migrated from deep conventional hydrocarbon reservoirs, very likely through the present fault systems. Therefore, it is assumed that ﬂuid ﬂow transports dissolved CH 4 upward and out of the deeper overpressurized reservoirs via the existing polygonal fault system and then forms the GH accumulations in the Kugmallit–Mackenzie Bay Sequences. We investigate the feasibility of this mechanism with a thermo– hydraulic–chemical numerical model, representing a cross section of the Mallik site. We present the ﬁrst simulations that consider permafrost formation and thawing, as well as the formation of GH accumulations sourced from the upward migrating CH 4 -rich formation ﬂuid. The simulation results show that temperature distribution, as well as the thickness and base of the ice-bearing permafrost are consistent with corresponding ﬁeld observations. The primary driver for the spatial GH distribution is the permeability of the host sediments. Thus, the hypothesis on GH formation by dissolved CH 4 originating from deeper geological reservoirs is successfully validated. Furthermore, our results demonstrate that the permafrost has been substantially heated to 0.8–1.3 °C, triggered by the global temperature increase of about 0.44 °C and further enhanced by the Arctic Ampliﬁcation effect at the Mallik site from the early 1970s to the mid-2000s. sub-permafrost GH accumulations are entirely different. Therefore, the present study focused on the numerical validation of the proposed formation mechanism of GH accumulations at the Mallik site to test further hypotheses on GH formation and investigate the impact of climate change on the stability history of sub-permafrost GH deposits.


Arctic Permafrost-Associated Gas Hydrates
Gas hydrates (GHs) are ice-like crystals storing enormous quantities of gas molecules such as methane (CH 4 ), and are possibly found in marine and permafrost-associated sediments, including gas-free and water-saturated Arctic sands [1]. GHs can form when favorable p-T conditions are met, and hydrate-forming gas is supplied either in free gas or supersaturated dissolved form. In Arctic permafrost-associated sediments, a series of continuously and uniformly distributed GH intervals of high CH 4 saturation were observed in core specimens sampled within the gas hydrate stability zone (GHSZ) [2]. The specifics of the hydrate-forming CH 4 supply (i.e., mass flow and physical state as well as other involved species) and their transport mechanisms [3] affect the heterogeneous characteristics of GH occurrences (i.e., their thickness and spatial distribution in the deposit). However, the hydrodynamics of the subsurface system related to fluid flow and CH 4 transport, bridging the pool of deep hydrocarbon gases and shallower GH occurrences, cannot be intuitively identified [4,5]. Therefore, numerical simulation is an irreplaceable method to study the hydrodynamics of subsurface systems, especially for sub-permafrost-associated GH-bearing strata.  [17]. A portion of the plain has been inundated by sea by thawing continental ice sheets at high northern latitudes. The thermogenic gaseous CH 4 migrates along the faults (red dashed lines) from deeper geological units (i.e., source rock) and accumulates in the upper sequences (e.g., Kugmallit-Richards Seq., Seq. = sequence, not to scale), which act as conventional gas reservoirs. Then, the accumulated CH 4 is transported upward into shallower sequences (e.g., Iperk-Mackenzie Bay Seq.), supporting the formation of GH occurrences.
After experiencing subaerial exposure during glacial episodes, permafrost-associated GHs [18] occur within and below the permafrost and are thus defined as intra-permafrost and sub-permafrost GHs, respectively. Dallimore and Collett [19] studied a sample containing intra-permafrost GHs recovered at a depth of 336.4 m within the ice-bonded permafrost near the Taglu gas accumulation in the MD, which proved intra-permafrost GH occurrence for the first time. Although the source gas composition of GH-bearing sediments varies, the CH 4 source of intra-permafrost GHs is likely derived from the degradation of in situ organic matter by microbial activity (methanogenesis) at low temperatures [19]. Instead, the primary source of sub-permafrost GH accumulations is CH 4 migrating from deep conventional hydrocarbon deposits through specific tectonic features (i.e., faults) [4] acting as fluid flow pathways. This hydrate-forming source gas originates from the deeper sedimentary sequences (source rocks of feed gas) and has been formed by thermogenic activities [20,21], such as thermal cracking [22]. Thermogenic CH 4 has not only been found within Arctic circumpolar GH accumulations in the Beaufort Sea and Russia [8], but also at the thermogenic basins lying between upper continental slopes and deep water sags, such as the Gulf of Mexico [23] and the South China Sea [24,25].

Review of Previous Numerical Studies on the Spatial-Temporal Evolution of Permafrost and Permafrost-Associated GHs in the Mackenzie-Beaufort Region
Present knowledge on the paleo-evolution of permafrost is mostly based on numerical experiments with the implementation of simple 1D and 2D geothermal models [13,15,[26][27][28][29], considering thermal effects to describe the ice-water phase transition (i.e., latent heat). Recently, Frederick and Buffett [7] established a more sophisticated multiphase fluid flow model to investigate the role of taliks as a potential pathway for gaseous CH 4 venting and predicted the degradation of Arctic permafrost-associated GH reservoirs in terms of their contribution to the global CH 4 budget [7]. Additionally, they employed their models to verify if submarine groundwater discharge laterally transports dissolved CH 4 to the Mackenzie-Beaufort (MB) shelf to form GHs [30,31].
Using the numerical simulator TEMP/W [32], Taylor et al. [26] established several 1D geothermal models to compare the sensitivity of relict permafrost and subsea GHs to climate change. Taylor et al. [27] further employed a 2D geothermal model to study the permafrost history under the inferred paleo-environment of the MB shelf and slope. In 2012, Majorowicz et al. [28] developed a series of 1D geothermal models to simulate profiles of permafrost and GHSZ bases, and validated these models by matching numerical results to field data. In addition, they studied GH and permafrost stability histories in the context of climate change [15] in the MB area. Moreover, Majorowicz et al. [29] utilized 1D and 2D geothermal models to predict the talik, permafrost and permafrost-associated GH histories in the MB region and concluded that the primary control in talik formation is lithology.

Comparison of Numerical Validations of Different Source Gas Migration Mechanisms in the Context of GH Formation and Accumulation
From the previous quantitative studies reviewed by You et al. [33], CH 4 hydrate formation models are categorized into six classifications by their respective gas migration mechanisms: M1-local biogenesis with/without diffusion; M2-local diffusion of dissolved gas; M3-advection with upward fluid flow and/or diffusion; M4-CH 4 recycling; M5-buoyancy-driven CH 4 gas flow; M6-in situ conversion of gaseous CH 4 accumulations into GH deposits as the GHSZ base shifts downward below the top of the gaseous CH 4 accumulations. Moreover, this classification should be extended by including the mechanism of GH formation via CH 4 transported in the dissolved state (M7), as suggested by Frederick and Buffett [30,31] and by Egeberg and Dickens [34].
Almost all explored marine CH 4 hydrate-bearing sediments can be described by the first five abovementioned mechanisms. M6 was defined to interpret the observed CH 4 hydrate distributions acquired from the sub-permafrost sediments at the Mount Elbert prospect [36,37] (cf. Table 1). In addition, M1 can be applied to describe the intra-permafrost hydrate occurrences discussed in [19]. Through comparisons of the GH accumulation time period and its spatial hydrate saturation (S h ) distribution (Table 1), the formation rate of high-saturation GHs at the GMGS3-W19 site [38] is presumed to be about 2-3 orders of magnitude higher than of those at the GMGS1-SH2 and GC 955 sites [39]. This difference in GH formation rates implies that the GH-forming gas at the GMGS3-W19 and GMGS1-SH2 sites is transported in the free-gas (M5) and dissolved (M7) states, respectively. Furthermore, it demonstrates that the notable complexity of GH-bearing environments and the diversity of GH formation mechanisms exist even in adjacent regions in approximately 14 km distance [40] in the Pearl River Mouth Basin, South China Sea. In comparison, the coastal sub-permafrost GH fields of the Alaskan (e.g., Mallik site [41]) and Canadian (e.g., Mount Elbert site [36,42]) Beaufort Sea likely experienced similar heterogeneous distributions of different GH formation mechanisms. To date, most well-studied sub-permafrost GH accumulations are located along the Arctic shoreline and contain high CH 4 hydrate saturations in the pore space of the respective sediments (e.g., Mallik [41] and Mount Elbert [36,42]). Although M4 [43] and M5 [35] can explain and predict the formation of CH 4 hydrate deposits with thick, high-saturation GH intervals in marine environments, they do not apply to concentrated terrestrial subpermafrost GH accumulations. Boswell et al. [42] compared the GH occurrence zone characteristics of the Mallik and Mount Elbert sites and confirmed that the interpretation of M6 does not apply to the Mallik site, whose intermediate and bottom GH intervals result from the lack of gaseous CH 4 supply. Instead, Frederick and Buffett [30,31] proposed a new hypothesis considering fluid circulation below the permafrost transporting dissolved CH 4 into the GHSZ (referred to as M7 in this study) and promoting the formation of an 800 m-thick permafrost-associated CH 4 hydrate occurrence with GH saturations of about 3%. Their results reasonably explain the existence of intra-permafrost GHs around the Taglu gas field as observed by Dallimore and Collett [19]. However, they do not achieve an agreement with the field data of the sub-permafrost GH intervals observed at the Mallik scientific drilling sites [2,9] (cf. Appendix A.1). Consequently, the spatio-temporal evolution of the permafrost and permafrost-associated CH 4 hydrate systems does not match the field observations without assuming that the dissolved CH 4 -rich fluid flows through the underlain fault systems into the targeted GHSZ of the Mallik site.

Hypothesis and Main Objectives of the Present Study
Conclusively, the mechanism controlling the formation of sub-permafrost CH 4 hydratebearing accumulations has not yet been sufficiently well understood. Thus, knowledge gaps still exist in terms of the interactions and interrelationships between the tectonic settings, as well as permafrost and permafrost-associated GH systems. The present study investigates the hypothesis that fluid flow transports dissolved CH 4 upward and out of the overpressurized zones via pre-existing polygonal fault systems, serving as feed gas migration pathways. Then, the dissolved CH 4 forms GH accumulations in the upper sequences, where the combination of large anticline systems, fault throws and sandy sediments act as hydrocarbon traps ( Figure 1). To our knowledge, a conceptual model has not yet been developed or quantitatively investigated based on the proposed formation mechanism of a sub-permafrost GH accumulation. Thus, the primary objective of this study is to validate the proposed mechanism by comparing the simulated spatio-temporal profiles of permafrost and permafrost-associated GH systems against temperature profiles and GH distributions observed at the Mallik site. Initially, the sub-permafrost GH occurrences in the MB area were inferred from seismic data obtained from conventional hydrocarbon exploration in the mid-1980s. Subsequently, the GH accumulations at the Mallik site have been subject to a series of international scientific joint investigations over the last decades. To date, four GH research and four industry exploration wells (Mallik P-59, L-38, J-37 and A-06) [44,45] at the Mallik site have encountered the underlain sub-permafrost GH sediments. Overall, two production research (Mallik 2L/5L-38) and two monitoring or water-reinjection wells (Mallik 3L/4L-38) have been drilled throughout the permafrost and its underlain GH intervals in the course of the numerous scientific Mallik field programs, as shown in Table A1.

Materials and Methods
According to previous studies [8,46], pan-Arctic permafrost-associated GHs are not as climate-vulnerable as the shallow marine GHs on the upper continental slopes, the subglacial GHs and the terrestrial plateau permafrost-associated GHs [14]. Although the relict permafrost degradation in the MD area has been continuously promoted by the contemporary climate change, leading to an increased probability of further GH destabilization, no significant amount of free gas has been monitored directly below the GHSZ base at the Mallik site [45,47]. Consequently, the decomposition of sub-permafrost GHs at the Mallik site has not yet been observed, as the GHSZ bottom is the region most sensitive to temperature changes induced by climate warming, as reported by Mestdagh et al. [46] and Tréhu et al. [48].

Available Seismic Data
Reflection seismic imaging is a feasible geophysical exploration technique for identifying GH deposits. High-saturation GH deposits usually correlate with observable, reverse-polarity bottom simulating reflectors (BSRs), high resistivity and seismic velocity (i.e., low acoustic travel-time differences) [46]. The appearance of BSRs indicates the high impedance contrast between the GH-bearing sediments with high seismic velocity and the underlying free or dissolved gas-charged sediments with low seismic velocity [49]. As a result, the BSR is regularly associated with the base of the GHSZ, and the higher the GH saturation (S h ), the more pronounced any BSR anomalies are [50].
In Figure 2, the seismic profile 85987 [45] crosses the Mallik anticline and runs to the north. Numerous normal faults penetrate the crest of this anticline with high dipping angles. Beneath and along the crest of the Mallik anticline, the seismic transect A-A' shows a number of high-amplitude reflectors indicated by the orange shading in Figure 2b, implying the sub-permafrost GH sediments. In addition, the phase boundary (p-T) marks the GHSZ base on the northern flank of the Mallik anticline. As plotted in Figure 2b, the Mallik anticline is bounded to the southwest trend by a large high dipping angle normal fault, which extends down into deeper stratigraphic sequences.

Uncertainties and Limitations of Seismic Reflection Measurements
The field parameters of seismic imaging were not optimal for evaluating shallow GH-related features, since these campaigns aimed at exploring deep conventional oil and gas deposits [45]. Additionally, ice-bearing permafrost and GH occurrences have similar seismic and electrical properties [31]. Therefore, BSRs cannot be expected within the ice-bearing permafrost [51], and it is nearly impossible to determine intra-permafrost GH by any standard surface-based geophysical observation technique. Despite these limitations, the reprocessing of industry seismic data [52] provides valuable insights into the dominating sedimentological and structural features on the distribution of inferred sub-permafrost GH accumulations.
In terms of the seismic characteristics of the Mallik site, the previously qualitative mapping of GH reservoirs, as well as the coarse resolution of the GH reservoir boundaries and thicknesses could not meet industrial production needs. In this context, Bellefleur et al. [48] addressed the lateral and depth-dependent geologic variations and lateral extent of the GH occurrences at further distances to the wells. In addition, Huang et al. [53] demonstrate that attenuation of seismic energy may be primarily attributed to scattering from small-scale heterogeneities and highly attenuate leaky mode propagation of seismic waves through larger-scale heterogeneities in sediments rather than the intrinsic attenuation of GH-bearing sediments.

Characteristics of the Phase Boundary at the Hydrocarbon Trap Bottom
The upper Cretaceous-Tertiary basin-fill and pre-Cenomanian rifted strata contain significant fossil fuel resources in the MD area. The combinations of large anticlines, fault-related structures and porous sandy sediments produce efficient hydrocarbon traps. The anticlines are large-amplitude structures with up to 3 km of relief and generally have rounded hinges [54], as shown in Figure 2b. They are thus commonly inferred to be structural traps containing the bulk of potential hydrocarbons. Most anticlines are asymmetric, with the steep branches usually facing basinwards, and many are cored with north-or northeast-verging thrust faults [9].
At the Mallik site, the observed high-saturation GH intervals occupy > 80% of the pore space of the unconsolidated clastic sediments, located at depths of 900-1100 mbgl and localized at the crest of the Mallik anticline [9]. Moreover, the seismic characteristics of the top, middle and bottom GH zones, identified by the A, B and C zones, respectively, were examined to assess the localized geologic influence on the spatial GH distribution. The phase boundary at the bottom of zone C coincides with the GHSZ base with no free gas below, further indicating a lack of gaseous CH 4 supply [41] at the Mallik site, which suggests that the formation mechanism M7 is valid here instead of M6 (gaseous CH 4 accumulations in situ conversion into GH deposits). As a comparison, Behseresht and Bryant [36] attributed a different S h distribution but similar vertical variation of the S h profile observed at the Alaskan Milne Point area to the discontinuous supply of gaseous CH 4 during the formation of two GH intervals via M6 instead of M7. The fact in support of this conclusion is that the mean value of S h significantly declined from ca. 60% to 15% at the bottom of the lower GH interval, which agrees with the GH distribution pattern generated by M6.
However, three well-log inferred GH zones [45] were confirmed by on-line mudgas monitoring [55], but the previous interpretation of the BSR as a free-gas zone below 1112 mbgl (the assumed GHSZ base shown in Figure 2b) was not identified in the mud gas analysis of the Mallik 5L-38 well. This observation is consistent with the seismic measurements [56], as well as the velocity and attenuation tomograms of the Mallik 5L-38 well [57], showing no clear indication of a free gas phase below the GH zone C. Moreover, Bellefleur et al. [47] improved the resolution of the surface 3D seismic data by compensation for attenuation effects of permafrost and GH-bearing sediments. They demonstrated the occurrence of a 5 m-thick free-gas interval located at ca. 75 m below the GHSZ base of the Mallik L-38 well and reported the absence of free gas within any GH intervals and directly underneath the p-T boundary at the Mallik site. According to MacKay et al. [49], super-saturated dissolved gas in the pore fluid may also produce a similar BSR event as the presence of a free-gas phase, which may explain the abovementioned BSR near the p-T boundary at the Mallik site. As a result, these findings support our assumption of CH 4 in dissolved form as the feed gas to form hydrate at the crest of the Mallik Anticline.

Faults as Feed Gas Migration Pathways
Faults are internally complex volumetric zones, playing a critical role as primary conduits for the convection of fluids beneath the permafrost, facilitating GH formation, gas migration and tectonics [58]. For example, many substantial high dipping angle faults occur in the major Taglu Fault Zone (TFZ) adopted from [17,59], as presented in Figure 2a, representing barriers to lateral fluid flow across the fault plane. However, fault damage zones may serve as preferential flow paths for the vertical migration of CH 4 -rich fluids originating at depth, causing temperature anomalies [59].
Some large high-dipping angle faults are depicted at the bottom of Figure 2b, with main GH accumulations above these faults, indicating that some faults are efficient pathways for CH 4 -bearing formation fluids. However, it is not yet fully understood to what extent these faults control the subsurface temperature distribution and CH 4 -rich fluid migration near the TFZ and how they impact the heterogeneity of GH accumulations at the Mallik site. Using the fault geometries ( Figure 3) derived from Figure 2b, the role of normal faults contributing to the formation of highly-saturated GH accumulations is investigated in the present study. Furthermore, we determine whether these faults are hydraulically conductive to fluid flow by altering their hydraulic properties (i.e., porosity and permeability) employed in the numerical model.

Origin of GH Feed Gases as well as Characteristics of the GHs and their Host Sediments
Characterization of the GH research wells in the MD via geochemical analyses, such as the Mallik 5L-38 well, shows that gases from sequences shallower than 500 mbgl are considered biogenic [60]. At depths of 550-850 mbgl, gases are mixtures of biogenic and thermogenic origin, while gases from more than an 890 mbgl depth, including the GH interval at depths of 890-1108 mbgl, are purely thermogenic. Generally, organic-and gas-geochemical analyses show that thermogenic gases probably migrated from thermallymatured sediments below 5000 mbgl, whereas low-maturity gases were generated from in situ organic matter, particularly the lignite layers of the Kugmallit Sequence [61]. Geochemical studies suggest that thermogenic gases have migrated up-dip along faults and stratigraphic boundaries, were adsorbed by lignite-rich strata and then intermixed with local microbial gases. In addition, archaeal DNA and microbial analyses suggest exceedingly low methanogenic archaeal populations and activities in deep sediments [62], supporting our simplification of the feed gas as thermogenic gas (CH 4 ) in the numerical simulations.
Based on core samples from the Mallik 2L/5L-38 wells, petrophysical analyses indicate GH accumulations of varying saturations in sand-rich sediments overlain by moderate to good fluid flow barriers [63]. As a result, the assumption of impermeable boundaries for caprocks is reasonable, given the relatively higher intrinsic permeability of the GHbearing sediments [64]. This also suggests that GHs occurred as a pore-filling type within a frame-supported pore-structure texture [65]. The Mallik GH deposits are controlled by the lithology of their host sediments, separated by fine-grained, silty facies and clay. Within the cored intervals located at depths from 890 to 1108 mbgl, coarse-grained sandstone sediments contain abundant GHs, whereas fine-grained sediments, i.e., the siltstone and mudstone, act as fluid flow barriers, showing minor to absent GH saturations.
The GH feed gas may have migrated into the GHSZ in numerous ways [66]. As a result, the migration of gas dissolved in formation fluids is required to form high-saturation GHs [50]. However, there is no valid evidence for the presence of a free gas phase within the GH zone (also referenced as the mixing layer) or directly below the GHSZ base, as discussed in Section 2.3. Therefore, our hypothesis on dissolved CH 4 -rich fluid migrating out of the overpressurized zone along preferential fluid flow paths (i.e., faults) [67][68][69] to the host sediments within the GHSZ is further supported, given the sediment below zone C shows a lack of free gas presence at the Mallik site.

Numerical Model Implementation 2.3.1. Mathematical Model and Modeling Assumptions
A framework of equations of state (EOS) for equilibrium CH 4 hydrate formation has been integrated with the flow and transport simulator TRANSE [70], referred to as T plus H (TRANSE + Hydrate) in a previous study [71]. In the current study, an EOS module describing the reversible processes of water-freezing and ice-thawing [10] has been implemented and integrated with T plus H to investigate the interactive processes of permafrost formation and degradation, as well as sub-permafrost GH formation.
To preserve the solution accuracy of the non-linear system of partial differential equations, we assume the following simplifications to maintain computational efficiency and numerical convergence requirements: 1.
The assumptions proposed by Li et al. [71] are considered, while temperature (−20-45°C) and pressure ranges (0-20 MPa) were extended in the updated EOS. Volume changes during ice-fluid and vice versa phase transitions are neglected; 2.
GH accumulations were categorized into various classes in previous simulation studies, depending upon varying boundary conditions. According to Moridis and Collett [72], Class-II deposits are defined as GH sediments overlain by impermeable rocks and underlain by aquifers without any presence of a free gas phase. In the numerical history matching for a six-day depressurization test, Uddin et al. [64,73] defined the three GH zones at the Mallik 2L-38 well as Class-II deposits. Moreover, data from specimen analysis suggests that processes of solute migration and GH formation at the Mallik site occurred in a semi-closed hydrodynamic system [74]. Therefore, this study employs the assumption of the presence of semi-closed Class-II GH sediments to investigate the formation of the sub-permafrost GH deposits overlain by impermeable upper sediments comprising permafrost and structural traps; 3.
Beneath the GH sequences at the Mallik site, fluid migration pathways are represented by faults of which the widths are unknown but assumed to be 60 m. The dissolved CH 4 -rich fluid flows into the modeling domain via the two high dipping-angle Faults 3 and 4, shown at the bottom of Figure 3; 4.
The depositional processes of the stratigraphic sequences and coastal erosion in the MD region, changes in sea level and tectonic subsidence or uplift since the late Pleistocene showed negligible influences on the evolution of sub-permafrost GH accumulations. Thus, these are not considered in the present modeling work. Figure 3 shows the lithology of the considered seismic cross section at the Mallik site as implemented into the numerical model, as well as its geometry and dimensions (12,

Model Parametrization
The hydrothermal properties of the sediments at the Mallik site were determined using an iterative history-matching procedure resulting in the parameters summarized in Table 2. The initially uniformly distributed hydrothermal properties of the permafrost and GH-bearing sediments change with the local increase in S h and ice saturation (S ice ), i.e., effective porosities and permeabilities will decrease in turn, as explained in Appendix A.2. Table 2. Hydrothermal properties of the permafrost and GH-bearing sediments as well as other parameters used in the present study.

Parameters Value Unit Reference
Effective permeability of I-M sequence κ x = 10 −5 , κ z = 10 −7 Darcy Assumed Intrinsic permeability of K-R sequence κ x = 5, κ z =  [71] According to the log interpretation and results from core analyses undertaken at the Mallik 2L and 5L-38 wells, the GH-bearing sediments are mainly composed of sands interbedded by shales, whereas the average effective porosity is between 0.24-0.4 and the mean permeability is about 2.9 mDarcy. Moreover, the highly saline residual pore fluid sampled at the GH intervals implies that fluid migration in the GH-bearing sediments at the Mallik site is slow due to the permeability reduction by GH accumulation and the in situ hydraulic properties (i.e., presence of interbedding shales).
In Table 3, the subaerial surface temperature is implemented at the top model boundary based on the arithmetic mean of the paleo-climate evolution derived from Taylor et al. [27] and Majorowicz et al. [15]. The listed initial fluid flow rate is the optimal input parameter derived from validating the respective simulated subsurface temperature distributions against the reported temperature profiles and anomalies in the MD region in [59] via an iterative modeling workflow. In addition, the thick permafrost and GH sequences may have persisted for about 1.0 Ma [15], but offshore permafrost formed during the Wisconsinan glacial period [9] when the Beaufort Shelf was exposed by seaward regression. However, the influence of sea-level changes is negligible in view of the deep buried onshore GH deposits residing under the thick permafrost at the Mallik site, and thus sea-level changes are not considered in the present study. Instead, a change in subaerial surface temperature from −16°C to −8.3°C is considered after the first stage at a simulation time of 997 ka to take into account the paleoclimate evolution. For that purpose, parameters reflecting average values reported in the literature are employed for model parametrization, as listed in Tables 2 and 3. Unless otherwise stated, the employed parameters are not changed after the implementation of the transitional boundary conditions.  Figure 4a shows the assessed boundaries of the seismically-inferred distribution of the two significant GH accumulations at the Taglu and Mallik sites on Richards Island, based on the seismic data reprocessed by Collet et al. [45] and initially acquired by the oil and gas industry in the mid-1980s. As depicted in Figures 2b and 4b, sub-permafrost GH deposits were seismically-inferred and numerically confirmed at the crest of the Mallik anticline with a connection to the deeply buried hydrocarbon GH feed gas via Faults 2-7 illustrated in Figure 3. Figure 4a shows that the GH deposit at the Mallik site is laterally continuously distributed and covers an area (green shaded) of ca. 51 km 2 [45]. The portion of the delineated GH deposit with the dark green shading denotes an area of higher GH saturations. Figure 4b suggests that GHs mainly accumulate in the upper host sediments of the Mallik anticline crest. As interpreted from Figure 5b, the simulation results (cf. Figure 5 and Table 4) are acquired from the Mallik J-37 well and projected locations of the industrial exploration wells Mallik P-59, L-38 and A-06 onto the 2D seismic transect A-A' to allow for model validation. Our simulation results reach a very good agreement with the well-logged and seismic observations [45] from the Mallik wells P-59, L-38, J-37 and A-06, which reveal highly variable S h of the Mallik GH deposit throughout the horizontal profiles shown in Figure 4a. The simulated S h profile presented in Figure 5b also confirms that the Mallik GH occurrences are highly concentrated around the location of the Mallik L-38 well, illustrated by the dark green shading in Figure 4a. Moreover, the simulated S h profiles (cf. Figure 5a,c,d) are consistent with the respective seismic-inferred GH concentrations from the Mallik P-59, J-37 and A-06 wells, located within the light green shaded delineation in Figure 4a.  Note: The well-log-inferred data are adopted from Collett et al. [45] and Dallimore et al. [9]. Well-log-inferred results that could not be derived from field reports [9,45] are represented by dashes. Figure 5e demonstrates that the S h profile of the GH interval is highly heterogeneous in line with the lithology along the borehole profile. It also indicates that coarse-grained sediments consisting of sand and pebbles contain abundant GH, whereas fine-grained sediments made of shale and silt host little to no GHs. This supports the conclusion that the S h distribution is lithologically controlled, as previously stated by Matsumoto et al. [80]. According to Katsube et al. [65], analysis of the petrophysical data (i.e., porosity, permeability and pore-size distribution) of the core samples acquired at the Mallik 5L-38 well support that relationship between S h and the lithological properties. It is known that silty and clayey sediments with relatively low porosities (0.25-0.3) have medium to very low permeabilities (generally 0.1-10 mDarcy). Overall, sediments of poor reservoir quality (i.e., low porosity and permeability) show a significantly lower potential to host GHs compared to those of good reservoir quality (i.e., high porosity and permeability) in the MD region.

Regional Distribution of Permafrost and Sub-Permafrost GH Accumulations
According to Boswell et al. [42], the intervals with high GH saturations (S h > 50%) exhibit high intrinsic permeabilities (often 1-5 Darcy), which are only known for sand-type and pebble-type sediments of high porosity (0.3-0.4). Consequently, a uniform porosity is employed in the simulations with varying permeabilities, rendering permeability the dominant factor of the lithology variations undertaken in the present study.
In Figure 5e, the relation between sediment characteristics and S h indicates that varying permeabilities in our simulations are the principal local constraint on pore occupancy within the GH columns. This also implies that the three GH zones occurring within the sanddominated intervals at depths of ca. 905-930 mbgl, 940-995 mbgl and 1065-1112 mbgl are preserved under the caprocks composed of shale and silt, interbedded with coal at depths of about 895-905 mbgl, 930-940 mbgl and 995-1070 mbgl, respectively. The presumption that fractures within the caprock connect those adjacent sand-rich sediments and supply dissolved CH 4 -rich fluid to facilitate GH formation until these clog the fractures is also supported by the previous findings. In contrast, the simulated S h profile at the Mallik L-38 well near the Mallik 5L-38 observation well demonstrates a relatively uniform GH distribution within the assumed homogeneous sandy host sediment, as shown in Figure 5b.
Throughout Table 4, the comparison of simulated and observed data implies that most of the numerically determined parameters, including the bases of the ice-bearing permafrost and peak GH saturations, match very well with the field observations with negligible deviations of up to 2.3%. In addition, all the simulated GHSZ bases are within the uncertainty range of the corresponding observations which amounts to about 9%. Although most of the simulated total GH interval thicknesses are consistent with the observed data with minor deviations of up to 8.5%, a relatively high deviation of ca. 26% is outlined between the observed and simulated results at the location of the Mallik A-06 well which was projected on the respective cross section. This deviation may be directly related to the fact that the Mallik A-06 well exhibits the greatest distance to that cross section, as presented in Table A1. Furthermore, it should be noted that the actual projection of the Mallik A-06 well is located outside of the green shaded area with a distance of about 500 m, as shown in Figure 4a.
In the MD, exhaustive studies have not yet been completed to survey the distribution of permafrost-associated GH. So far, the GH distribution observed by geophysical methods may be attributable to the fact that these measurements are easier to employ to explore large-scale and highly saturated GH deposits in the vicinity of hydrocarbon deposits, as noted by Frederick and Buffett [31]. As a result, it may be assumed that small-scale and less saturated GH deposits exist outside the green shaded areas in Figure 4a. Moreover, given the lack of extensive field observations on S h for each GH interval at the Mallik P-59, J-37 and A-06 wells, it is not feasible to conduct a quantitative analysis, but only qualitative comparisons of the peak saturations based on the limited available information.

Observed and Simulated Temperatures in the Vicinity of the Mallik L-38 Well
Subsurface p-T conditions control the stability of permafrost and sub-permafrost GH accumulations, while permafrost stability is much less susceptible to pressure perturbations (e.g., caused by regressive/transgressive marine cycles) but is more sensitive to temperature changes. Generally, a change in temperature at the ground surface will disrupt the subsurface temperature profile, which is originally in thermal equilibrium with the corresponding basal heat flow below the permafrost and the GH-bearing sandy columns. The transition of the temperature profile to a new thermal equilibrium condition is shown in Figure 6a, as the result of implementing transitional boundary conditions at a simulation time of 3 kaBP to mimic the paleo-climate evolution from the Late Holocene Epoch on [15,27] (cf. Table 3). As the arithmetic mean of subaerial surface temperature shifts from −16°C to −8.3°C, the simulation continues for a subsequent period of 3 ka until the final simulation time of 1 Ma is reached.  [81], respectively. Sequence boundaries are modified after [2]; (b) map of the near-surface ground temperature distribution in the Mackenzie Delta, in which the representative ground temperature data was collected from hydrocarbon exploration wells developed in the late 1960s and early 1970s, modified after [9]. With some modifications based on the map presented by Dallimore and Collett [9], Figure 6b shows the representative near-surface temperatures acquired at several hydrocarbon exploration wells. As indicated in Figure 6b, the current near-surface ground temperature ranges from −8°C to −9°C at the Mallik site, which matches the arithmetic mean of the subaerial surface temperatures of −8.3°C listed in Table 3. However, the subpermafrost parts (depth > 600 mbgl) of the simulated and observed temperature profiles demonstrate that the GH intervals (depth > 900 mbgl and < 1100 mbgl) are still located above the GHSZ base. This supports the conclusion that even though submergence of the shelf and coastal retreat occurred with rising sea levels since the Holocene, there has not been enough time to trigger significant permafrost degradation [9] and GH dissociation by the increase in temperature. Based on the observations and simulation results, Table 5 lists the depths of the GHSZ and ice-bearing permafrost bases. Table 5. Comparison of DTS-logged data [81] with the simulated key parameters of the permafrost and GHSZ.
When comparing the profiles in Figure 6a, the simulated sub-permafrost temperature profile of the Mallik L-38 well is consistent with the DTS-logged observations, while the simulated temperature at the Mallik J-37 well deviates from the observations of the Mallik 3L/4L/5L-38 wells by almost 2.0 K at the GHSZ base. This suggests an almost 100 m-deeper simulated base of the GHSZ at the Mallik J-37 well compared to the logged GHSZ bases at the Mallik 3L/4L/5L-38 wells and the simulated GHSZ base at the Mallik L-38 well, as indicated in Table 4. Furthermore, the relatively large variation between the location of the Mallik J-37 to L-38 wells matches the prediction of the subsurface temperature field presented by Chen et al. [59]. It has been demonstrated that subsurface temperature variation could reach about 5.0 K/km along the lateral direction at the given depth of 1100 mbsl. As addressed in Section 2.1.1, Mallik J-37 was a legacy well from previous industrial drilling activities, probably conducted in the 1970s. Unfortunately, temperature data for the Mallik J-37 well are not available in any existing open-access datasets and publications. Therefore, the simulation results are compared to the adjacent DTS-logged observations of the scientific research wells Mallik 3L-38, 4L-38 and 5L-38 in Figure 6a.
However, the intra-permafrost intervals of the simulated and observed temperature profiles demonstrate that the observed near-surface temperatures (depth < 50 mbgl) are higher than in the numerical predictions by up to 0.8 to 1.3 K. This can be accounted as the combined consequence of the global warming above the pre-industrial levels before 2004 [82] and an overestimation of the ice content within the pore space in the permafrost intervals [83] as discussed in the following section.

Discussion
As addressed in Section 1.3, notable variations in GH formation mechanisms have been identified in the vicinity of the same GH-bearing basins, even though their host sediments have similar lithologic properties. Likewise, although sediment properties of the Mallik and Mount Elbert sites are identical, the formation mechanisms of their sub-permafrost GH accumulations are entirely different. Therefore, the present study focused on the numerical validation of the proposed formation mechanism of GH accumulations at the Mallik site to test further hypotheses on GH formation and investigate the impact of climate change on the stability history of sub-permafrost GH deposits. Figure 2b shows that some potential GH-bearing intervals, inferred from the 2D seismic response and depicted in orange shading, accumulate along the Kugmallit-Mackenzie Bay Sequences boundary (KMB), and appear as several scattered belts above the KMB. However, since the scattered belt-like GH intervals above 890 mbgl (Figure 2b) were not confirmed by well logging obtained at the Mallik site, the corresponding sediments are parametrized to be impermeable caprocks as listed in Table 3. Likewise, Figure 4b indicates that one simulated GH accumulation is distributed along Fault 6 and above the Faults 7 and 8 with a thickness up to ca. 300 m. However, it is not found in the seismic interpretation shown in Figure 2b. This inconsistency suggests either a likely undiscovered GH accumulation located ca. 2 km south from the Mallik J-37 well near the Taglu GH accumulation or that Fault 7 was not hydraulically active since the Late Pleistocene. However, it is not possible to form any definite conclusions without additional field observations from the boreholes that penetrate the potential GH-bearing sediments in this area due to the fact that this region has not been the target of previous hydrocarbon exploration and scientific investigation activities. Therefore, it merits further attention to promote the mapping precision of GH resources, as additional field data become available.

Factors Impacting the GH Saturation Distribution
As discussed in Section 3.2, lithology and fault architecture control the spatial evolution of the S h profiles within the simulated and observed GH intervals. In addition, their common local constraint on pore occupancy (i.e., maximum S h ) within the GH columns is the sediment permeability. Although Figure 5e reflects the lithological heterogeneity along the GH-bearing column at the location near the Mallik L-38 well, it is not feasible to integrate this lithologic variation into the model geometry in the form of permeability variations due to the lack of field data covering the area presented in Figure 2a.
The variations in the sequence boundary depths between the projected locations of the Mallik L-38 and Mallik 2L/5L-38 wells are not negligible. As shown in Figures 5e  and 6a, the KMB of the Mallik 5L-38 well is at a depth of ca. 930 mbgl, while the KMB at the projected location of the Mallik L-38 well is at ca. 900 mbgl, as plotted in Figures 2b and 3. This limitation is expected to be overcome by implementing a 3D regional model geometry for subsequent simulation studies based on the available 2D seismic profiles, interpreted cross sections and well data.
In addition, Figure 4a also illustrates that some faults (Faults 12-17), located directly beneath the caprock and above Faults 3 and 4 contain the highest hydrate saturations (S h ≈ 92%) of the whole Mallik GH deposit. For comparison, the nearby region of Faults 12-17 contains S h < 70%. This regional S h variation strongly correlates with the intrinsic permeability difference between the faults and sediments (Table 2). Therefore, this clearly shows that the simulated and observed heterogeneities of S h distribution at the Mallik site are mainly dominated by structural geology and lithology.

Impact of Climate Change on Permafrost Heating
Permafrost heating mainly refers to the temperature increment of the ice-bearing sediment addressed in the present study. Permafrost has been significantly heated during the past 3-4 decades [84]. Smith et al. [85] conclude that heating rates are generally below 0.3 K per decade in sub-Arctic regions for warmer discontinuous permafrost with nearsurface ground temperatures close to 0°C. However, for cooler continuous permafrost with near-surface ground temperatures below −2°C, heating rates could be up to ca. 1 K per decade in the high-latitude Arctic areas. The Arctic Amplification (AA) phenomenon [86,87] well documents that the Arctic is warming at a twofold to a threefold rate of the global average, with some typical AA heating gradient trends discussed by Biskaborn et al. [86]. Their study demonstrates that the annual amplitude temperature of permafrost near the ground surface increased globally by 0.29 ± 0.12 K during 2007-2016, whereas an increase by 0.39 ± 0.15 K has been observed for the continuous permafrost such as the Mallik and Mount Elbert sites as listed in Table 6. Note: The average global temperature increment is obtained from the observation shown in Figure 7a, which is ca. 0.44 K from 1969 to 2004. Figure 6a shows that near-surface ground temperatures deviate from 0.8-1.3 K between the simulated and observed temperature profiles at depths < 50 mbgl. As interpreted in Section 3.2, this deviation could be caused by overestimating the ice content within the permafrost because the pore volume of the ice-bonded permafrost is not often fully saturated with ice [88] but contains large amounts of organic matter and minerals instead. Hence, these materials will not experience the consumption of latent heat in the process of permafrost-heating. This may cause the observed climate change-driven permafrost degradation to be substantially faster than in previous predictions [89], even without taking thermokarst-inducing processes [83] and AA effects [86,90] into account.
The projections of the current ground surface temperatures vary from ca. −7.0°C (Mallik 3L/5L-38 wells) to ca. −7.5°C (Mallik 4L-38 well), derived from the extrapolation of the DTS-logged temperature profiles, shown in Figure 6a. This inferred ground temperature range also approximately matches the DTS-measured near-surface temperature of −7.8°C at the Alaskan Ignik Sikumi field in 2009, which was recovered from drilling-induced temperature disturbance and cementing effects after 128 days of well completion [91]. Furthermore, the near-surface ground temperature distribution (Figure 6b) is inferred from the dataset of hydrocarbon exploration wells drilled between the late 1960s and early 1970s [9,92] According to Burn and Kokelj [92], the near-surface ground temperatures were generally increased by 2.0 K from the initial estimate of −8 to −9°C (Figure 6b) in the early 1970s to the observed −6 to −7°C (Figure 7b) at the Mallik site in the mid-2000s. During this period, global warming triggered a temperature increase by 0.55 K from ca. 0.29°C in 1968 to ca. 0.84°C in 2007 [82] (Figure 7a), indicating that the estimated permafrost heating intensity of the Mallik site is over three times higher than the global warming average, confirming the AA effect.
As shown in Table 6, the impact of global warming on the near-surface ground temperature of the continuous permafrost at the Mallik site is estimated to be in the range of 0.88-1.32 K from the late 1960s through the mid-2000s, considering the abovementioned AA trend. This forecasted temperature increment is in perfect agreement with the abovementioned deviation of 0.8 to 1.3 K between the observations and simulations, as shown in Figure 6a. This suggests that the temperature increase caused by climate change is the main reason for the deviations between the simulated and observed near-surface temperatures at the Mallik site. Moreover, it also implies that the bias sources are broader than the impact of climate change, owing to interactions between soil erosion, wildfire [93], vegetation and ground ice content [85]. According to Min et al. [94], the anthropogenic influence on the Arctic Sea ice became detectable already in the early 1990s. Thus, it can be deduced that global warming has significantly driven permafrost heating since the early 1970s at the Mallik site. Unfortunately, unlike glaciers and snow covers, the lack of data on the evolution of permafrost in temporal and spatial terms cannot yet be effectively compensated by remote sensing [95]. Consequently, our inference on the permafrost heating time period requires validation by additional field observations.

Conclusions
In this study, we have developed and applied a thermo-hydro-chemical model to study the spatio-temporal evolution of permafrost and the formation of sub-permafrost GH accumulations, facilitated by the upward migration of dissolved CH 4 -rich fluids flowing out of an overpressurized zone along the fault systems at the Mallik site. The simulated temperature profiles, the base of the ice-bearing permafrost, the thickness of the hydrate intervals and the peak S h are in very good agreement with corresponding field observations. Consequently, the above presented simulation results support us to make the following conclusions: 1.
The feed gas (thermogenic CH 4 ) transport mechanism proposed in this study has been validated with regard to its dissolved migration state and highly dipping faults as migration pathways. Simultaneously, our simulations prove the general feasibility of the previously addressed laboratory-scale CH 4 hydrate formation method [71] at field scale for the sub-permafrost GH accumulations at the Mallik site. 2.
GH-rich accumulations are generally formed in favorable geothermal environments, preserving thick GHSZ and occur in areas with moderate tectonic deformation intensity. Hence, the varying GH saturations observed along the well intervals can be attributed to the variability of the physical properties of the host sediments and lithology at the Mallik site.

3.
Our simulation results exhibit an undiscovered GH accumulation located ca. 2 km south of the Mallik J-37 well along the transect A-A', given that Fault 6 ( Figure 3) has been hydraulically active since the Mid-Pleistocene. Moreover, Faults 3 and 4 ( Figure 3) may still be hydraulically active in the present day, and thus deserve further attention and investigation in view of their role in the stability of the GH accumulations.

4.
At the Mallik site, the sub-permafrost GH accumulations did not release significant CH 4 in view of the global carbon budget under the contemporary global warming events [8] until the mid-2000s, because they are well preserved below the thick ice-bounded permafrost within the undisturbed GHSZ. Our numerical simulations evidently support field observations on permafrost heating at the Mallik site, induced by the AA-enhanced pan-Arctic climate warming, which has been observed since the early 1970s.
The presented numerical modeling framework can be applied for the calibration of geophysical measurements and the validation of the interpretation of BSRs derived from seismic reflection profiles at the reservoir scale. This framework is capable of investigating and projecting the potential gas hydrate enrichment and resource density within the target Arctic sandy sediments. Moreover, it supports improving the understanding of sub-permafrost groundwater migration and the evolution of hydraulic conductivity of geologic fault systems. In addition, the model can be used to quantify the AA phenomenon on the evolution of the mean annual near-surface temperature distribution in pan-Arctic permafrost regions. The 2D model used in this study will be extended to a 3D model in the near future, using available 2D seismic cross sections and well data to further quantitatively investigate the role of the detected geologic faults in the formation of the GH deposit at the Mallik site.  Acknowledgments: The first author would like to thank Benjamin Nakaten (GFZ German Research Centre for Geosciences) for providing the model pre-processing tool modelator, used to implement the seismic profiles into the 2D model geometry in the present study. The authors acknowledge the constructive comments of the anonymous reviewers, which supported improving the quality of the manuscript.

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

Nomenclature:
The following nomenclatures are used in this manuscript: The mass balance, energy balance and continuity equations employed in the scope of the present numerical simulation are presented by Li et al. [71] and Kempka [70]. With introducing the frozen state of pore fluid (i.e., ice) as an additional component, the following equations of state (EOS) were integrated into TRANSE. The volume fraction relations of interstitial pore components are defined as: where ϕ is intrinsic porosity, S h stands for the CH 4 hydrate saturation.
The terms of φ f , φ h and φ ice refer to the effective porosity available for the flow of the mobile component (i.e., fluid), as well as the pore volume fraction occupied by hydrate and ice, respectively. The constraint φ h + φ f + φ ice = 1 implies that pore space is fully saturated. Additionally, Θ denotes the unfrozen fraction of fluid within the pore space and is generally assumed to be a temperature-dependent piece-wise function in the interval of phase transition [96], such as In Equation (A4), T describes the mixture temperature in the pore space; T L is the temperature of the thawing/freezing point (liquidus, usually −1.5°C in the MB region [27]); and w denotes the phase transition interval, i.e., w = T L − T S (usually w = 1°C in the MB region [27]), where T S is the temperature of the frozen point (solidus).
The derivative of the above partition function Θ is integrated into the EOS to improve the convergence of the numerical solution; that is For fluid-saturated, frozen and CH 4 hydrate-bearing sediments, the energy balance equation is updated by replacing the term of effective porosity φ with φ h , then re-written as In Equation (A6), c pr is the specific heat capacity of the immobile components, including the sediment matrix, ice and hydrate, which is modified by introducing the ice component to the formula and taking the latent heat of ice formation, L [11], into account: In Equation (A7), c s , c ice and c h are the heat capacities of the sediment matrix, ice and CH 4 hydrate, respectively, while L is the specific latent heat of ice formation (cf. Table 3). The following terms are updated by adding the ice component to each formula.
The average thermal conductivity of the immobile and mobile components, λ a , is expressed as where λ f is the heat conductivity of the pore fluid (mobile components), while the heat conductivity of the immobile components, λ r , is described by where λ s , λ ice and λ h are heat conductivities of the sediment matrix, ice and CH 4 hydrate, respectively. The density of the immobile components, ρ r , is where ρ s , ρ ice and ρ h are the densities of the sediment matrix, ice and CH 4 hydrate, respectively.