Geologic controls on the genesis of the Arctic permafrost and sub-permafrost methane hydrate-bearing system in the Beaufort–Mackenzie Delta

The Canadian Mackenzie Delta exhibits a high volume of proven sub-permafrost gas hydrates that naturally trap a significant amount of deep-sourced thermogenic methane (CH4) at the Mallik site. The present study aims to validate the proposed Arctic sub-permafrost gas hydrate formation mechanism, implying that CH4-rich fluids were vertically transported from deep overpressurized zones via geologic fault systems and formed the present-day observed GH deposit since the Late Pleistocene. Given this hypothesis, the coastal permafrost began to form since the early Pleistocene sea-level retreat, steadily increasing in thickness until 1 Million years (Ma) ago. Data from well logs and 2D seismic profiles were digitized to establish the first field-scale static geologic 3D model of the Mallik site, and to comprehensively study the genesis of the permafrost and its associated GH system. The implemented 3D model considers the spatially heterogeneously distributed hydraulic properties of the individual lithologies at the Mallik site. Simulations using a proven thermo-hydro-chemical numerical framework were employed to gain insights into the hydrogeologic role of the regional fault systems in view of the CH4-rich fluid migration and the geologic controls on the spatial extent of the sub-permafrost GH accumulations during the past 1 Ma. For > 87% of the Mallik well sections, the predicted permafrost thickness deviates from the observations by less than 0.8%, which validates the general model implementation. The simulated ice-bearing permafrost and GH interval thicknesses as well as sub-permafrost temperature profiles are consistent with the respective field observations, confirming our introduced hypothesis. The spatial distribution of GHs is a result of the comprehensive interaction between various processes, including the source-gas generation rate, subsurface temperature, and the hydraulic properties of the structural geologic features. Overall, the good agreement between simulations and observations demonstrates that the present study provides a valid representation of the geologic controls driving the complex permafrost-GH system. The model’s applicability for the prediction of GH deposits in permafrost settings in terms of their thicknesses and saturations can provide relevant contributions to future GH exploration and exploitation.


Introduction
The Mackenzie Delta (MD) is a river-mouth depocentre, and the second largest Arctic delta (Forbes et al., 2022). It is the most economically accessible area along the Arctic coast of the Beaufort-Mackenzie Delta Basin (Dixon et al., 2019). The basin is an essential component of the Canning-Mackenzie deformed assessment unit (Houseknecht et al., 2020), also recognized as the Beaufort-Mackenzie tectono-sedimentary element (Chen et al., 2021). A high amount of methane stored as sub-permafrost gas hydrate (GH) resource has been deterministically appraised in the MD (Osadetz and Chen, 2010). GHs are ice-like crystalline solids consisting of hydrate-forming gases, such as methane (CH 4 ), trapped within the water molecules forming cageshaped structures (Sloan Jr and Koh, 2007). Permafrost refers to subsurface sediments exhibiting temperatures below 0°C for at least two consecutive years, regardless of the sediment composition (Woo, 2012).
Decades of industrial exploration (Taylor and Judge, 1976;Judge et al., 1981) and scientific research Dallimore et al., 2005;Dallimore et al., 2012) in this petroliferous region have produced a broad spectrum of geoscientific data and knowledge on the basin. The GH composition found at the Mallik sites, including the Mallik P-59, L-38, 2L-38, 3L-38, 4L-38, 5L-38, J-37, and A-06 wells (see Figure 1), is dominated by thermogenic CH 4 , which migrated from deep conventional hydrocarbon deposits located in the Taglu and Richards sequences illustrated in Figure 2 (Lorenson et al., 1999;Lorenson et al., 2005). According to the hypothesis addressed by Chen et al. (2008), the CH 4 -rich fluid is migrating upward out of the deeper hydraulic overpressure zone (Hu et al., 2018) through geologic gas-source faults to form CH 4 hydrates in the shallower Kugmallit Sequence (Figure 1) over geological times. When the CH 4 -rich fluid enters the Methane Hydrate Stability Zone (MHSZ), the dissolved CH 4 is converted from the supersaturated solution into the hydrate phase (Kashchiev and Firoozabadi, 2002). This hypothesis on the sub-permafrost GH formation mechanism is to be validated by the 3D numerical model developed in the scope of the present study.
A significant amount of GH was formed and preserved in the Kugmallit Sequence, mainly consisting of sandstones (Figure 2A) due to the presence of the Mallik anticline which is sealed by the Mackenzie Bay Sequence with a lithology mainly composed of shale (Huang, 2009). The GH presence was observed in seismic responses and well logs which are summarized in Table 1. In the far-field of the anticline, the pore space within the sandy sediment is occupied by CH 4 -loaded fluids rather than highly-saturated with GHs due to the lack of a geologic trap and the absence of CH 4 super-saturated formation fluids (Huang, 2009).
In view of previously undertaken modeling activities, earlier basin-wide 3D geothermal hydrocarbon flow models (Kroeger et al., 2008;Kroeger et al., 2009) to identify the origin of the hydrocarbon province and recent lithosphere-scale 3D structural models (Sippel et al., 2013;Sippel et al., 2015) on the temperature and maturity history did not investigate the genesis of GHs. Nevertheless, these models were an important basis for the present study, aiming to quantitatively assess the spatial distribution of GH occurrences in relation to the potential CH 4 migration pathways. Furthermore, available geothermal models focusing on permafrost evolution (Majorowicz et al., 2012a;Majorowicz et al., 2012b;Taylor et al., 2013;Majorowicz et al., 2015;Overduin et al., 2019) were exclusively realized in one and two dimensions and validated by only a limited number of bottom-hole temperature data acquired from the basin, neglecting the presence of sub-permafrost GHs. To date, Frederick and Buffett (2014), Frederick and Buffett (2015), Frederick and Buffett (2016) have presented a more sophisticated multiphase fluid flow simulator to account for permafrostassociated GH formation via dissolved CH 4 . In addition to GH formation, their simulator was also applied to investigate GH dissociation under present and future global warming events (Frederick and Buffett, 2014). Although their 2D simulation approaches could predict the intra-permafrost GH formation with GH saturations (S h ) of less than 3.5%, they did not simulate the formation of highly saturated sub-permafrost GH intervals. Besides, their approach is not capable of quantifying the influence of structural geological heterogeneity (i.e., anticlines, faults) on GH distribution.
As mentioned above, the integrated permafrost and subpermafrost GH reservoir system has not yet been sufficiently investigated. To promote the understanding of this complex system, we previously employed a 2D hydrogeological model (Li et al., 2022a) to study the spatial-temporal patterns of permafrost and the genesis of sub-permafrost GH accumulations facilitated by upward-migration of CH 4 -rich formation fluids. The proposed mechanism for the formation of highly saturated sub-permafrost GH via dissolved CH 4 -rich fluid has been quantitatively validated by comparison to observed temperature profiles and seismic-inferred GH distributions.
Although the previous 2D model (Li et al., 2022a) used basic structural geologic data, a more detailed mapping of the model geometry and comprehensive spatial structural analyses are required to refine the regional model. For instance, it has already been suggested before that many adjacent structures are known but have yet to be adequately investigated for their GH migration and trapping potentials (Lane, 2002). In addition, 3D geological models are significantly more suited for developing and visualizing geological knowledge than 2D cross-sections in highly heterogeneous subsurface environments (Thornton et al., 2018). By Frontiers in Earth Science 02 frontiersin.org

FIGURE 1
Overview map of the study area with the trend of employed seismic profiles, wells, and regional fault zones in the Ivik-Mallik-Taglu area. Map created using a base map from OpenStreetMap; Spatial reference: WGS 84/UTM zone 8N, modified after Chabab and Kempka (2023b).
applying various methods to reprocess the 3D seismic-reflection data acquired from the Mallik 5L-38 site in 2002 (Brent et al., 2005), the complexity of the sub-permafrost GH-bearing sediments has been visualized and demonstrated in a series of studies conducted by Bellefleur et al. (2006), Bellefleur et al. (2007), Riedel et al. (2009. These studies report that local heterogeneity within sedimentary rocks significantly influences the GH distribution at the Mallik site. Therefore, the need for implementing a simulation based on a 3D static model for quantifying the geologic controls on GH distribution is substantial. For instance, it is not feasible to parameterize specific faults near the model boundary as hydraulically impermeable in a 2D simulation, because this would turn a semi-closed model into a closed system, while a 3D model can easily overcome this obstacle by extending the simulation domain accordingly. Additionally, the projection of the Mallik wells onto a 2D transect does not necessarily represent the in-situ geothermal conditions and structural geology at the well locations, leading to less reliable simulation results. For a 3D model, the additional dimension allows a full spatial representation of the hydraulic dynamics and role of geologic faults in terms of the upward-migrating CH 4 -rich fluid, and thus a qualitative and quantitative assessment of the geologic controls on the spatial extent of heterogeneous sub-permafrost GH accumulations. Conclusively, the previous 2D model does not suffice to comprehensively evaluate the permafrost and GH-bearing system, but rather provides qualitative information on the general validity of the previously introduced hypothesis. Therefore, we present a 3D structural model covering approximately 38 onshore-km 2 of the MD in the present study. In addition to well logs (Table 1), published seismic interpretations Dallimore et al., 2005) were employed to construct a full-scale 3D model of the Mallik anticline GH trap. Our model is used to quantitatively assess the structural geologic controls provided by the hydraulically conductive fault and anticline systems as well as hydrocarbon traps on the temporal and spatial development of the sub-permafrost GH occurrences since the Late Pleistocene.   Ashford et al. (2012). The details on Mallik J-37, A-06, P-59 and L-38 are adopted from field reports, such as coordinates (Osadetz et al., 2005), drilling period (Taylor and Judge, 1976;Judge et al., 1981), and heat flux (Majorowicz and Smith, 1999). The information on Mallik 2L-38, 3L-38, 4L-38 and 5L-38 are adopted from field reports, such as coordinates (Ashford et al., 2012) and drilling periods .
2 Methods and input data 2.1 Geological setting and 3D structural geological model The Beaufort MD Basin is an important petroleum province due to its abundant petroleum resources. The Taglu gas field discovery in 1971 led to exploratory efforts primarily focusing on oil during the 1970s to the mid-1980s . For example, the Mallik L-38 well was drilled into a fault-bounded anticline. Initially, GH deposits were considered a drilling hazard in the exploration of deeper petroleum prospects. Once assumed to be rare, GHs are now expected to occur globally in vast volumes. Recent research activities have been aimed at advancing the potential for cost-effective usage of GH deposits as a reliable alternative energy resource and bridging technology during the energy transition from fossil fuels to renewables. In the MD, the first GH sample was cored from the permafrost at the Taglu field (Dallimore and Collett, 1995), and deliberate GH studies were inaugurated in 1998 by drilling the Mallik 2L-38 well , followed by the development of the Mallik 3L-38, 4L-38, and 5L-38 wells in 2002, as listed in Table 1.

Geology of the Beaufort Sea and Mackenzie Delta region
The MD is part of a rifted continental margin basin in the Canadian Arctic that was formed during the Early Cretaceous. Deposits in the delta and Beaufort Sea continental shelf cover the Frontiers in Earth Science 04 frontiersin.org geological period from the Paleozoic to the Holocene. Westwardthickening Paleozoic rocks, intersected by numerous faults, are overlain by the main post-rift basin-fill, a thick sedimentary succession formed by deltaic processes that marks a regional unconformity at the boundary between the Upper Cretaceous and pre-Cretaceous rocks (Dixon and Dietrich, 1988;Dixon et al., 1992;Collett et al., 1999). Marine organic-rich muds deposited during the Upper Cretaceous transition into younger Cenozoic deltaic sandstones and shales result from a series of sedimentation cycles of the progressing Laramide orogeny southwest of the present coastline (Dixon and Dietrich, 1988;Dixon et al., 1992;Miller et al., 2005). Modern deltaic sediments and older fluvial and glacial deposits cover the Mesozoic and Cenozoic strata, whose thickness tends to increase northwards towards the coast to 12-16 km below sea level . Major unconformities terminate each transgressive-regressive sequence within the Cenozoic strata (Dixon and Dietrich, 1988;Dixon et al., 1992;Collett et al., 1999;Miller et al., 2005;Chen et al., 2021). The GH-containing Cenozoic sequence within the MD in the Ivik-Mallik-Taglu area encompasses four distinct lithostratigraphic units. The Eocene Richards Sequence consists of fine-grained distal pro-delta and delta-slope deposits that comprise mostly mudstones and siltstones. These are unconformably overlain by coarser deltafront and delta-plain deposits from the Oligocene Kugmallit Sequence, which exhibits the thickest sediment deposition among the present lithostratigraphic units (Dixon and Dietrich, 1988;Dixon et al., 1992). Predominantly outer shelf and deep-water deposits comprising mudstones and siltstones define the Late Oligocene to Middle Miocene Mackenzie Bay Sequence, which lies unconformably above the Kugmallit Sequence and older strata. Near-shore sand and gravel facies transitioning to mud and silt deposits on the continental slope are characteristic of the thick and comparatively less deformed deposits of the Pliocene to Pleistocene Iperk Sequence, which unconformably covers underlying sequences and truncates most of the early to mid-Cenozoic structural features at its base (Dixon and Dietrich, 1988;Dixon et al., 1992;Osadetz et al., 2005). Cenozoic tectonics with extensional deformation phases, following the Mesozoic rifting due to the opening of the Canada Basin and compressional deformation during the Early Eocene to Late Miocene have overprinted the basin-fill in the MD. Normal, thrust and strike-slip faulting of the underlying strata together with a predominantly north-west trending folding resulted in several large-amplitude anticlinal systems, such as the Mallik anticline, which serve as efficient structural traps for the accumulation of GHs. The required methane most likely migrated upward from deep conventional hydrocarbon reservoirs through the existing fault systems (Dixon and Dietrich, 1988;Dixon et al., 1992;Brent et al., 2005). According to the sediment thicknesses, the bulk of reserves rests under the Beaufort Sea and about one-third are located onshore Osadetz et al., 2005).

3D structural geological model implementation
The structural geological 3D static model covers an area of 10 km × 20 km. With an areal size of 12.7 km × 3 km, the numerical model encompasses the central geological model area, including the Mallik anticline, major F1 Fault, and Mallik wells ( Figure 1). The surfaces of the Iperk/Mackenzie Bay, Mackenzie Bay/Kugmallit, and Kugmallit/Richards Sequence boundaries were generated using six high-resolution 2D seismic reflection profiles (84942, 85251, 85984, 85987, 85988, and 86055), which have been commissioned by Imperial Oil Ltd. in 1984 as part of the seismic survey across the Ivik-Mallik-Taglu area and reprocessed by Collett et al. (1999). Further, three 2D seismic reflection profiles from Brent et al. (2005) were used to identify lithological contacts and fault zones. In the first instance, the reflection profiles were georeferenced to the UTM grid system (Spatial reference: EPSG Projection 32608 -WGS 84/UTM zone 8N). Two-way traveltimes of the sequence boundaries were digitized from the 2D seismic reflection profiles using the Petrel™ software package (Schlumberger, 2012) and converted to elevation depth based on average checkshotand VSP-derived velocities in the Mallik area between 2.25 km/s and 3.0 km/s (Brent et al., 2005). Vertical fault lines of six major regional fault zones were derived from the seismic profiles. Fault strikes were adopted from the two-way traveltime structure map from Brent et al. (2005) for Faults F1-F6. Two additional fault zones at the top of the Mallik anticline structure were interpreted from the 2D seismic profiles in Collett et al. (1999) for Fault F7. The derived set of polylines was then used to calculate the elevation depth of the surfaces and to generate the fault zones using the convergent interpolation algorithm (Schlumberger, 2012).
Available log data from Collett et al. (1999) on the depth of the relevant sequence boundaries from eleven wells in the Ivik-Mallik-Taglu area, all drilled by Imperial Oil Ltd. in the early 1970s, were used for depth adjustment of the processed surfaces (Figure 3). In addition, reflection profiles and wells outside the study area were taken into account to implement stratification inclination trends correctly. For a more detailed description of the model implementation, the reader is kindly referred to Chabab and Kempka (2023b).

Numerical model implementation 2.2.1 Numerical modeling assumptions
A framework of equations of state for equilibrium CH 4 hydrate formation and the reversible processes describing permafrost aggrading/degrading, i.e., water-freezing and ice-melting (Mottaghy and Rath, 2006), has been coupled with the open-source flow and transport simulator (TRANSPORTSE) developed by Kempka (2020) and verified for complex reactive transport problems (Kempka et al., 2022). This coupling was referred to as T plus H (TRANSPORTSE + Hydrate) in our previous studies (Li et al., 2022a;. Specimen analyses  suggest that migration of dissolved CH 4 and GH formation occurred in a semiclosed hydrodynamic system in the vicinity of the Mallik 5L-38 well. The GH intervals at the Mallik 2L-38 were defined as Class-II GH deposits (Uddin et al., 2012;Uddin et al., 2014). According to Moridis and Reagan (2011), Class-II GH deposits are GH intervals capped by impermeable rocks and underlain by aquifers without the presence of a free gas phase. Here, the implementation of impermeable permafrost and caprock sequences located above the  Chabab and Kempka (2023b). For a more detailed view of the employed 2D seismic profiles, the reader is kindly referred to Collett et al. (1999); Brent et al. (2005).
GH-bearing sediments introduces a Class-II GH deposit system in the present simulation.
Hydrocarbon fluids from the deeper sequences are likely to migrate upward through present fault systems (Xia et al., 2022), gas chimneys (McNeil et al., 2011), mud diapirs, mud volcanoes (Zhang et al., 2021), and unconformities (Levell et al., 2010). However, gas chimneys and unconformities, mud diapirs, and mud volcanoes have not yet been reported near the onshore subpermafrost GH reservoirs in the MD. Therefore, the fault systems detected by seismic surveys are regarded as primary conduits for the circulation of fluids beneath the permafrost, facilitating the migration of CH 4 dissolved in these fluids and the formation of GH accumulations (Hillman et al., 2020).
According to the analysis of core samples from the GH-bearing intervals undertaken by Lorenson et al. (1999), Lorenson et al. (2005), the feed gas composition of the GHs derived from the methane-to-hydrocarbon molecular gas ratios is composed of approximately >99.5% thermogenic CH 4 . This supports the model implementation based on the GH-forming gas being solely composed of CH 4 . In addition, the vitrinite reflectance of the GHbearing sediments found in the Mallik 2L-38 well is relatively low with <0.28% (Snowdon, 1999). This indicates that the maturity of the organic matter within the sediment is too low to result in methanogenic activities. Therefore, the feed gas trapped in the GH intervals must have been generated from deeper, thermally matured sequences, probably at depths >5,000 mbgl (Lorenson et al., 1999;Lorenson et al., 2005). This gas migrated upward into the overlying sequences including the Taglu and Richards Sequences, and was then trapped within conventional natural gas reservoirs until its secondary migration via geologic faults supplied the feed gas to shallower geologic units such as the Kugmallit Sequence ( Figure 2).
As demonstrated in our previous study (Li et al., 2022a), thick hydrate-bearing sediment intervals with high S h are commonly associated with deeper conventional petroleum reservoirs, connected to extensional faults. They are closely correlated to the high spatial density of the number of structural elements (i.e., faults and anticlines). Hence, many high-angle dipping faults (F1, F2, F3, and F7, see Figure 4) occurring in the Taglu fault zone (Dixon et al., 2019) have the potential to be preferential pathways for the vertical migration of CH 4 -rich fluids (Chen et al., 2008). However, a quantification of the extent to which these faults control the subsurface temperature regime and CH 4 -rich fluid migration, and the heterogeneous distribution of GH accumulations with high S h has not yet been undertaken. Furthermore, altering the permeabilities of the faults in the numerical model allows for the identification of hydraulically active and closed faults, acting as fluid flow paths or barriers.

Numerical model geometry and parametrization
The input dataset derived from the 3D static geological model covers a volume of 20 km × 10 km × 2 km, and considers seven faults (F1 to F7) and two sedimentary sequence boundaries, including the Iperk-Mackenzie Bay (I-M) and Mackenzie Bay-Kugmallit (M-K) ones. Based on these pre-processed data, a 3D numerical model geometry (Figure 4) of the Mallik site was elaborated for the first time. Subsequently, that model was further applied in the numerical investigation on the geologic controls determining the temporal and spatial formation of sub-permafrost GH accumulations. Figure 4 shows the spatial extent of the employed numerical 3D model for the present simulation study, consisting of three Sequences (Iperk, Mackenzie Bay, and Kugmallit Sequences) and four geologic fault zones (F1, F2, F3, and F7). The grid elements located at the P3-P5 and P2-P8 planes represent open flow boundaries, implemented by constant chemical species concentrations and p-T conditions (Dirichlet boundary conditions, Figure 4B). Moreover, the grid elements located at the P1-P3 plane represent the ground surface and are applied as impermeable boundary with constant p-T conditions (Neumann "no-flow" boundary condition). Figure 4C shows the grid elements of the I-M Sequences, representing the above-mentioned impermeable boundary, acting as impermeable anticline sediment that overlies the GH trapping layer. In addition, constant pressure and regionallydependent heat flux ( Table 2) conditions apply for the grid elements on the P5-P7 plane (Figure 4D,F), except for the grid elements belonging to F1. The grid elements at the bottom of F1 (Figure 4F) serve as an inlet for the CH 4 -rich fluid, which flows into the model domain at a constant rate, and are parametrized with constant temperature conditions, as listed in Table 2. The applied initial temperature of the constant temperature condition is determined from the permafrost formation simulation study, which lasts about 0.6 Ma prior to the subsequent simulation on sub-permafrost GH formation for 1 Ma.
As depicted in Figure 4, the study domain with a spatial extent of 12,720 m × 3,000 m × 1,500 m, is discretized by 148,500 (66 × 30 × 75) grid elements. For optimum computational efficiency, grid elements along the y-direction are discretized by variable sizes of 480 m, 240 m, and 120 m (Figure 4F), while a constant discretization of 100 m and 20 m along the x-and z-directions is applied, respectively.
The subsurface temperature regime in the MD has been affected by low subaerial temperatures since the Middle Pliocene (ca. 3 MaBP). Meanwhile, thin permafrost may have progressively formed as the Arctic subaerial temperature decreased below the freezing point and heat flux dropped back to the regional mean value of 55 mW m -2 from its peak of 80 mW m -2 triggered by the Late Miocene uplift event (Kroeger et al., 2008). As a result, the permafrost bottom accelerated to deepen and eventually pushed the freezing point isotherm downward to the base of the Iperk Sequence since the early Pleistocene sea-level retreat (Figure 5A), whereby the permafrost thickness kept growing until 1 Ma ago. Under the impermeable section provided by the thick permafrost (Chuvilin et al., 2022), super-saturated dissolved CH 4 could be stored in form of GHs and prevented from escaping through upper sequences into the atmosphere. Overall, Figure 5 implies that there may be a time window allowing for the formation of thick permafrost during the past 1.6 Ma to 1.0 Ma, and the subsequent genesis of sub-permafrost GH occurrences since 1.0 MaBP. The time window derived from recent studies Hansen et al., 2013) is consistent with the GH formation timing (Collett, 1993), assuming that GH occurrences found on the Alaska coastal plain of the circum-Beaufort Sea are likely younger than 1.65 Ma.
According to the active source rock depth limit (Pang et al., 2020), the basin's heat flow generally controls the maximum burial depth for source rocks to generate and expel hydrocarbons produced by thermal cracking of kerogen (Behar et al., 1997). Since 1.6 MaBP, the Beaufort MD Basin exhibits an approximate heat flow of 55 mW m -2 ( Figure 5B) in agreement with the moderate-heat-flow (40-60 mW m -2 ) basin suggested by Pang et al. (2020). It implies that the upper limit of hydrocarbon generation, migration and accumulation is below 2,500 mbgl, conforming to the assumption  that the GH-forming thermogenic gases originated from the Taglu Sequence ( Figure 2B). Above 2,500 mbgl, the accumulated hydrocarbons initially migrated from the underlying source rocks, coinciding with our previously validated GH formation mechanism (Li et al., 2022a).
The boundary conditions and hydrothermal properties of the sediments representing the Mallik site were determined by applying an iterative history matching procedure and are summarized in Tables 2, 3. The model employs the parameters reflecting a general average of values reported in the literature. According to well log interpretations and core analyses of the Mallik 2L-38 and 5L-38 wells, the effective porosity is 0.24-0.4 with a mean permeability of 2.9 mD. Fluid pressure is hydrostatic, and the water table is assumed to coincide with the ground surface. The initially relatively homogeneous geophysical property distribution within the GHbearing sediments is altered significantly, as effective composite porosities and permeabilities decreased with the increment in S h .

Results and Discussion
Our simulation results confirm that CH 4 originating from the deep conventional hydrocarbon reservoirs started to flow into the target sandy sediments via Fault F1 in the dissolved state since 1 MaBP. As a hydraulically active preferential pathway, Fault F1 allows the upward migration of CH 4 into the MHSZ located within the Kugmallit Sequence below the Mallik anticline crest. Concurrently, Fault F2 acts as an active hydrological conduit for lateral CH 4 -rich fluid flow, while Fault F3 is as permeable as the sandy sediments of the Kugmallit Sequence. In contrast, the normal (F7) and reverse faults (F4 and F6 in Figure 1) are impermeable hydrological barriers. During the past 1 Ma, the super-saturated amount of dissolved CH 4 was instantly consumed by GH formation, triggered by the decrement in CH 4 solubility as fluid migrated into shallower depths, where p and T decreased. As shown in Figures 4A,C, the Iperk Sequence is penetrated by the top of Fault F1, which supplies super-saturated dissolved-CH 4 formation fluids that facilitate the formation of intra-permafrost GHs (Dallimore and Collett, 1995) until Fault F1 is clogged by the local decrease in porosity. However, the genesis of intra-permafrost GH is beyond the focus of the present study. Consequently, the upper part of Fault F1 intersecting with the Mackenzie Bay Sequence was parametrized as impermeable.
3.1 Simulated permafrost, sub-permafrost GH interval, and subsurface temperature distributions Figure 6A shows the simulated sub-permafrost S h distribution after a simulation time of 1 Ma, implying that the spatial density of the number of structural elements positively correlates with the heterogeneous spatial distribution of sub-permafrost GH accumulations. The densest GH enrichment surrounds Fault F1 below the anticline crest and extends coastward, where subsurface flow discharges into the sea. Further, Figure 6A indicates the locations of the Mallik well, for which the present 3D simulation results are compared against our previous 2D numerical simulation results, as illustrated in Figures 6B-D, 6E. According to the sediment analysis conducted by Medioli et al. (2005) for Mallik 5L-38, GHs occur mostly in the pore space within the thick sand layers at 886-1,108 m depth ( Figure 6F). The indicated thickness of the GH interval at Mallik 5L-38 of 222 m deviates from that observed at Mallik L-38 by approximately 3% (Table 4). Therefore, the short distances (ca. 117-252 m) between the industrial well Mallik L-38 and scientific wells Mallik 2L-38, 3L-38, 4L-38, and 5L-38 (Table 1) are negligible. Consequently, the simulation results acquired from the Mallik L-38 well equivalently represent its four neighboring wells, since their close proximity is represented by one grid element in the numerical model.
For >87% of the Mallik wells summarized in Table 1, the permafrost thickness predicted by the 3D model agrees with the well logs and seismic observations by >99.2% (8 m) ( Table 4), which validates the overall permafrost evolution model implementation. Around the Mallik anticline, the ice-bearing permafrost and GH   Frozen point of pore fluid −2.5°C Taylor et al. (2013) interval thicknesses as well as sub-permafrost temperature profiles simulated by the 3D model are consistent with their respective field observations at the Mallik P-59, L-38, 2L-38, 3L-38, 4L-38, 5L-38, and J-37 wells. This confirms our hypothesis on the Canadian Arctic sub-permafrost GH formation mechanism and the validity of the input parameters listed in Table 2 and 3. As plotted in Figures 6C,D, the S ice and S h profiles simulated by means of the present 3D model are in excellent agreement with our previously presented 2D simulation results for the Mallik wells L-38 and J-37. Deviations between the 2D and 3D simulation results regarding the S h increase with the projection distances from the Mallik P-59 and A-06 wells to the Mallik J-37 well, as indicated in Figures 6B,E. The 2D simulation results on the Mallik P-59 and A-06 wells were acquired from the projected well locations onto the transect of the seismic profile 85987 (Li et al., 2022a), which intersects the Mallik J-37 well. It should be noted that the 2D projected location of the Mallik A-06 well is not located in the seismic and well-logging domains  but at a distance of ca. 500 m (Li et al., 2022a). To date, the published lithology, S h , and subsurface temperature data observed at the Mallik P-59, L-38, J-37, and A-06 wells are limited due to the resolution of the industrial exploration that was targeted at the deep oil reservoirs in the 1970s (Taylor and Judge, 1976). For this reason, the lack of observed peak S h is insufficient to validate the simulated data listed in Table 4. Thus, it merits further efforts to improve the mapping resolution of the geophysically interpreted GH resources as new field data become available, as previously discussed by Bellefleur et al. (2012). Figure 6F shows that there is an evident correlation between the well-logged highly heterogeneous S h with regard to lithology along the borehole profile of the Mallik 5L-38 well ( Table 1). This indicates that coarse-grained layers consisting of permeable sands and pebbles host abundant GHs, whereas fine-grained sediments with nearly no permeability (Katsube et al., 2005), including shales and silts, contain significantly fewer GHs . This relationship implies the conclusion that the S h distribution is lithologically controlled (i.e., by porosity and permeability) according to Jenner et al. (1999); Matsumoto et al. (2005). In contrast to Figure 6F, the simulated 2D and 3D S h intervals at the Mallik L-38 well demonstrate a relatively uniform S h distribution within the Kugmallit Sequence, which is parameterized with homogeneous permeability (Figure 6C). Since homogeneous porosities and pore-size distributions were employed in all simulations, it is suggested that permeability is the local constraint on the pore occupancy limit by GH. However, exhaustive surveys on the distribution of sub-permafrost GH resources have not yet been undertaken at the Mallik site. Furthermore, our simulations are subject to the poorly investigated petrophysical properties regarding the spatial permeability distribution in the reservoir. Reliable data supporting the parameterization of a heterogeneous permeability distribution in the reservoir sequences of interest cannot be derived from these. Table 4 shows that the ice-bearing permafrost thicknesses and peak S h match the field observations with negligible deviations ranging from −1.9% to 0.8%, and 2.2%, respectively. All simulated depths of the MHSZ base are within the tolerance range of the corresponding observations. In addition, the simulated total thicknesses of the GH intervals near the Mallik anticline match their respective observations with minor deviations of only 2.6%-4.7%, which is close to the employed simulation resolution of 1.3% along the vertical direction. However, our model overestimates the thicknesses of the sub-permafrost hydrate-bearing sediment layer at the Mallik A-06 and P-59 wells. This probably indicates that the reservoir quality as determined by porosity, shale-to-sandstone ratio, and intrinsic permeability (Boswell et al., 2011;Boswell et al., 2020) in the far-field of the Mallik anticline, is much poorer than the available near-field data suggests.
Another aspect noted by Frederick and Buffett (2016) is that the thick GH intervals with high S h are easier to delineate, whereas GH intervals with low S h are not likely to be detected by geophysical measurements. For example, GHs have been documented where no Bottom Simulating Reflector (BSR) was observed (Paull et al., 1996). Consequently, a BSR is a strong indicator for the presence of at least low-S h GHs, but the lack of a BSR does not necessarily imply the absence of GHs, as stated by Tréhu et al. (2006). Thus, it is very likely that the industrial exploration of hydrocarbon resources may have underestimated the far-field thicknesses of the identified GH intervals.
The 3D simulation results reveal the highly variable S h of the Mallik GH deposit throughout its vertical and horizontal directions. The largest economically developable GH occurrences are concentrated under the northern flank of the Mallik anticline along Fault F1, as illustrated in Figure 6A. In contrast to the scientific wells, such as Mallik 5L-38, it is not feasible to conduct quantitative analyses but qualitative comparisons of the peak S h at the Mallik A-06, J-37, and P-59 wells, given the limited observations on S h derived from the boreholes. Despite the much less concentrated S h of < 5% predicted by the numerical simulations at the Mallik A-06 well, the far-field remains worthwhile for further investigation in the scope of new drilling campaigns.
The subsurface p-T conditions dominate the stability of the permafrost and the sub-permafrost GH occurrences, while permafrost stability is much more susceptible to global warming than pressure perturbations caused by sea-level fluctuations. Generally, a change in temperature at the ground surface will disrupt subsurface temperature regimes originally equilibrated with the basal heat flux below the permafrost and the GH-bearing sandy sediment. For permafrost at depths of <600 mbgl, the temperature profiles simulated by the 2D model and those logged by DTS demonstrate that the present-day permafrost is warmer by 0.8-1.3 K compared to the 2D numerical predictions, as analyzed and discussed in our previous work (Li et al., 2022a). According to representative near-surface temperatures collected from the exploration wells, the present-day near-surface ground temperature ranges from −8 to −9°C at the Mallik wells (Burn and Kokelj, 2009). Without implementing the previously proposed transitional boundary condition to mimic the climate changes since the Late Holocene (Li et al., 2022a), the present 3D model reproduces the paleo-geothermal conditions of the Arctic permafrost very well.
As shown in Figure 7, the sub-permafrost temperature profiles (depth >600 mbgl) at the Mallik L-38 and J-37 wells simulated by the 3D model are consistent with the DTS observations adopted from the Mallik 3L-38, 4L-38 and 5L-38 wells. In contrast, the temperature at the Mallik J-37 well simulated by the 2D model deviates from the DTS observations by almost 2 K. Although this relative variation may fit the 2D prediction of the regional temperature constraint by Chen et al. (2008), the authors suggested that subsurface temperature variations could reach about 5 K/km along the lateral direction at a depth of 1,100 mbgl. This may be attributed to the uncertainty arising from plotting well data linked to line charts on the temperature cross sections (Chen et al., 2008) or changes in the sub-permafrost geothermal environment under contemporary climate settings and fault activities induced by historical seismicity (Hyndman et al., 2005). On the other hand, the consistency between the 3D-simulated and DTS-logged temperature profiles within the sub-permafrost GH intervals (860 mbgl < depth < 1,100 mbgl) suggests that the Late Holocene climate warming has not yet significantly altered the thermal GH interval environment.

Geologic controls on the subsurface temperature, permafrost, and GH interval distribution
The spatial and temporal evolution of the vertically averaged S h and cumulative GH thicknesses is presented in Figure 8. At first, the upward migrating CH 4 -rich fluid flows into the sandy Kugmallit Sequence within the MHSZ via the preferential pathway Fault F1, which is overlain by the Mallik anticline crest. After a simulation time of 0.1 Ma (Figure 8A), the average S h of the accumulated GHs reaches ca. 21% at the upper part of Fault F1. As GHs continue to form up to a simulation time of 1 Ma, more GHs accumulate along Faults F2 and F3 towards the coast in the MD. This observation implies that the simulated CH 4 -rich fluid flows towards the coast, joining the groundwater discharge to the sea, which contradicts the 2D model results indicating that the CH 4 -rich fluid Note: The well-log-inferred data are adopted from Majorowicz and Smith (1999); Collett et al. (1999); Dallimore et al. (2005). The dashes represent field observation data that could not be derived from field reports Majorowicz and Smith, 1999;Dallimore et al., 2005).

FIGURE 7
Distributed Temperature Sensing (DTS)-logged and simulated temperature profiles. DTS-observations at the Mallik 3L-38, 4L-38, and 5L-38 wells were obtained after 622, 605, and 575 days of their final cementation during the second post-field DTS survey (2003/09/19-21), after Henninges et al. (2005). 2D simulation results were adopted from our previous model (Li et al., 2022a), and the 3D temperature profiles extracted from the present simulations after the sub-permafrost GH formation of 1 Ma. Sequence boundaries are modified after .
migrates inland. This also emphasizes the influence of lithological controls on the GH distribution and S h enrichment. As mentioned earlier, the simulation results reach a very good agreement with the well logs and seismic interpretations  from the Mallik P-59, L-38, J-37, and A-06 wells (Table 4), which exhibit substantially variable S h distributions along both lateral and vertical directions. Figure 8C shows that the GH deposit at the Mallik site is laterally heterogeneously distributed and covers an area of ca. 9 km 2 , where an average vertical S h above 30% persists after a simulation time of 1 Ma. Moreover, Figure 8D indicates that the GH deposit encompasses an area of thick GH-bearing sediments (thickness >250 m). Red shading outlines the reservoir with the thickest GHbearing sediments (thickness >500 m) located below the Mallik anticline, which are penetrated by the Mallik J-37 well. However, the majority of the GH-bearing intervals only contains negligible S h < 1% along the Mallik J-37 well (Figure 6D), as shown in Figure 8C.
One important trend reproduced by the 3D model is that GHs mainly accumulate below the crest of the anticline trap bounded by the Faults F1 and F3 (Figures 8C,D). The reason for this distribution pattern lies in variations in the lithology and geothermal environment, indicating differences in sediment permeability as well as specific heat conductivity and capacity, mainly controlled by the shale and sand contents of the respective sedimentary units. Besides, the GH occurrences with the highest concentrations cover a region of ca. 3 km 2 along the northern flank of the Mallik anticline and east of the connection line between the J-37 to L-38 wells, conditioned by the overlain I-M Sequences which act as hydraulic barriers. Since neither seismic nor borehole data are available to confirm the presence or absence of these highly concentrated GH accumulations along the anticline flank bounded by Fault F1, it highlights the feasibility of the employed 3D model to predict yet undiscovered GH resources.

Conclusion
In the course of the present study, we developed a field-scale static structural geological model of the Mallik anticline based on available well and seismic data. In addition, we employed a thermo-hydro-chemical simulator, previously validated against Mallik field data, to study the geologic controls on the genesis of permafrost and sub-permafrost GH accumulations. The simulation results show that the calculated sub-permafrost temperature profiles, thicknesses of permafrost and hydrate intervals, as well as peak S h are consistent with field observations. Simultaneously, the locally predicted GH distributions match the observed GH occurrences. Therefore, the proposed feed gas (thermogenic CH 4 ) migration and sub-permafrost GH accumulation mechanisms have been validated by the present study, with the following conclusions reached: 1. The complex stratigraphic and structural controls on the heterogeneity of the GH distribution at reservoir scale have been quantitatively confirmed by iterative history matching of the permeability of the underlying faults and lithological units in the Mallik anticline. Since the Late Pleistocene, the major normal fault F1, serving as an active hydraulic conduit, permits the upward migration of hydrate-forming gas in the dissolved state into the MHSZ within the reservoir sequences below the Mallik anticline crest. Furthermore, the normal fault F2 acts as a hydrological conduit for lateral CH 4 -rich fluid flow, while the normal fault F3 is as permeable as the sandy GH-bearing sediments. The normal (F7) and reverse faults (F4, F6) are most likely hydrological barriers.
2. At the Mallik site, the GH occurrences mainly develop in the sandy sediments within the Kugmallit Sequence, whereas the Iperk and Mackenzie Bay Sequences act as seals (barriers for CH 4 transport). The stability of sub-permafrost GHs is mainly determined by the present geothermal regime and hydrodynamic conditions. Concentrated GH accumulations under the Mallik anticline crest are preserved by the favorable thermal conditions around Faults F1, F2 and F3, accompanied by a high spatial density of structural elements. 3. By comparing the simulated results against corresponding welllog data, the presented numerical framework confirms the solid quality of the established static geological model elaborated on the basis of published interpretations of the seismic reflection profiles. The introduced simulation framework can be applied to support the interpretation of such geophysical measurements.
The presented results allow for the assessment and projection of potential GH enrichments at the Mallik site and provide a new perspective on present Arctic sub-permafrost hydrocarbon systems. Moreover, they contribute to the understanding of sub-permafrost groundwater flow, demonstrating that the preferred flow direction of the CH 4 -rich fluid is coastward to join the groundwater discharge in the Canadian MD. Our model also shows excellent potentials for investigating the genesis of similar integrated natural systems and the evolution of permafrost under various climatological events in the North American coastal region of the Beaufort Sea. The validated simulation framework can be applied to other worldwide GH deposits to contribute to the assessment of global GH resources.
Frontiers in Earth Science 13 frontiersin.org

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
Conceptualization, TK and ZL; methodology, ZL, EC, and TK; software, ZL, EC, and TK; verification and validation, ZL; formal analysis, ZL; investigation, ZL; writing-original draft preparation, ZL and EC; writing-review and editing, ZL, TK, EC, ES and JS; visualization, ZL and EC; supervision, TK; All authors have read and agreed to the published version of the manuscript.

Funding
The first author thanks for the financial support provided by China Scholarship Council (Grant No. 201806930024). The authors acknowledge the financial support granted for the article processing charge by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) within the funding program "Open Access Publikationskosten"-Project Number 491075472.