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

Abstract

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₄), which migrated from deep conventional hydrocarbon reservoirs, very likely through the present fault systems. Therefore, it is assumed that fluid flow transports dissolved CH₄ 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 first simulations that consider permafrost formation and thawing, as well as the formation of GH accumulations sourced from the upward migrating CH₄-rich formation fluid. The simulation results show that temperature distribution, as well as the thickness and base of the ice-bearing permafrost are consistent with corresponding field 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₄ 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 Amplification effect at the Mallik site from the early 1970s to the mid-2000s.

1. Introduction

1.1. Arctic Permafrost-Associated Gas Hydrates

Gas hydrates (GHs) are ice-like crystals storing enormous quantities of gas molecules such as methane (CH₄), 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₄ saturation were observed in core specimens sampled within the gas hydrate stability zone (GHSZ) [2]. The specifics of the hydrate-forming CH₄ 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₄ 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.

Permafrost, with near-surface temperatures below 0 °C for a minimum of two years [6], widely exists in the Canadian Arctic, and thus also in the Mackenzie Delta (MD) region. As illustrated by Figure 1, high-latitude terrestrial permafrost in general comprises a seasonally-thawed active layer underlain by perennially-frozen ground, often accompanied by thick thaw bulbs (taliks) overlain by thermokarst lakes [7,8]. Three factors mainly control the thermally equilibrated thickness of the permafrost and subsurface temperature gradient: (1) the mean annual temperature at the base of the active layer, (2) the thermal conductivity of sediments and (3) the basal heat flow [9]. During the recent glacial–interglacial history and regressive/transgressive marine cycles, permafrost forms or degrades, while the latent heat of the water–ice phase changes is consumed or released over considerable amounts of time [10,11]. In general, permafrost evolution is accompanied by the formation of permafrost-associated GHs, including the shallow pan-Arctic continental shelves and coastal plains [12,13] and the Qinghai–Tibet Plateau [14]. Permafrost and permafrost-associated GHs in the MD may exist since about 1.0 MaBP. Moreover, the current subsea relict permafrost formed terrestrially during late Pleistocene low sea-level stands [15]. Under contemporary global warming events, the relict permafrost degradation is continually promoted [7], whereby thick ice-bonded permafrost may prevent GHs from contributing greenhouse gases to the global carbon budget [16].

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₄ 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₄ 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₄ 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].

1.2. 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₄ venting and predicted the degradation of Arctic permafrost-associated GH reservoirs in terms of their contribution to the global CH₄ budget [7]. Additionally, they employed their models to verify if submarine groundwater discharge laterally transports dissolved CH₄ 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.

1.3. 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₄ 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₄ recycling; M5—buoyancy-driven CH₄ gas flow; M6—in situ conversion of gaseous CH₄ accumulations into GH deposits as the GHSZ base shifts downward below the top of the gaseous CH₄ accumulations. Moreover, this classification should be extended by including the mechanism of GH formation via CH₄ transported in the dissolved state (M7), as suggested by Frederick and Buffett [30,31] and by Egeberg and Dickens [34].

According to the migration distance of the source gas, these mechanisms can be further classified into short-range [3,5] (M1, M2, M3 and M6), medium-range (M4) and long-range migration mechanisms (M5 [5,33,35] and M7 [31]). The major type of source gas origin used in the abovementioned CH₄ hydrate formation and accumulation models can be mainly either thermogenic (M4, M5 and M7), mixed (M6) or biogenic (M1, M2, M3, M4 and M5).

Almost all explored marine CH₄ hydrate-bearing sediments can be described by the first five abovementioned mechanisms. M6 was defined to interpret the observed CH₄ hydrate distributions acquired from the sub-permafrost sediments at the Mount Elbert prospect [36,37] (cf.

Table 1. Summary of numerical models used to investigate GH accumulation mechanisms and comparisons of deposit features characterizing the evolution pattern of GH occurrences.

Drilling Site	GMGS1-SH2	GMGS3-W19	Beaufort Shelf	GC 955	Mount Elbert-01
Location	Shenhu area, slope of South China Sea	Shenhu area, slope of South China Sea	Canadian Beaufort Shelf	Green Canyon, Gulf of Mexico	Alaskan North Slope
States of source gas	Dissolved CH4	Gaseous CH4	Dissolved CH4	Gaseous CH4	Gaseous CH4
GH layer conditions	GH	GH	Permafrost-associated GH	GH underlain by free CH4 layer	Sub-permafrost GH
Saturation distribution along GH intervals	0–25% (60–250 mbsf), peak at 48% (210–220 mbsf)	0–45% (130–160 mbsf)	∼ 3%	79–93% (413–442 mbsf)	50% in average with peak at 75% (615–628 and 650–662 mbgl)
Simulation time period	1.5 Ma	30 ka	100 Ma	5.6 ka	Not applied
Model dimension	1D	1D	1D/2D	1D/2D	1D
GH accumulating mechanism	GH sedimentation and formation (in situ methanogenesis and M3)	M5	Submarine groundwater discharge transporting dissolved CH4 (M7)	M5	M6
Reference	[39]	[38]	[30,31]	[5]	[36]

). 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₄ 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₄ hydrate deposits with thick, high-saturation GH intervals in marine environments, they do not apply to concentrated terrestrial sub-permafrost 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₄ supply. Instead, Frederick and Buffett [30,31] proposed a new hypothesis considering fluid circulation below the permafrost transporting dissolved CH₄ into the GHSZ (referred to as M7 in this study) and promoting the formation of an 800 m-thick permafrost-associated CH₄ 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₄ hydrate systems does not match the field observations without assuming that the dissolved CH₄-rich fluid flows through the underlain fault systems into the targeted GHSZ of the Mallik site.

1.4. Hypothesis and Main Objectives of the Present Study

Conclusively, the mechanism controlling the formation of sub-permafrost CH₄ hydrate-bearing 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₄ upward and out of the overpressurized zones via pre-existing polygonal fault systems, serving as feed gas migration pathways. Then, the dissolved CH₄ 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.

2. Materials and Methods

2.1. Geological Settings at the Mallik Site

2.1.1. GH Research and Industry Exploration Wells

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. Coordinates of the Mallik wells and their distances to the Mallik J-37 well serving as the projection origin.

Wellsite	Coordinates (Latitude, Longitude)	Distance	Reference
Mallik J-37	69°26${}^{\prime }$38${}^{″}$ N, 134°38${}^{\prime }$23${}^{″}$ W	-	[44]
Mallik 2L-38	69°27${}^{\prime }$40.7${}^{″}$ N, 134°39${}^{\prime }$30.4${}^{″}$ W	2034 m
Mallik 3L-38	69°27${}^{\prime }$38.3${}^{″}$ N, 134°39${}^{\prime }$41.6${}^{″}$ W	2047 m
Mallik 5L-38	69°27${}^{\prime }$39.5${}^{″}$ N, 134°39${}^{\prime }$38.3${}^{″}$ W	2066 m
Mallik A-06	69°25${}^{\prime }$01${}^{″}$ N, 134°30${}^{\prime }$16${}^{″}$ W	6071 m
Mallik P-59	69°28${}^{\prime }$49${}^{″}$ N, 134°42${}^{\prime }$45${}^{″}$ W	4940 m

Note: Mallik J-37 well directly crosses the transect A-A’ of the 2D seismic profile 85987 adapted from [45], whereas all other wells are projected onto the transect A-A’.

.

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].

2.1.2. 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.

2.1.3. 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.

2.1.4. 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₄ supply [41] at the Mallik site, which suggests that the formation mechanism M7 is valid here instead of M6 (gaseous CH₄ 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₄ 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 mud-gas 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₄ in dissolved form as the feed gas to form hydrate at the crest of the Mallik Anticline.

2.2. Conceptual Model for Validating the Formation Mechanism of Sub-Permafrost GH Deposits

2.2.1. 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₄-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₄-bearing formation fluids. However, it is not yet fully understood to what extent these faults control the subsurface temperature distribution and CH₄-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.

2.2.2. 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 thermally-matured 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₄) 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 GH-bearing 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₄-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.

2.3. Numerical Model Implementation

2.3.1. Mathematical Model and Modeling Assumptions

A framework of equations of state (EOS) for equilibrium CH₄ hydrate formation has been integrated with the flow and transport simulator TRANSE [70], referred to as T_plusH (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_plusH 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:

• 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;
• 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;
• 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₄-rich fluid flows into the modeling domain via the two high dipping-angle Faults 3 and 4, shown at the bottom of Figure 3;
• 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.

2.3.2. Model Geometry

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,000 m in horizontal and 2000 m in vertical direction). It comprises two sequences with the respective lithology markers L = 0 and L = 1, and 17 faults (L = 2–18). The lateral model boundaries are represented by Dirichlet boundary conditions in terms of transported species concentrations and p-T terms to mimic an infinite aquifer. Moreover, the upper model boundary is represented by a Dirichlet boundary with fixed p-T conditions, overlying the aforementioned impermeable sediments. Finally, constant pressure and basal heat flow are applied at the model bottom, excluding the grid elements representing Faults 3 and 4. These faults represent source elements for the inflowing CH₄-rich formation fluid, with constant temperature conditions and a pre-determined fluid flow rate. The grid element thicknesses along the lateral and vertical directions equal 60 m and 10 m, respectively. In total, the 2D modeling domain consists of $200\times 200=40,000$ elements.

2.3.3. 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. 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	${\kappa }_x={10}^{-5}$, ${\kappa }_z={10}^{-7}$	Darcy	Assumed
Intrinsic permeability of K-R sequence	${\kappa }_x=5$, ${\kappa }_z=1$	Darcy	Assumed
Intrinsic permeability of Faults 2–18	${\kappa }_x=50$, ${\kappa }_z=10$	Darcy	Assumed
Intrinsic porosity of sediment matrix	0.3	-	[45]
Salinity of pore fluid	10	kg m−3	[45]
Initial pore pressure gradient	$1.012\times {10}^4$	Pa m-1	[75]
Density of sediment grain	2650	kg m−3	[45]
Hydration number	6.1	-	[76]
Thermal conductivity of wet sediment	2.45	W m−1 K−1	[75]
Thermal conductivity of CH4 hydrate	0.68	W m−1 K−1	[77]
Specific heat of sediment matrix	830	J kg−1 K−1	[77]
Specific heat of CH4 hydrate	2100	J kg−1 K−1	[77]
Diffusion coefficient	$1.0\times {10}^{-10}$	m2 s−1	[71]
Density of inhibitor (NaCl)	2160	kg m−3	[78]
Compressibility of porous medium	$1.0\times {10}^{-10}$	Pa−1	[71]

. 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.

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. List of applied initial and boundary conditions in the simulation.

Parameter	Value	Unit	Reference
Simulation time of the first stage	997	ka	Modified after [15]
Inflowing fluid rate	$6.0\times {10}^{-10}$	kg m−3 s−1	Assumed
Dissolved CH4 concentration	2.1	kg m−3	Assumed
Subaerial surface temperature since 1 MaBP	−16	°C	Modified after [27]
Basal heat flow	57	mW m−2	[27]
Freezing point of pore fluid	−1.5	°C	[27]
Frozen point of pore fluid	−2.5	°C	[27]
Specific latent heat of water-ice phase change	333.6	kJ kg−1	[77]
Transitional boundary condition used since 3 kaBP
Simulation time of the second stage	3	ka	Modified after [15]
Current subaerial surface temperature	−8.3	°C	Modified after [27]

, 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 Table 2 and Table 3. Unless otherwise stated, the employed parameters are not changed after the implementation of the transitional boundary conditions.

3. Simulation Results

3.1. Regional Distribution of Permafrost and Sub-Permafrost GH Accumulations

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 Figure 2b and Figure 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² [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. Comparison of the well-log-inferred data with the simulated key parameters of the GHSZ and permafrost distributions as well as GH saturation profiles.

Parameter/Well	P-59		L-38		J-37		A-06
	Observed	Modeled	Observed	Modeled	Observed	Modeled	Observed	Modeled
Depth of ice-bearing permafrost base (mbgl)	638	630	605	600	615	620	632	630
Depth of GHSZ base (mbgl)	1200 ± 100	1240	1100 ± 100	1130	1300 ± 100	1220	1300 ± 100	1280
Total thickness of sub-permafrost GH intervals (m)	156.8	150	213.1	230	292.2	300	182.8	230
Peak Sh within GH intervals (%)	-	5	90	91	-	15	-	18

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.

) 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.

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.

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 sand-dominated 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₄-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.

3.2. 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.

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 sub-permafrost 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. Comparison of DTS-logged data [81] with the simulated key parameters of the permafrost and GHSZ.

Parameter (mbgl)/Well	Observed			Simulated
	3L-38	4L-38	5L-38	J-37	L-38
Depth of ice-bearing permafrost base	599.4 ± 3.5	604.4 ± 3.5	600.4 ± 3.5	620	600
Depth of GHSZ base	1104 ± 3.5	1105 ± 3.5	-	1210	1130

Note: Observations of Mallik 3L/4L/5L-38 wells are originally obtained from [81] with the units of ‘mKB’, which were converted into ‘mbgl’ by compensating the relative depth to rotary kelly-bushing which is ca. 4.6 m above the ground level. The dash represents field observation data that could not be derived from [81].

lists the depths of the GHSZ and ice-bearing permafrost bases.

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.

4. 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.

4.1. Factors Impacting the GH Saturation Distribution

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.

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 Figure 5e and Figure 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 Figure 2b and Figure 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.

4.2. 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 near-surface 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. Comparison of permafrost heating and its projection on the increased near-surface annual amplitude temperature of permafrost at the Mallik site from 1969 to 2004.

Mode	Permafrost Heating Rate	Permafrost Type	Projected Temperature Increment	Reference
Biskaborn	0.39 ± 0.15 (K/decade)	Continuous permafrost	1.365 ± 0.525 K	[86]
Smith	0.3 (K/decade)	Discontinuous permafrost	1.05 K	[85]
AA	Two to threefold rate of the global temperature increment average	Discontinuous/continuous permafrost	0.88–1.32 K	[86,87]

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], while the DTS-logged temperature profiles were recorded in September 2003 (Figure 6a). The above DTS-logged ground temperature range of −7.0 to −7.5 °C (2003) is consistent with the near-surface ground temperatures recorded from 2004 to 2007 at the Mallik site (Figure 7b), considering a non-negligible AA-enhanced temperature increment of ca. 0.2 to 0.3 K within these three years.

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.

5. 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₄-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:

• The feed gas (thermogenic CH₄) 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₄ hydrate formation method [71] at field scale for the sub-permafrost GH accumulations at the Mallik site.
• 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.
• 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.
• At the Mallik site, the sub-permafrost GH accumulations did not release significant CH₄ 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.