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 (CH₄) at the Mallik site. The present study aims to validate the proposed Arctic sub-permafrost gas hydrate formation mechanism, implying that CH₄-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 CH₄-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.

1 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₄), trapped within the water molecules forming cage-shaped 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., 1999; 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₄, 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₄-rich fluid is migrating upward out of the deeper hydraulic overpressure zone (Hu et al., 2018) through geologic gas-source faults to form CH₄ hydrates in the shallower Kugmallit Sequence (Figure 1) over geological times. When the CH₄-rich fluid enters the Methane Hydrate Stability Zone (MHSZ), the dissolved CH₄ 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.

FIGURE 1

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

FIGURE 2

FIGURE 2. (A) Schematic and (B) chronology of the sub-permafrost GH-bearing system at the Mallik site with an active source rock generating thermogenic CH₄, its geological elements (trap, seal, reservoir, and source intervals), and processes (migration, accumulation, and preservation), modified after Li et al. (2022a), according to Chen et al. (2008); Kroeger et al. (2008). (Seq. = sequence, not to scale).

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₄-loaded fluids rather than highly-saturated with GHs due to the lack of a geologic trap and the absence of CH₄ super-saturated formation fluids (Huang, 2009).

TABLE 1

TABLE 1. Overview of scientific and industrial exploration wells at the Mallik site.

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₄ 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 permafrost-associated GH formation via dissolved CH₄. 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 sub-permafrost 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₄-rich formation fluids. The proposed mechanism for the formation of highly saturated sub-permafrost GH via dissolved CH₄-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 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), Bellefleur et al. (2012). 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₄-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² of the MD in the present study. In addition to well logs (Table 1), published seismic interpretations (Dallimore et al., 1999; 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.

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 (Dallimore et al., 1999). 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 (Dallimore et al., 1999), followed by the development of the Mallik 3L-38, 4L-38, and 5L-38 wells in 2002, as listed in Table 1.

2.1.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 geological period from the Paleozoic to the Holocene. Westward-thickening 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 (Collett et al., 1999). 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 delta-front 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 (Dallimore and Collett, 1999; Osadetz et al., 2005).

2.1.2 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 checkshot- and 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).

FIGURE 3

FIGURE 3. Domain and geometry of the implemented 3D Mallik static geological model with major sequence boundaries. (A) Elevation depth of Iperk/Mackenzie Bay and (B) Mackenzie Bay/Kugmallit sequence boundaries, (C) implemented regional faults zones of the Ivik–Mallik–Taglu area, and (D) sequence boundaries digitized from two-way traveltime of employed seismic profiles. Vertical exaggeration by 5 times, spatial reference: WGS 84/UTM zone 8N, modified after 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).

2.2 Numerical model implementation

2.2.1 Numerical modeling assumptions

A framework of equations of state for equilibrium CH₄ 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_plusH (TRANSPORTSE + Hydrate) in our previous studies (Li et al., 2022a; b).

Specimen analyses (Wright et al., 2005) suggest that migration of dissolved CH₄ 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 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 sub-permafrost 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₄ 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₄. This supports the model implementation based on the GH-forming gas being solely composed of CH₄. In addition, the vitrinite reflectance of the GH-bearing 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₄-rich fluids (Chen et al., 2008). However, a quantification of the extent to which these faults control the subsurface temperature regime and CH₄-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.

FIGURE 4

FIGURE 4. (A,B) Overview of the model geometry, including the fault systems (F1, F2, F3, and F7) as well as Iperk, Mackenzie Bay, and Kugmallit Sequences. (C,D,E) Model geometry dimensions along the x-z, x-y, and y-z planes, respectively. (F) Lateral grid discretization employed in the present simulations.

2.2.2 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 regionally-dependent 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₄-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.

TABLE 2

TABLE 2. Initial and boundary conditions applied in the present simulations.

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₄ 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 (Dallimore et al., 1999; 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.

FIGURE 5

FIGURE 5. (A) Global mean sea-level estimations relative to the Late Holocene sea-level since 3 MaBP (Hansen et al., 2013). (B) Time series of Arctic subaerial temperature (Arctic T) and basal heat flow in the past 3 Ma, adopted from Kroeger et al. (2008).

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 (Dallimore et al., 1999) 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 GH-bearing sediments is altered significantly, as effective composite porosities and permeabilities decreased with the increment in S_h.

TABLE 3

TABLE 3. List of hydrothermal properties of the permafrost and GH-bearing sediments, and other parameters used in the presented simulations.

3 Results and Discussion

Our simulation results confirm that CH₄ 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₄ 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₄-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₄ was instantly consumed by GH formation, triggered by the decrement in CH₄ 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₄ 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.

FIGURE 6

FIGURE 6. (A) Simulated field-scale sub-permafrost GH accumulation distribution at a simulation time of 1 Ma (see Figure 4 for model geometry). Simulated ice-bearing permafrost saturation (S_ice) and sub-permafrost gas hydrate saturation (S_h) profiles obtained for the respective locations of the (B) Mallik P-59, (C) L-38, (D) J-37, and (E) A-06 wells. (F) Lithology and S_h logs from the Mallik 5L-38 well after Bauer et al. (2005), presented with permission from Natural Resources Canada under the Open Government License—Canada version 2.0.

TABLE 4

TABLE 4. Comparison of well logs (observed) with simulated permafrost and GH intervals.

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 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 (Collett et al., 1999) 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 (Medioli et al., 2005). 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.

FIGURE 7

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 Dallimore et al. (1999).

3.2 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₄-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₄-rich fluid flows towards the coast, joining the groundwater discharge to the sea, which contradicts the 2D model results indicating that the CH₄-rich fluid 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 (Collett et al., 1999) 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 8

FIGURE 8. Plan view of temporal series of (A,C) average GH saturations (S_h), and (B,D) total GH thicknesses within the GH-bearing intervals following the sub-permafrost GH formation after (A,B) 0.1 Ma, and (C,D) 1 Ma.

Figure 8C shows that the GH deposit at the Mallik site is laterally heterogeneously distributed and covers an area of ca. 9 km², 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 GH-bearing 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² 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.

4 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₄) 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₄-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₄ 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 well-log 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₄-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 world-wide GH deposits to contribute to the assessment of global GH resources.

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.

Acknowledgments

The first author would like to thank Benjamin Nakaten (GFZ German Research Centre for Geosciences) for providing the model pre-processing tool Geomodelator, used to implement the grid of the numerical model in the present study. The authors acknowledge the constructive comments of the reviewers, which supported improving the quality of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ashford, D., Numasawa, M., Martin, C., Yamamoto, K., Dallimore, S., Wright, J., et al. (2012). “Overview of engineering and operations activities conducted as part of the JOGMEC/NRCan/Aurora Mallik 2007–2008 gas hydrate production research well program, part b: 2008 field program,” in Scientific results from the JOGMEC/NRCan/Aurora Mallik 2007-2008 gas hydrate production research well program, Mackenzie Delta, Northwest Territories, Canada (Geological Survey of Canada), 53–65. Bulletin 601. doi:10.4095/292082

CrossRef Full Text | Google Scholar

Bauer, K., Pratt, R. G., Weber, M. H., Ryberg, T., Haberland, C., and Shimizu, S. (2005). “Mallik 2002 cross-well seismic experiment: Project design, data acquisition, and modelling studies,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 14. Bulletin 585. doi:10.4095/220866

CrossRef Full Text | Google Scholar

Behar, F., Vandenbroucke, M., Tang, Y., Marquis, F., and Espitalie, J. (1997). Thermal cracking of kerogen in open and closed systems: Determination of kinetic parameters and stoichiometric coefficients for oil and gas generation. Org. Geochem. 26, 321–339. doi:10.1016/S0146-6380(97)00014-4

CrossRef Full Text | Google Scholar

Bellefleur, G., Riedel, M., and Brent, T. (2006). Seismic characterization and continuity analysis of gas-hydrate horizons near Mallik research wells, Mackenzie Delta, Canada. Lead. Edge 25, 599–604. doi:10.1190/1.2202663

CrossRef Full Text | Google Scholar

Bellefleur, G., Riedel, M., Brent, T., Wright, F., and Dallimore, S. R. (2007). Implication of seismic attenuation for gas hydrate resource characterization, Mallik, Mackenzie Delta, Canada. J. Geophys. Res. Solid Earth 112, B10311. doi:10.1029/2007JB004976

CrossRef Full Text | Google Scholar

Bellefleur, G., Riedel, M., Huang, J., Saeki, B., Milkereit, T., Ramachandran, K., et al. (2012). Scientific results from the JOGMEC/NRCan/Aurora Mallik 2007–2008 gas hydrate production research well program, Mackenzie Delta, Northwest Territories, Canada. Ottawa: Geological Survey of Canada, Bulletin 601. 601, 107–123. doi:10.4095/292086

CrossRef Full Text | Google Scholar

Boswell, R., Rose, K., Collett, T. S., Lee, M., Winters, W., Lewis, K. A., et al. (2011). Geologic controls on gas hydrate occurrence in the Mount Elbert prospect, Alaska North Slope. Mar. Geophys. Res. 28, 589–607. doi:10.1016/j.marpetgeo.2009.12.004

CrossRef Full Text | Google Scholar

Boswell, R., Hancock, S., Yamamoto, K., Collett, T., Pratap, M., and Lee, S.-R. (2020). “6 – natural gas hydrates: Status of potential as an energy resource,” in Future energy. Editor T. M. Letcher 3rd edn. (Elsevier), 111–131. doi:10.1016/B978-0-08-102886-5.00006-2

CrossRef Full Text | Google Scholar

Brent, T. A., Riedel, M., Caddel, M., Clement, M., Collett, T. S., and Dallimore, S. R. (2005). “Initial geophysical and geological assessment of an industry 3D seismic survey covering the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate research well,” ARRAY(0x55aa16400708 in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada. Editors S. R. Dallimore, and T. S. Collett (Ottawa: Geological Survey of Canada), 14. Bulletin 585. doi:10.4095/220886

CrossRef Full Text | Google Scholar

Burn, C. R., and Kokelj, S. V. (2009). The environment and permafrost of the Mackenzie Delta area. Permafr. Periglac. 20, 83–105. doi:10.1002/ppp.655

CrossRef Full Text | Google Scholar

Chabab, E., and Kempka, T. (2023a). 3D structural geological model of the Mallik Anticline, Canadian Arctic. Potsdam: GFZ Data Services. [Dataset]. doi:10.5880/GFZ.3.4.2023.001

CrossRef Full Text | Google Scholar

Chabab, E., and Kempka, T. (2023b). 3D structural geological model of the Mallik Anticline, Canadian Arctic (report), Tech. rep. Scientific Technical Report STR - Data; 23/01. doi:10.48440/gfz.b103-23010

CrossRef Full Text | Google Scholar

Chen, Z., Osadetz, K. G., Issler, D. R., and Grasby, S. E. (2008). Hydrocarbon migration detected by regional temperature field variations, Beaufort-Mackenzie Basin, Canada. AAPG Bull. 92, 1639–1653. doi:10.1306/07300808011

CrossRef Full Text | Google Scholar

Chen, Z., Dietrich, J., Lane, L. S., Li, M., and Dixon, J. (2021). Beaufort–Mackenzie tectono-sedimentary element. Geol. Soc. Lond. Mem. 57. doi:10.1144/M57-2016-3

CrossRef Full Text | Google Scholar

Chuvilin, E., Bukhanov, B., Yurchenko, A., Davletshina, D., Shakhova, N., Spivak, E., et al. (2022). In-situ temperatures and thermal properties of the East Siberian Arctic shelf sediments: Key input for understanding the dynamics of subsea permafrost. Mar. Geophys. Res. 138, 105550. doi:10.1016/j.marpetgeo.2022.105550

CrossRef Full Text | Google Scholar

Collett, T. S., Lee, M. W., Dallimore, S., and Agena, W. F. (1999). “Seismic- and well-log-inferred gas hydrate accumulations on Richards Island,” in Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 544, 357–376. Bulletin 544. doi:10.4095/210782

CrossRef Full Text | Google Scholar

Collett, T. S. (1993). Natural gas hydrates of the prudhoe Bay and kuparuk river area, north slope, Alaska. AAPG Bull. 77, 793–812. doi:10.1306/BDFF8D62-1718-11D7-8645000102C1865D

CrossRef Full Text | Google Scholar

Dallimore, S. R., and Collett, T. S. (1995). Intrapermafrost gas hydrates from a deep core hole in the Mackenzie Delta, Northwest Territories, Canada. Geology 23, 527. doi:10.1130/0091-7613(1995)023⟨0527:IGHFAD⟩2.3.CO;2

CrossRef Full Text | Google Scholar

Dallimore, S. R., and Collett, T. S. (1999). “Regional gas hydrate occurrences, permafrost conditions, and Cenozoic geology, Mackenzie Delta area,” in Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 544, 31–43. Bulletin 544. doi:10.4095/210746

CrossRef Full Text | Google Scholar

Dallimore, S. R., Uchida, T., and Collett, T. S. (1999). Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada. Tech. rep. Canada. doi:10.4095/210723

CrossRef Full Text | Google Scholar

Dallimore, S. R., and Collett, T. S. (2005). Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada. Tech. Rep., 585. Canada. doi:10.4095/220702

CrossRef Full Text | Google Scholar

Dallimore, S. R., Yamamoto, K., Wright, J. F., and Bellefleur, G. (2012). Scientific results from the JOGMEC/NRCan/Aurora Mallik 2007–2008 gas hydrate production research well program, Mackenzie Delta, northwest territories, Canada, 601. Geological Survey of Canada, 291. Bulletin 601. doi:10.4095/291751

CrossRef Full Text | Google Scholar

Dixon, J., and Dietrich, J. R. (1988). “The nature of depositional and seismic sequence boundaries in Cretaceous-Tertiary strata of the Beaufort-Mackenzie Basin,” in Sequences, stratigraphy, sedimentology: Surface and subsurface; proceedings of A Canadian society of petroleum geologists technical meeting (Canadian Society of Petroleum Geologists), 63–72. Memoir no. 15.

Google Scholar

Dixon, J., Dietrich, J. R., and McNeil, D. H. (1992). Upper cretaceous to Pleistocene sequence stratigraphy of the beaufort-mackenzie and banks island areas, northwest Canada. Tech. rep. Canada. doi:10.4095/133237

CrossRef Full Text | Google Scholar

Dixon, J., Lane, L., Dietrich, J., McNeil, D., and Chen, Z. (2019). “Chapter 17 – geological history of the late cretaceous to cenozoic beaufort-mackenzie basin, arctic Canada,” in The sedimentary basins of the United States and Canada. Editor A. D. Miall 2nd edn. (Elsevier), 695–717. doi:10.1016/B978-0-444-63895-3.00017-6

CrossRef Full Text | Google Scholar

Forbes, D., Craymer, M., James, T., and Whalen, D. (2022). Subsidence drives habitat loss in a large permafrost delta, Mackenzie River outlet to the Beaufort Sea, Western Arctic Canada. Can. J. Earth Sci. 59, 914–934. doi:10.1139/cjes-2021-0127

CrossRef Full Text | Google Scholar

Frederick, J. M., and Buffett, B. A. (2014). Taliks in relict submarine permafrost and methane hydrate deposits: Pathways for gas escape under present and future conditions. J. Geophys. Res. Solid Earth 119, 106–122. doi:10.1002/2013JF002987

CrossRef Full Text | Google Scholar

Frederick, J. M., and Buffett, B. A. (2015). Effects of submarine groundwater discharge on the present-day extent of relict submarine permafrost and gas hydrate stability on the Beaufort Sea continental shelf. J. Geophys. Res. Solid Earth 120, 417–432. doi:10.1002/2014JF003349

CrossRef Full Text | Google Scholar

Frederick, J. M., and Buffett, B. A. (2016). Submarine groundwater discharge as a possible formation mechanism for permafrost-associated gas hydrate on the circum-Arctic continental shelf. J. Geophys. Res. Solid Earth 121, 1383–1404. doi:10.1002/2015JB012627

CrossRef Full Text | Google Scholar

Hansen, J., Sato, M., Russell, G., and Kharecha, P. (2013). Climate sensitivity, sea level and atmospheric carbon dioxide. Philos. Trans. R. Soc. A 371, 20120294. doi:10.1098/rsta.2012.0294

CrossRef Full Text | Google Scholar

Henninges, J., Schrötter, J., Erbas, K., and Huenges, E. (2005). “Temperature field of the Mallik gas hydrate occurrence – implications on phase changes and thermal properties,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 14. Bulletin 585. doi:10.4095/220890

CrossRef Full Text | Google Scholar

Henninges, J. (2005). Thermal properties of gas-hydrate-bearing sediments and effects of phase transitions on the transport of heat deduced from temperature logging at Mallik, NWT, Canada. Doctoral thesis. Berlin: Technische Universität Berlin, ehemalige Fakultät VI - Bauingenieurwesen und Angewandte Geowissenschaften. doi:10.14279/depositonce-1204

CrossRef Full Text | Google Scholar

Hillman, J. I. T., Crutchley, G. J., and Kroeger, K. F. (2020). Investigating the role of faults in fluid migration and gas hydrate formation along the southern Hikurangi Margin, New Zealand. Mar. Geophys. Res. 41, 8. doi:10.1007/s11001-020-09400-2

CrossRef Full Text | Google Scholar

Houseknecht, D. W., Bird, K. J., Garrity, C. P., and Survey, U. S. G. (2020). “Geology and assessment of undiscovered oil and gas resources of the Amerasia Basin Province, 2008,” in The 2008 circum-arctic resource appraisal. Editors T. E. Moore, and D. L. Gautier (U.S. Geological Survey Professional), 33. Paper 1824. doi:10.3133/pp1824BB

CrossRef Full Text | Google Scholar

Hu, K., Dietrich, J., Chen, Z., and Hannigan, P. (2018). Petrophysical analyses of hydrocarbon reservoirs and overpressure zones in Tertiary deep-marine strata in the northern Beaufort-Mackenzie Basin, Arctic Canada. Tech. rep. doi:10.4095/306960

CrossRef Full Text | Google Scholar

Huang, J.-W. (2009). Seismic imaging of gas hydrate reservoir heterogeneities. Doctoral thesis. Canada, Toronto: University of Toronto.

Google Scholar

Hyndman, R., Cassidy, J., Adams, J., Rogers, G., and Mazzotti, S. (2005). Earthquakes and seismic hazard in the Yukon-Beaufort-Mackenzie. CSEG Rec. 5, 32–66.

Google Scholar

Jenner, K., Dallimore, S., Clark, I., Paré, D., and Medioli, B. (1999). “Sedimentology of gas hydrate host strata from the JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well,” in Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada,), 544, 57–68. Bulletin 544. doi:10.4095/210748

CrossRef Full Text | Google Scholar

Judge, A. S., Taylor, A. E., Burgess, M., and Allen, V. S. (1981). Canadian geothermal data collection: Northern wells, 1978-80. Tech. rep. doi:10.4095/8275

CrossRef Full Text | Google Scholar

Kashchiev, D., and Firoozabadi, A. (2002). Driving force for crystallization of gas hydrates. J. Cryst. Growth 241, 220–230. doi:10.1016/S0022-0248(02)01134-X

CrossRef Full Text | Google Scholar

Katsube, T. J., Dallimore, S. R., Jonasson, I. R., Connell-Madore, S., Medioli, B. E., Uchida, T., et al. (2005). “Petrophysical characteristics of gas-hydrate-bearing and gas-hydrate-free formations in the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 14. Bulletin 585. doi:10.4095/220729

CrossRef Full Text | Google Scholar

Kempka, T., Steding, S., and Kühn, M. (2022). Verification of TRANSPORT Simulation Environment coupling with PHREEQC for reactive transport modelling. Adv. Geosci. 58, 19–29. doi:10.5194/adgeo-58-19-2022

CrossRef Full Text | Google Scholar

Kempka, T. (2020). Verification of a Python-based TRANsport Simulation Environment for density-driven fluid flow and coupled transport of heat and chemical species. Adv. Geosci. 54, 67–77. doi:10.5194/adgeo-54-67-2020

CrossRef Full Text | Google Scholar

Kroeger, K. F., Ondrak, R., di Primio, R., and Horsfield, B. (2008). A three-dimensional insight into the Mackenzie Basin (Canada): Implications for the thermal history and hydrocarbon generation potential of Tertiary deltaic sequences. AAPG Bull. 92, 225–247. doi:10.1306/10110707027

CrossRef Full Text | Google Scholar

Kroeger, K. F., di Primio, R., and Horsfield, B. (2009). Hydrocarbon flow modeling in complex structures (Mackenzie Basin, Canada). AAPG Bull. 93, 1209–1234. doi:10.1306/05220908141

CrossRef Full Text | Google Scholar

Lane, L. S. (2002). Tectonic evolution of the Canadian Beaufort Sea–Mackenzie delta region: A brief review. CSEG Rec. 27, 49–56.

Google Scholar

Levell, B., Argent, J., Doré, A. G., and Fraser, S. (2010). Passive margins: Overview. Geol. Soc. Lond. Pet. Geol. Conf. Ser. 7, 823–830. doi:10.1144/0070823

CrossRef Full Text | Google Scholar

Li, Z., Spangenberg, E., Schicks, J. M., and Kempka, T. (2022a). Numerical simulation of coastal sub-permafrost gas hydrate formation in the Mackenzie Delta, Canadian Arctic. Energies 15, 4986. doi:10.3390/en15144986

CrossRef Full Text | Google Scholar

Li, Z., Spangenberg, E., Schicks, J. M., and Kempka, T. (2022b). Numerical simulation of hydrate formation in the LArge-Scale Reservoir Simulator (LARS). Energies 15, 1974. doi:10.3390/en15061974

CrossRef Full Text | Google Scholar

Lorenson, T. D., Whiticar, M. J., Waseda, A., Dallimore, S. R., and Collett, T. S. (1999). “Gas composition and isotopic geochemistry of cuttings, core, and gas hydrate from the JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well,” in Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 143–164. Bulletin 544.

CrossRef Full Text | Google Scholar

Lorenson, T. D., Whiticar, M. J., Collett, T. S., Dallimore, S. R., and Dougherty, J. A. (2005). “Complete gas composition and isotopic geochemistry from the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well: Cuttings, core, gas hydrate, and production testing results,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 19. Bulletin 585. doi:10.4095/220744

CrossRef Full Text | Google Scholar

Majorowicz, J., and Smith, S. (1999). “Review of ground temperatures in the Mallik field area: A constraint to the methane hydrate stability,” in Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 45–56. Bulletin 544. doi:10.4095/210747

CrossRef Full Text | Google Scholar

Majorowicz, J., Osadetz, K., and Safanda, J. (2012a). Gas hydrate formation and dissipation histories in the northern margin of Canada: Beaufort-Mackenzie and the Sverdrup Basins. J. Geol. Res. 2012, 1–17. doi:10.1155/2012/879393

CrossRef Full Text | Google Scholar

Majorowicz, J., Safanda, J., and Osadetz, K. (2012b). Inferred gas hydrate and permafrost stability history models linked to climate change in the Beaufort-Mackenzie Basin, Arctic Canada. Clim. Past. 8, 667–682. doi:10.5194/cp-8-667-2012

CrossRef Full Text | Google Scholar

Majorowicz, J., Osadetz, K., and Safanda, J. (2015). Models of talik, permafrost and gas hydrate histories–Beaufort Mackenzie Basin, Canada. Energies 8, 6738–6764. doi:10.3390/en8076738

CrossRef Full Text | Google Scholar

Matsumoto, R., Tomaru, H., Chen, Y. F., Lu, H., and Clark, I. D. (2005). “Geochemistry of the interstitial waters of the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 14. Bulletin 585. doi:10.4095/220764

CrossRef Full Text | Google Scholar

McNeil, D. H., Dietrich, J. R., Issler, D. R., Grasby, S. E., Dixon, J., and Stasiuk, L. D. (2011). “A new method for recognizing subsurface hydrocarbon seepage and migration using altered foraminifera from a gas chimney in the beaufort-mackenzie basin,” in Shale tectonics (American Association of Petroleum Geologists). doi:10.1306/13231315M933425

CrossRef Full Text | Google Scholar

Medioli, B., Wilson, N., Dallimore, S., Paré, D., Brennan-Alpert, P., and Oda, H. (2005). “Sedimentology of the cored interval, JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 21. Bulletin 585. doi:10.4095/220726

CrossRef Full Text | Google Scholar

Miller, R. D., Hunter, J. A., Doll, W. E., Carr, B. J., and Collett, T. S. (2005). “High-resolution seismic imaging of the gas hydrate stability zone at the Mallik L-38 research site,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 14. Bulletin 585. doi:10.4095/220888

CrossRef Full Text | Google Scholar

Moridis, G. J., and Reagan, M. T. (2011). Estimating the upper limit of gas production from Class 2 hydrate accumulations in the permafrost: 1. Concepts, system description, and the production base case. J. Pet. Sci. Eng. 76, 194–204. doi:10.1016/j.petrol.2010.11.023

CrossRef Full Text | Google Scholar

Moridis, G. J., Kowalsky, M. B., and Pruess, K. (2005). HydrateResSim users manual: A numerical simulator for modeling the behavior of hydrates in geologic media. Tech. rep. Berkeley, CA (United States): Lawrence Berkeley National Lab.

Google Scholar

Mottaghy, D., and Rath, V. (2006). Latent heat effects in subsurface heat transport modelling and their impact on palaeotemperature reconstructions. Geophys. J. Int. 164, 236–245. doi:10.1111/j.1365-246X.2005.02843.x

CrossRef Full Text | Google Scholar

Osadetz, K. G., and Chen, Z. (2010). A re-evaluation of Beaufort Sea-Mackenzie Delta Basin gas hydrate resource potential: Petroleum system approaches to non-conventional gas resource appraisal and geologically-sourced methane flux. Bull. Can. Pet. Geol. 58, 56–71. doi:10.2113/gscpgbull.58.1.56

CrossRef Full Text | Google Scholar

Osadetz, K. G., Morrel, G. R., Dixon, J., Dietrich, J., R., Snowdon, L. R., Dallimore, S. R., et al. (2005). “Beaufort-mackenzie basin: A review of conventional and nonconventional (gas hydrate) petroleum reserves and undiscovered resources,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 19. Bulletin 585. doi:10.4095/220353

CrossRef Full Text | Google Scholar

Overduin, P. P., Schneider von Deimling, T., Miesner, F., Grigoriev, M. N., Ruppel, C., Vasiliev, A., et al. (2019). Submarine permafrost map in the Arctic modeled using 1-D transient heat flux (SuPerMAP). J. Geophys. Res. Oceans 124, 3490–3507. doi:10.1029/2018JC014675

CrossRef Full Text | Google Scholar

Pang, X., Jia, C., Zhang, K., Li, M., Wang, Y., Peng, J., et al. (2020). The dead line for oil and gas and implication for fossil resource prediction. Earth Syst. Sci. Data 12, 577–590. doi:10.5194/essd-12-577-2020

CrossRef Full Text | Google Scholar

Paull, C., Matsumoto, R., and Wallace, P. (1996). Proceedings, initial reports, ocean drilling program, leg 164, gas hydrate sampling on the blake ridge and carolina rise. College Station: Texas A and M University. ODP.

Google Scholar

Riedel, M., Bellefleur, G., Mair, S., Brent, T. A., and Dallimore, S. R. (2009). Acoustic impedance inversion and seismic reflection continuity analysis for delineating gas hydrate resources near the Mallik research sites, Mackenzie Delta, Northwest Territories, Canada. Geophysics 74, B125–B137. doi:10.1190/1.3159612

CrossRef Full Text | Google Scholar

Ripmeester, J., Lu, H., Moudrakovski, I., Dutrisac, R., Wilson, L., Wright, F., et al. (2005). “Structure and composition of gas hydrate in sediment recovered from the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well, determined by X-ray diffraction and Raman and solid-state nuclear magnetic resonance spectroscopy,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 6. Bulletin 585. doi:10.4095/220808

CrossRef Full Text | Google Scholar

Schlumberger (2012). Petrel seismic-to-evaluation software. version 2011.2.7. [Dataset].

Google Scholar

Sippel, J., Scheck-Wenderoth, M., Lewerenz, B., and Kroeger, K. F. (2013). A crust-scale 3D structural model of the Beaufort-Mackenzie Basin (Arctic Canada). Tectonophysics 591, 30–51. Basin Dynamics. doi:10.1016/j.tecto.2012.10.030

CrossRef Full Text | Google Scholar

Sippel, J., Scheck-Wenderoth, M., Lewerenz, B., and Klitzke, P. (2015). Deep vs. shallow controlling factors of the crustal thermal field – insights from 3D modelling of the Beaufort-Mackenzie Basin (Arctic Canada). Basin Res. 27, 102–123. doi:10.1111/bre.12075

CrossRef Full Text | Google Scholar

Sloan, E. D., and Koh, C. A. (2007). Clathrate hydrates of natural gases. 3rd edn. Boca Raton, FL, USA: CRC Press. doi:10.1201/9781420008494

CrossRef Full Text | Google Scholar

Snowdon, L. R. (1999). “Methane and carbon dioxide gas-generation kinetics, JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well,” in Scientific results from JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 544, 125–141. Bulletin 544. doi:10.4095/210772

CrossRef Full Text | Google Scholar

Taylor, A. E., and Judge, A. S. (1976). Canadian geothermal data collection: Northern wells, 1975. Tech. rep. doi:10.4095/8415

CrossRef Full Text | Google Scholar

Taylor, A. E., Dallimore, S. R., Hill, P. R., Issler, D. R., Blasco, S., and Wright, F. (2013). Numerical model of the geothermal regime on the Beaufort Shelf, arctic Canada since the last interglacial. J. Geophys. Res. Solid Earth 118, 2365–2379. doi:10.1002/2013JF002859

CrossRef Full Text | Google Scholar

Thornton, J. M., Mariethoz, G., and Brunner, P. (2018). A 3D geological model of a structurally complex Alpine region as a basis for interdisciplinary research. Sci. Data 5, 180238. doi:10.1038/sdata.2018.238

PubMed Abstract | CrossRef Full Text | Google Scholar

Tréhu, A. M., Ruppel, C., Holland, M., Dickens, G. R., Torres, M. E., Collett, T. S., et al. (2006). Gas hydrates in marine sediments: Lessons from scientific ocean drilling. Oceanography 19, 124–142. doi:10.5670/oceanog.2006.11

CrossRef Full Text | Google Scholar

Uddin, M., Wright, J. F., Dallimore, S. R., and Coombe, D. (2012). “Gas hydrate production from the Mallik reservoir: Numerical history matching and long-term production forecasting,” in Scientific results from the JOGMEC/NRCan/Aurora Mallik 2007-2008 gas hydrate production research well program, Mackenzie Delta, Northwest Territories, Canada (Geological Survey of Canada), 601, 261–289. Bulletin 601. doi:10.4095/292094

CrossRef Full Text | Google Scholar

Uddin, M., Wright, F., Dallimore, S., and Coombe, D. (2014). Gas hydrate dissociations in Mallik hydrate bearing zones a, b, and c by depressurization: Effect of salinity and hydration number in hydrate dissociation. J. Nat. Gas. Eng. 21, 40–63. doi:10.1016/j.jngse.2014.07.027

CrossRef Full Text | Google Scholar

Waite, W. F., Santamarina, J. C., Cortes, D. D., Dugan, B., Espinoza, D. N., Germaine, J., et al. (2009). Physical properties of hydrate-bearing sediments. Rev. Geophys. 47, RG4003. doi:10.1029/2008RG000279

CrossRef Full Text | Google Scholar

Woo, M.-k. (2012). Permafrost hydrology. Berlin, Heidelberg: Springer. doi:10.1007/978-3-642-23462-0

CrossRef Full Text | Google Scholar

Wright, J. F., Dallimore, S. R., Nixon, F. M., and Duchesne, C. (2005). “In situ stability of gas hydrate in reservoir sediments of the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well,” in Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie delta, northwest territories, Canada (Geological Survey of Canada), 585, 11. Bulletin 585. doi:10.4095/220742

CrossRef Full Text | Google Scholar

Xia, Y., Yang, J., Chen, Y., Lu, S., Wang, M., Deng, S., et al. (2022). A review of the global polygonal faults: Are they playing a big role in fluid migration? Front. Earth Sci. 9. doi:10.3389/feart.2021.786915

CrossRef Full Text | Google Scholar

Zhang, W., Liang, J., Liang, Q., Wei, J., Wan, Z., Feng, J., et al. (2021). Gas hydrate accumulation and occurrence associated with cold seep systems in the northern south China sea: An overview. Geofluids 2021, 1–24. doi:10.1155/2021/5571150

CrossRef Full Text | Google Scholar

Nomenclature

Abbreviations

GH gas hydrate

mbgl meters below the ground level

Ma million years

MaBP million years before present

MD Mackenzie Delta

MHSZ methane hydrate stability zone

p pressure

S_h gas hydrate saturation

S_ice ice saturation

T temperature

κ_x horizontal permeability along the x-axis

κ_y horizontal permeability along the y-axis

κ_y vertical permeability