Coupled Basin and Hydro-Mechanical Modeling of Gas Chimney Formation: The SW Barents Sea

Abstract

Gas chimneys are one of the most intriguing manifestations of the focused fluid flows in sedimentary basins. To predict natural and human-induced fluid leakage, it is essential to understand the mechanism of how fluid flow localizes into conductive chimneys and the chimney dynamics. This work predicts conditions and parameters for chimney formation in two fields in the SW Barents Sea, the Tornerose field and the Snøhvit field in the Hammerfest Basin. The work is based on two types of models, basin modeling and hydro-mechanical modeling of chimney formation. Multi-layer basin models were used to produce the initial conditions for the hydro-mechanical modeling of the relatively fast chimneys propagation process. Using hydro-mechanical models, we determined the thermal, structural, and petrophysical features of the gas chimney formation for the Tornerose field and the Snøhvit field. Our hydro-mechanical model treats the propagation of chimneys through lithological boundaries with strong contrasts. The model reproduces chimneys identified by seismic imaging without pre-defining their locations or geometry. The chimney locations were determined by the steepness of the interface between the reservoir and the caprock, the reservoir thickness, and the compaction length of the strata. We demonstrate that chimneys are highly-permeable leakage pathways. The width and propagation speed of a single chimney strongly depends on the viscosity and permeability of the rock. For the chimneys of the Snøhvit field, the predicted time of formation is about 13 to 40 years for an about 2 km high chimney.

1. Introduction

Fluid flow tends to be localized in space and time at all length scales within the Earth: from the deep mantle to the shallow subsurface [1,2,3,4]. Focused fluid flows have shaped the planet during its 4.5 billion years history [4,5]. Their variety on Earth is enormous: they often have similar geometry but different origins caused by various geological processes and mechanisms. The understanding of focused fluid flow is a significant challenge because of the non-linear nature of the problem. Today, seismic gas chimneys are among the most common manifestations of the focused fluid flows in sedimentary basins worldwide [6,7,8]. Understanding gas chimney formation and evolution are essential for oil and gas exploration and secure long-term storage of CO₂ [9,10,11].

Gas chimneys are vertical fluid escape pathways through sealing rock and overburden rock strata. They are observed in sedimentary basins worldwide [7]. Their morphology is non-unique in terms of scales, seismic signatures, and shapes [12,13]. Interpretations of deteriorated seismic signals provide their present-day location, while circular seafloor craters (pockmarks) indicate the locations in the present and past [8,14]. There are two fundamental mechanisms for the formation of a gas chimney. One is hydraulic fracturing driven by overpressure when fluid pressure exceeds the lithostatic. However, this mechanism would be limited to the upper 1 km of the overburdened rocks since the pressure does not allow for boiling pore fluids [15]. Another mechanism is related to ascending porosity wave propagation through poro-visco-elastoplastic media [16]. That is influenced by rock viscosity, porosity, permeability, and initial subsurface architecture.

Most current research on gas chimneys discusses geological processes leading to chimney formation, e.g., [1,4,17,18,19]. However, only a very few theoretical works with modeling on the reservoir and basin scales are presented. The numerical modeling of the gas leakage and chimney forming mechanisms is very challenging. Tasianas et al. [20] modeled the gas leakage along pre-existing conducting geological structures on the reservoir scale. They considered the flow along three types of structures preliminary set up in the model: faults, fracture network penetrating the reservoir’s caprock, and high permeable fluid conduits. However, they did not consider the processes leading to the formation of these structures. Ostanin et al. [21] considered the basin-scale mechanism of gas leakage by faults sealed during glacial loading and conductive during ice retreat. Duran et al. [22] used the basin modeling approach to estimate the volume of leaked gases from the reservoir to the seafloor. They assumed that flow pathways were formed by a capillary failure of the seal when the upward-directed buoyancy pressure, generated by the lower density of hydrocarbons (HC), plus any excess overpressure in the reservoir exceeds the capillary resistance pressure of the seal. Iyer et al. [15] considered a vent complex for gas leakage. Vent complexes are associated with large igneous provinces globally and are similar to gas chimneys. The vents are considered conduits for fluid and gas phases generated within the contact aureole of the associated magmatic intrusion. Iyer et al. [15] used a simplified ad-hoc modeling approach, assuming that rock permeability increases with increasing pressure and returns to the initial value as fluid pressure drops using the same function.

Few models [16,23,24] attempted to explain chimney formation using the porosity wave mechanism. The porosity wave can transfer fluid through geological layers with low permeability in localized areas with high porosity. The character of this process is spontaneous and self-propagating. Mckenzie [25] initially suggested this type of model to explain the migration of melt in partially molten rocks. Afterward [26], the same model was used for the understanding of sedimentary compaction. Yarushina and Podladchikov [16] created a simple model for porous fluid flow in a deformable matrix. The model can capture the range of rheological responses within Earth’s lithosphere. This model can predict the behavior of fluid-rock systems during sedimentary compaction at a wide range of temperatures and time scales. The model was limited by hydrostatic compaction and decompaction. Later [23], the effects of shear stresses on (de)compaction processes were considered. Based on this model, Rass et al. [24] presented a 3D simulation of the spontaneous development of high-permeability pathways through two geological layers from a fluid-enriched source region. These regions reproduce natural seismic pipes or chimneys. However, these models are based on synthetic datasets and do not consider geological processes on a basin-scale preceding the chimney formation.

This work addresses the formation of chimney-like focused fluid flow structures in the South-Western (SW) Barents Sea by considering coupled basin- and reservoir-scale modeling. We focused on the Snøhvit and Tornerose fields of the Hammerfest Basin in the SW Barents Sea. The objects’ choice is related to two factors: (1) a good understanding of basin and petroleum systems evolution and (2) the presence of gas chimneys traced from the main reservoir of Stø Formation in the Snøhvit field. Hence, the causal relationships of reasons, duration, and location of gas chimneys are determined [20]. In addition, this object is of great interest to CO₂ sequestration. The basin modeling approach provides the chimney location’s thermal, petrophysical, and structural characteristics with seismic data analysis. The followed reservoir-scale hydro-mechanical modeling of gas leakage associated with the porous rock deformations uses outputs from the basin model as an initial setup for simulations.

2. Geological Settings and Petroleum Systems

The geological setting of the epicontinental SW Barents Sea is well studied and described in works, e.g., [26,27,28,29,30], and briefly summarised here, noting the main geological events. SW Barents Sea has experienced a relatively deep basin burial history. The basin accounts for several tectonic episodes of lithosphere extension accompanied by subsidence and basin formation from Devonian to Paleocene-Eocene time [26]. The basin tectonic history begins with the Caledonian Orogeny covering the territory from Greenland to the Norwegian shelf [31]. Extensional episodes from Devonian to Early Cenozoic time with accompanied subsidence formed the segregated SW Barents Sea basins with block-faulted structures presented by highs and troughs [26,32]. Compressional events formed the fault complexes along the boundaries of major structural elements [33]. The Cenozoic period is introduced by uplifts and erosions in the basins and on structural highs, while at the shelf edge, thick progradation structures were deposited [34]. The deep erosion phase began in the Oligocene due to tectonic uplift after the Norwegian-Greenland Sea opening [35]. The rapid decrease of a global temperature in the northern hemisphere occurred at ~2.5 Ma and was accompanied by a regional glaciation [36]. The Pliocene to Pleistocene uplift and erosion associated with the isostatic response of cycle loading of ice-cap determined the late-stage evolution of the SW Barents Sea [34]. The Hammerfest Basin is one of the SW Barents Sea primary oil and gas regions. The basin is enclosed to the south by the Troms–Finnmark Fault Complex, to the west by the Ringvassøy–Loppa Fault Complex, while to the north and east by the Loppa High and the Bjarmeland Platform, respectively (Figure 1). The Hammerfest basin initially was structurally continuous with the Loppa High [35,37]. The Carboniferous extension caused the tilting of the Loppa High and Hammerfest basin to the Early Permian. The isolated sea deposition environment with anoxic episodes formed the Kobbe formation (Fm) as a possible source rock.

During the Late Triassic, the formed accommodation space during subsidence was infilling by sediments. The Late Triassic siliciclastic Snadd Fm is considered a possible source rock. Tectonic subsidence, faulting reactivations, eustatic sea-level changes, and sediments’ deposition controlled the evolution of the margin from the Late Triassic to the Middle Jurassic. Simultaneously, sea-level changes formed the foreshore environment to offshore and estuarine sub-environments strongly influenced by wave processes and tidal effects [39]. The deposition formed the Tubåen Fm and Stø Fm as the main reservoir of the Snøhvit and Tornerose fields. During Middle Jurassic to Late Jurassic, global sea-level rise formed the marine environment for shale sedimentation, mainly the primary source rock Hekkingen Fm [29,35,40]. Transgression sequence in the central part of the Hammerfest Basin during the Middle Paleocene initiated the marine sedimentation. Progradational sedimentation bodies sourced from the platform areas to NNE filled the basin during the Late Paleocene [22].

The Jurassic formations Stø and Tubåen are the primary reservoirs in the Hammerfest Basin according to the petroleum system analysis of Duran et al. [22,40]. The Snøhvit field reservoir is located on the west side of the transect, and the Tornerose discovery is on the east side (see yellow ovals in Figure 2b). The reservoir of Stø Fm in the Snøhvit area consists of natural gas with an oil leg, while Tornerose discovery accumulates pure gas [22,40].

Several gas chimneys traced from the primary reservoir of Stø Fm in the Snøhvit field were observed by different researchers [4,21,22,37] (Figure 2b). Mohammedyasin et al. [37] applied high-quality pre-stack time-migrated (PSTM) 3D seismic data interpretations (Figure 1b) to analyze the gas leakage pathways. Geometrically, the chimneys are determined as tubular-shaped (Chimney 1), cone-shaped (Chimney 2), and Christmas tree-structured (Chimney 3) (Figure 2b). Based on the acoustic masking, the lateral extension of the chimneys varies from ~1 to ~10 km with decreasing depth. Chimney 1 and Chimney 2 decrease laterally downward, while Chimney 3 has a Christmas tree geometry.

Ostanin et al. [21] and Duran et al. [22] consider the Pliocene–Pleistocene glaciation cycles a primary reason for hydrocarbon loss from Stø Fm reservoir. Nevertheless, Ostanin et al. [21] consider the basin fault system a pathway for gas and fluids leakage during ice unloading. In contrast, Duran et al. [22] consider mainly the seal capillary failure during the peak of ice loading. In both works, gas leakage events are associated with several deglaciation periods from 0.02 to 0.01 Ma.

3. Data and Methods

We consider both basin-scale geological processes that preceded the chimney formation and porous fluid flow coupled with deformation processes associated with relatively short and local chimney formation episodes. To capture these processes at different time and length scales, we combine two approaches: basin modeling and coupled hydro-mechanical modeling at the reservoir scale. First, basin modeling allows reproducing the burial history of the basin. It produces the properties of the basin, such as porosity, effective thermal properties, permeability, paleotemperatures, and organic matter maturity rank. Second, we perform coupled hydro-mechanical modeling at the target reservoir-sized domain based on basin modeling results which allows reproducing spontaneous development of gas chimneys. The basin modeling provides the geological setting and the initial conditions for the chimney formation, which is handled within the second simulation step and used as the initial setup for the reservoir-scale hydro-mechanical model.

Chimneys that form due to non-linear hydro-mechanical coupling occur over less than 1 Ma, and its spatial resolution varies from 1 m to 1 km. The non-linear hydro-mechanical coupled model generates fluid escape chimney-like structures cutting through the sedimentary units. The porosity and permeability evolve in time and space within the reservoir-sized domain. However, there is no feedback of these changes back to the basin modeling simulator. Therefore, the impact of the chimneys on the paleo-fluid pressure is not considered.

3.1. Basin Modelling

3.1.1. Burial Histories

Basin modeling aims at reproducing the burial history, subsidence rate, erosion period and rate, sedimentation rate, and thermal history of a sedimentary basin [41]. The modeling considers present-day stratigraphy, including the sedimentary units’ geometry, deposition time, and lithology. The first step in a basin simulation is a calibration of the burial history. It is done by estimating the net thicknesses of (porosity-free) rock in each formation. These net thicknesses remain constant throughout the burial history. The net thicknesses are used in a forward simulation of the burial history, which includes the computation of the porosity of the sedimentary units at each time step. There are two approaches to calibrating the net (porosity-free) thicknesses in the burial history. The so-called “back stripping” method works backward in time by stripping off layer by layer and decompacting the remaining layers. See details in [42]. The decompaction is done using porosity as a function of burial depth. The other approach of estimating the formation thicknesses is to run forward simulations as an iterative process. The net thicknesses are updated at the end of the burial history by comparing them with the present-day thicknesses. The second approach is more general than the first and allows for porosity functions that depend on depth and time, temperature, and pressure [43]. The backstripping approach is decoupled from the simulation of paleo-temperature and paleo-pressure. It is, therefore, decoupled from the forward simulation [44,45,46,47]. In this study, we use the Athy [48] porosity function of depth

$$\phi \left(z\right)={\phi }_0·e^{-z/B},$$

(1)

where $\phi$ is the porosity at the depth $z$, ${\phi }_0$ is surface porosity, B is the compaction length (see

Table 1. Complete list of symbols used in the work.

Parameter	Description	Value	Unit
ϕ	Porosity
ϕ0	Initial porosity
z	Burial depth		km
B	Compaction length scale		km
ρ	Density		kg/m3
ρm	Density of rock matrix		kg/m3
ρw	Density of water at 20 °C	1000	kg/m3
ρeff	Bulk density		kg/m3
Cρ	Specific heat capacity		J/kg/K
Cρeff	Bulk specific heat capacity		J/kg/K
Cρm	Specific heat capacity of rock matrix at 20 °C		J/kg/K
Cρw	Specific heat capacity of water	4182	J/kg/K
k	Thermal conductivity		W/m/K
keff	Bulk thermal conductivity		W/m/K
kr	Thermal conductivity of rock matrix		W/m/K
kw	Thermal conductivity of fluids		W/m/K
Q	Radiogenic heat production		W/m3
T	Temperature		K
t	Time		Ma or yr
${\rho }_f$	Pore-fluid density	1020	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
${\rho }_s$	Solid density	2040	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
${\mu }_f$	Fluid shear viscosity	$8\times {10}^{-4}$	$\mathrm{P}\mathrm{a}·\mathrm{s}$
${\eta }_{\varphi }$	Effective solid bulk viscosity		$\mathrm{P}\mathrm{a}·\mathrm{s}$
$k$	Permeability		${\mathrm{m}}^2$
$L_c$	Compaction length		m
$t_c$	Viscous compaction time		yr
$g$	Gravity constant	9.8	$\mathrm{m}/{\mathrm{s}}^2$
$k_{\varphi }$	Dynamic permeability		${\mathrm{m}}^2/\mathrm{P}\mathrm{a}/\mathrm{s}$
${\mu }_s$	Solid shear viscosity		$\mathrm{P}\mathrm{a}·\mathrm{s}$
$v^s$	Solid velocity		$\mathrm{m}/\mathrm{s}$
$v^f$	Fluid velocity		$\mathrm{m}/\mathrm{s}$
$q^D$	Darcy flux		$\mathrm{m}/\mathrm{s}$
$\overline{\rho }$	Total porosity averaged density		$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
$P_f$	Fluid pressure		$\mathrm{P}\mathrm{a}$
$\overline{P}$	Total pressure		$\mathrm{P}\mathrm{a}$
${\tau }_{ij}$	Stress deviator		$\mathrm{P}\mathrm{a}$
${\delta }_{ij}$	Kronecker delta

). Each lithology has its surface porosity and compaction length.

In the previous work [49], particularly for the same case study, we compared the present-day state of the basin models derived by both approaches. We used the BMT software [50] to reconstruct the basin history by the backstripping method, while we use the simulator TecMod 2019.1 to reconstruct the basin with iterative forward simulations. Works [44,47] give a detailed description of BMT and TecMod simulators, respectively. Both approaches produced similar results from the point of view of thermal regime and organic matter maturity. Thus, in the present work, we use only the results of the forward modeling approach.

3.1.2. Thermal Histories

Here we briefly review the temperature equation used in basin modeling. A more detailed presentation of heat flow in sedimentary basins can be found in the books of Hantschel and Kauerauf [41], Wangen [43], and Allen and Allen [51].

The temperature equation for transient heat conduction and radiogenic heat production is

$$\rho C_p\frac{\partial T}{\partial t}=\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)+Q,$$

(2)

where $\rho$ is the bulk density, C_p is the bulk specific heat capacity, k is the bulk thermal conductivity, T is the temperature, z is the vertical coordinate, t is the time, and Q is the radiogenic heat production. The temperature equation is for simplicity written in one dimension for the vertical direction. Notice that the temperature equation does not have a heat convection term. Heat convection may often be ignored in basin simulations, but not always [52].

The bulk thermal conductivity k_eff is obtained as the geometric mean

$$k_{eff}=k_r^{\left(1-\phi \right)}·k_w^{\phi },$$

(3)

of the rock matrix thermal conductivity $k_r$ and fluid heat conductivity $k_w$ [53,54].

The Sekiguchi [55] temperature correction has been applied to rock matrix thermal conductivities. The volumetric bulk heat capacity is obtained by as the arithmetic mean of the volumetric heat capacities of the rock matrix and the fluid [41]

$${\rho }_{eff}{C}_{p_{eff}}={\rho }_m{C}_{p_m}\left(1-\phi \right)+{\rho }_w{C}_{p_w}·\phi ,$$

(4)

where ${\rho }_m$ and ${\rho }_w$ are rock matrix and water densities, respectively, and $C_{p_m}$ and $C_{p_w}$ are rock matrix and water-specific heat capacities, respectively. Both of them include temperature corrections by Waples and Waples [56] and Somerton [57].

A two-dimensional basin has four boundaries. The vertical boundaries are considered isolated, which implies that they have zero heat flux. The surface temperature is assumed known as a function of time. The basin base can either have a given temperature thermal gradient or heat flow as a function of time. The thermal state of the base of the basin is the result of simulations of the entire lithosphere when the basin is fully integrated. If there are no data for the temperature or heat flow at the basin base, the lithosphere can be included in the basin simulation. The lithosphere can be included as a passive domain or as a dynamic domain that undergoes stretching. Lithospheric stretching produces thermal transients that may affect the basin over time spans of several tens of million years [43].

The basin’s thermal model must be calibrated to the data obtained from wells of present-day temperatures and thermal maturity parameters, e.g., vitrinite reflectance [41].

3.2. Hydro-Mechanical Modeling at the Reservoir Scale

Porosity waves were proposed as a mechanism generating focused fluid flow in sedimentary basins [16,23,24,58,59]. They result from the coupling of porous fluid flow and viscous matrix deformations and decompaction weakening [60]. Rock rheology influences the shape of the porosity waves, which can take the form of elongated drops or chimneys [58,60,61].

We use the coupled hydro-mechanical equations to model the spontaneous formation and propagation of high-porosity and permeability channels at the reservoir scale [16]. These equations describe the filtration of incompressible pore fluid coupled with the deformation of the viscous solid matrix.

Mass balances for the rock and pore fluid can be written as

$$\frac{\partial {\rho }_s\left(1-\phi \right)}{\partial t}+{\nabla }_i\left({\rho }_s\left(1-\phi \right)v_i^s\right)=0$$

(5)

$$\frac{\partial {\rho }_f\phi }{\partial t}+\nabla \left({\rho }_f\phi v_i^f\right)=0.$$

(6)

The momentum balances for the solid and fluid phases has a form

$$\left({\nabla }_j{\tau }_{ij}-\overline{P}{\delta }_{ij}\right)-g_i\overline{\rho }=0,$$

(7)

$$\phi \left(v_i^f-v_i^s\right)+\frac{k_{\varphi }}{{\mu }_f}{\nabla }_i\left(P_f+{\rho }_f{g}_i\right)=0,$$

(8)

where ${\rho }_s,{\rho }_f$ are solid and fluid densities, $\overline{\rho }=$ $\phi {\rho }_f+\left(1-\phi \right){\rho }_s$ is the total porosity-averaged density, $\phi$ is the porosity, $v_i^s,v_i^f$ are the components of solid and fluid velocities, $\overline{P}, P_f$ are the total and fluid pressures, ${\tau }_{ij}$ are the components of the stress deviator, $g_i$ is the component of the gravity acceleration vector, $k_{\varphi }=k_0{\left(\frac{\phi }{{\phi }_0}\right)}^3$ is the Carman-Kozeny relation (the non-linear porosity-dependent permeability), ${\mu }_f$—the pore fluid viscosity, $\phi \left(v_f-v_s\right)=q_i^D$ is the Darcy flux vector.

Within this work, we consider viscous rheology, which is described by the following relation

$${\tau }_{ij}=2{\mu }_s\left(\frac{1}{2}\left({\nabla }_i{V}_j+{\nabla }_j{V}_i\right)-\frac{1}{3}{\delta }_{ij}{\nabla }_i{V}_i\right),$$

(9)

where ${\mu }_s$ is the solid shear viscosity.

The viscous (de)compaction is described by the constitutive equation [16]

$${\nabla }_i{v}_i^s=-\frac{\overline{P}-P_f}{\left(1-\phi \right){\eta }_{\varphi }},$$

(10)

where ${\eta }_{\varphi }$ is the solid bulk viscosity, which is porosity dependent and reduced by decompaction weakening factor in the case of overpressure. Equation (10) is the closure relation for the governing system of Equations (5)–(10), which is used for the reservoir modeling of chimneys’ spontaneous formation and evolution.

Rock response is visco-elasto-plastic [62]. Elastic deformation gives an immediate response, while viscous deformation develops on a longer time scale. The relative importance of viscous and elastic terms is controlled by Deborah number introduced in [63]. In our case, the estimation of Deborah is based on the typical values of elastic moduli presented in [62]. Previous model studies show that viscous deformation is essential for forming focused fluid flow [60,64,65,66]. For Deborah numbers from the range of $1\times {10}^{-4}–1\times {10}^{-3}$ elastic properties don’t significantly impact the formation of focused fluid flow [63,67]. Therefore, in our simulations, we keep only viscous deformation which gives the first-order effect.

The process of focused fluid flow affected by viscous deformation is characterized by compaction length

$$L_c=\sqrt{\frac{k_{\varphi }{\eta }_{\varphi }}{{\mu }_f}}$$

(11)

and viscous compaction time

$$t_c=\frac{{\eta }_{\varphi }}{L_c\left({\rho }_s-{\rho }_f\right)g}$$

(12)

where $k_{\varphi }$ is the dynamic permeability, ${\mu }_f$ is the fluid viscosity, ${\eta }_{\varphi }$ is the effective bulk viscosity, $\Delta \rho$ is the difference in densities of the rock and fluid, $g$ is the gravity constant.

4. Modeling Results

4.1. Basin Modeling

4.1.1. Dataset from the Hammerfest Basin

A transect from the Hammerfest basin was simulated with TecMod2D 2019.1. The burial history begins with the deposition of the Ørret Fm in Permian, a little more than 250 Ma ago. Sixteen layers represent the stratigraphy from the Late Permian to the present-day seabed (Figure 2). In the model, each layer has its lithology (

Table A1. Material properties and lithological characteristics of modeled stratigraphic units (based on [22]).

			Physical Properties A
Layer	Color B	Lithology	ϕ0	B	ρ	k	Cρ	Q
				(km)	(kg/m3)	(W/m/K)	(J/kg/K)	(μW/m3)
Nordland Gp		Siltstone (organic lean)	0.55	1.96	2720	2.05	921	1
Torsk Fm		Shale (organic lean, typical)	0.70	1.20	2700	1.70	879	2
Kveite-Kviting fms.		Siltstone (organic lean)	0.55	1.96	2720	2.05	921	1
Kolmule Fm		Shale (organic lean, silty)	0.67	1.20	2700	1.77	879	2
Kolje Fm		Shale (typical)	0.70	1.20	2700	1.64	879	2
Knurr Fm		Shale (typical)	0.70	1.20	2700	1.64	879	2
Hekkingen Fm		Shale (organic rich, 8% TOC C)	0.70	1.20	2500	1.20	879	3
Fuglen Fm		Shale (organic lean, siliceous, typical)	0.70	1.20	2710	1.90	879	1
Stø Fm 01		Sandstone (typical)	0.41	3.23	2720	3.95	837	1
Stø Fm 02		Siltstone (organic lean)	0.55	1.96	2720	2.05	921	1
Nordmela Fm		Siltstone (organic lean)	0.55	1.96	2720	2.05	921	1
Tubåen Fm		Sandstone (clay poor)	0.42	3.33	2700	5.95	837	0
Fruholmen Fm		Siltstone (organic rich, 2–3% TOC)	0.55	1.96	2700	2.00	921	1
Snadd Fm		Siltstone (organic rich, 2–3% TOC)	0.55	1.96	2700	2.00	921	1
Kobbe Fm		Siltstone (organic rich, 2–3% TOC)	0.55	1.96	2700	2.00	921	1
Havert-Klappmys fms.		Shale (organic lean, silty)	0.67	2.33	2700	1.77	879	2
Ørret Fm		Siltstone (organic lean)	0.55	2.44	2720	2.05	921	1

^A. See symbols’ designation in complete list in Table 1. ^B Color codes correspond to the lithology units from Figure 3a. ^C Total Organic Carbon.

). There are no lateral changes of the lithologies, except for the main reservoir in the Stø Fm, which changes laterally from sandstone (Stø Fm 01) to siltstone (Stø Fm 02) at the lateral position of 28 km (Figure 3a). The lithology description of each layer and corresponding petrophysical properties are presented in Appendix A based on work [22]. The Oligocene-Miocene erosional event occurs between 29–15 Ma [22]. The eroded thicknesses range from 140 to 840 m, and the thicknesses are increasing from southwest to northeast. This erosional scenario corresponds to the average appraisal of eroded thickness. We follow this scenario since it produces the average porosity for the reservoir among scenarios with a deviation range not higher than ±2%. The glacial activity during Pliocene-Pleistocene is ignored in the model, even though the glaciation and/or deglaciation impact the pore fluid pressure. Nevertheless, we do not consider this event as a critical trigger for the chimney formation.

The lithosphere was included in the simulations. The initial thickness of the upper crust, lower crust, and mantle are 20 km, 20 km, and 100 km, respectively. The modeling of rifting processes assumes differential thinning of the crust and the mantle. The Hammerfest basin has experienced four rift phases, according to Reemst et al. [34] and Skogseid et al. [68]. The model includes only three of them (Triassic, Jurassic-Cretaceous, Palaeocene-Early Eocene) since the first rift phase (Devonian-Carboniferous) was earlier the first horizon in the model.

The boundary condition at the lithosphere base is a constant temperature of 1300 °C [69,70]. The basin surface is either subaerial or underwater. The surface is bounded by the temperature that varies between 10 to 25 °C through the Paleozoic-Pleistocene period [71]. Temperature does not go above 3 °C in the interglacial Pleistocene periods [72,73]. Therefore, the surface temperature is taken to be 6 °C at present [74].

The study uses TecMod’s automatically reconstructed paleo-water depth [47].

4.1.2. Calibration

The stratigraphy was automatically calibrated by the simulator [47,75]. Less than 10 forward simulations were needed to match the present-day layers’ thicknesses. The relative difference between the simulated and observed formation thicknesses is less than 1%.

The calculated thermal parameters were calibrated against temperature and vitrinite reflectance data from a series of wells. These wells are 7120/6-1, 7121/4-1, 7121/4-2, 7121/5-1, which are close to position 50 km on the profile, and to wells 7120/7-1, 7120/7-2, 7120/7-3, 7120/8-1; 7120/8-2 close to lateral position 29 km [23] (and links therein) (Figure 1b). The vitrinite reflectance was computed using the EASY-Ro model [76].

4.1.3. Results

According to Mohammedyasin et al. [37], three chimneys are observed only in the Snøhvit area tracing from the Stø Fm to the seafloor. For a chimney to leak out the gas in a reservoir, the gas must overcome the high capillary threshold of the seal. The Fuglen Fm is the primary seal of the reservoir in the Stø Fm [22,40]. The seal thickness of the Stø Fm in the Tornerose area is up to ten times higher than in the Snøhvit area (Figure 3a). The low thickness of the seal in the Snøhvit area makes it vulnerable for the Oligocene-Miocene erosional event and cycles of glaciation during the Pliocene-Pleistocene. The units above the seal are more coarse-grained, less compacted, and less cemented. These lithological units have properties that allow a chimney to form and to conduct the gas upward.

We observe the higher reservoir thickness with higher porosity of up to 18–22% in the area of Snøhvit field than in Tornerose field (Figure 3b). Thus we assume that this reservoir characteristic is the most optimal for chimney formation [49].

At Palaeocene, the maximum temperatures for the surface of Kobbe Fm source rock in the Snøhvit area reach ~185–190 °C, while the Tornerose area reaches not higher than 140–150 °C (Figure 4). A temperature higher than 170 °C indicates that the secondary cracking of oil to gas [77] was realized. Such thermal conditions contribute to the formation of thermogenic gas, increasing overpressure in the overlying rocks. Other source rocks, Snadd Fm and Hekkingen Fm, experienced lower peak temperatures than 170 °C for both fields (Figure 4).

Underneath the Snøhvit area, at the present day, the surface of the Kobbe Fm is located in the wet gas window, the surface of Snadd Fm—in the main oil window, while the Hekkingen Fm—predominantly in early oil window by Tissot and Welte [78] (see Figure 5). Underneath the Tornerose area, the temperatures and vitrinite reflectance values are significantly lower than in the Snøhvit area. Only the surface of Kobbe Fm is located in the main oil window. The oil in the Snøhvit reservoir area seems to have migrated from the Upper Jurassic Hekkingen Fm, which is in accordance with the petroleum system analysis of Duran et al. [22,40]. The high maturity of the Triassic Snadd Fm and Kobbe Fm source rocks is responsible for the high volumes of gas migrating to the Jurassic reservoirs.

4.2. Hydro-Mechanical Modeling at the Reservoir Scale

4.2.1. Initial Reservoir Model Based on Dataset from the Hammerfest Basin

To provide chimney formation simulation on the reservoir scale, we consider the target reservoir-sized domain (~$32\times 2.5$ km) of the basin model presented in Figure 6. This domain has a simplified structure of layers lumped by their lithology characteristics and porosity distribution given in the basin model.

Table 2. Properties of the layers in the base reservoir model.

No. of Layer	Lithology	Porosity	Dynamic Permeability	Effective Solid Bulk Viscosity
		(%)	$({\mathrm{m}}^2/\mathrm{P}\mathrm{a}/\mathrm{s})$	$(\mathrm{P}\mathrm{a}·\mathrm{s})$
1	Shale	15–18	$2\times {10}^{-12}–6\times {10}^{-12}$	$8\times {10}^{-15}–1\times {10}^{-16}$
2	Sandstone	18–22	$6\times {10}^{-12}–8\times {10}^{-12}$	$6\times {10}^{-15}–8\times {10}^{-15}$
3	Shale	9–13	$3\times {10}^{-13}–1\times {10}^{-12}$	$1\times {10}^{-14}–5\times {10}^{-15}$
4	Siltstone	18–60	$6\times {10}^{-12}–6\times {10}^{-10}$	$3\times {10}^{-15}–1\times {10}^{-16}$

summarises the ranges of porosity values corresponding to a specific layer of the target domain.

We use coupled hydro-mechanical model (5)–(10) presented in Section 3.2 to simulate the formation and propagation of high-porosity chimneys across the target region. The model does not take into account the fault system since, according to seismic investigations [37], the chimneys’ distribution zones are located far from faults. Thus, the faults are not the plumbing system for these three chimneys. Numerical implementation includes discretizing (5)–(10) using a finite-volume technique on a regular Cartesian grid in 2D in MATLAB. The computational domain has a resolution of $1048\times 256$ grid points. The mechanical part of (5)–(10) was solved using free-slip (no shear stress) boundary conditions on all sides of the computational domain. For the fluid flow part, we utilized no flux boundary conditions on all vertical sides of the computational domain, and fixed flux values at the bottom and top boundaries. Initial conditions for the numerical simulation are presented in Figure 7 and Figure 8, namely initial porosity, dynamic permeability, and bulk viscosity distribution. The initial porosity of the model (Figure 7a) corresponds to the values from the basin model. Figure 7b shows the initial dynamic permeability, $k_{\varphi }/{\mu }_f$, which varies in the range of ${10}^{-13}–{10}^{-11} {\mathrm{m}}^2/\mathrm{P}\mathrm{a}\cdot \mathrm{s}$ depending on the specific layer. The pore-fluid density is assumed to be $1020 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$. The solid is twice less buoyant than the pore-fluid. The fluid shear viscosity equals $8\times {10}^{-4} \mathrm{P}\mathrm{a}·\mathrm{s}$ [24].

Within this work, we performed a sensitivity analysis of the numerical model to the initial values of effective bulk viscosity, which is not very well constrained by the laboratory data [24]. For these purposes, we performed two simulations with different ranges of initial effective bulk viscosities. The considered values are taken from [24,62,80]. The initial effective bulk viscosity, ${\eta }_{\varphi }$, is presented in Figure 8 and has values in the range of ${10}^{14}–{10}^{16} \mathrm{P}\mathrm{a}\cdot \mathrm{s}$ in the first simulation and ${10}^{15}–{10}^{17} \mathrm{P}\mathrm{a}\cdot \mathrm{s}$ in the second one.

The process of chimneys formation and propagation affected by viscous deformation is characterized by compaction length and viscous compaction time given by Equations (11) and (12).

The fluid density is assumed to be twice less than the solid one. Since the target reservoir domain has a multi-layered structure, we have different values of compaction length (3–300 m and 10–1000 m in the first and the second simulation, respectively) depending on a layer considered in the system (Figure 9). Note, the characteristic time and length are used to non-dimensionalize the problem.

4.2.2. Chimney Modeling Results

Figure 10 and Figure 11 show the evolution of high-porosity chimneys produced in numerical simulations for two cases of lower and higher effective bulk viscosities of reservoir layers. We observe fluid flow episodes which form chimney-like patterns reflecting the multi-layered structure of the considered domain. Chimneys need about 13 years to reach the surface in case of lower effective bulk viscosity values (${10}^{14}–{10}^{16} \mathrm{P}\mathrm{a}\cdot \mathrm{s}$). For the other case, developed high porosity channels take the larger time to get to the surface, namely around 40 years. The existence of a low-permeable weak shale layer significantly delays the formation of focused fluid flow pathways, which is observed during both simulations. According to the simulation results, channels need about 8 and 25 years to rise through the shale layer in case of lower and higher effective bulk viscosities of reservoir layers, respectively. The resulting locations of high-porosity chimneys are represented in Figure 11a,b. Developed chimneys resulting from the simulation with higher effective bulk viscosities of reservoir layers are characterized by larger width, corresponding to the larger compaction length values. The numerical simulation results are compliant with the identified seismic chimney zones in the target reservoir-size domain shown in Figure 11c. In the middle of the considered domain, seismic data shows a wide chimney, while numerical results predict the formation of a large number of separate smaller channels. These smaller channels might not be resolved individually in seismic and might be seen as one large chimney instead. Results show that the formation of high-permeable channels is associated with topological gradient and the location of initial porosity anomalies. Besides, there is a tendency to accumulate porosity anomalies near lithological boundaries associated with the seal rock (Figure 10) with a following on porosity wave penetration through it.

5. Discussion

Based on basin modeling results, we can infer that chimney formation requires a combination of several factors. Firstly, the source rocks should be mature by being in the wet gas window. Hydrocarbon generation from overmature source rocks creates overpressure in the reservoir. Secondly, the reservoir rock should have a large capacity for long-term gas preservation under high pressure, which is the source for chimney formation. Thirdly, the seal should have low thicknesses (less than hundreds of meters), and additionally, it must lose its integrity by deformation from the tectonic movements and diagenesis.

At this step, we characterized the conditions of chimney zones formation from basin and petroleum system modeling points of view. Comparing the Tornerous and Snøhvit fields, it was found that the latter is in the most favorable conditions for the chimney formation. Next, for chimney modeling on the reservoir scale, we take from the basin model the area of the Snøhvit field. Thus, reservoir, seal rock and overburden rock geometry, and distribution of porosity are extracted.

The results of our basin modeling are different from the results presented in Duran et al. [22] and Ostanin et al. [21]. The previous works used a basal heat flow trend as the lower boundary condition to a basin bottom during the thermal history reconstruction (see Section 3.1.2). Thus, the temperature Equation (2) is solved only for a basin domain. The heat flow accounts only for a response from the lithosphere deformation but ignores the blanketing effect [81,82], when the deposition of “cold” sediments may lead to depressed vertical thermal gradients if the sedimentation rate is “large”. The work [46] demonstrated that using the lower thermal boundary condition that ignores the blanketing effect and thermal effect from erosion of the sediments leads to the estimation of HC generation up to 86%. Thus, the estimates of transformation ratio (TR), a mass of generated HC, and timing of expulsion, migration, and preservation for the analyzed source rocks Kobbe Fm, Snadd Fm, and Hekkingen Fm can contain similar errors despite satisfying well data calibration in works of Duran et al. [22] and Ostanin et al. [21]. Therefore, the gas leakage volumes can be overestimated too. The evidence of temperature overestimation over time without accounting for the blanketing effect is presented in Figure 12. In the Snøhvit area, the surface of Kobbe Fm by Ostanin et al. [21] has higher temperature values up to 45 °C, although the values are matched at present-day. The overestimation of temperature results in a higher maturity rank. For instance, the present-day value of Snadd Fm in the Snøhvit area is estimated as ~1.05 Ro% [21], while the estimates of the present work are close to 0.8 Ro%. Therefore, to obtain the consistent thermal history of the basin, it is required to account for lithosphere deformation during basin evolution and to set the lower boundary condition at the lithosphere bottom to consider the blanketing effect. Otherwise, the calibration of the thermal model against well data is only localized fitting.

Another point of the criticism is connected with the suggested mechanisms of chimney formation. The earliest petroleum system models [22,83] proposed the mechanism of gas leakage as capillary failure in the absence of faults. The capillary failure happened due to the pressure increasing during several glaciation cycles. However, Ostanin et al. [21] proven that the overpressure does not reach the fracture gradients for the Fuglen Fm and Hekkingen Fm seal rocks. Furthermore, in the present work, we showed that the source rocks formations’ temperature history was overestimated, leading to overestimating overpressure, inducing the capillary failure.

Thus, the basin model of Ostanin et al. [21] considers the primary mechanism of gas chimneys formation due to the fault system. Faults were assumed to be partially sealing or conductive during glacial unloading due to significant deviatoric stresses and isostatic compensation and closed during glaciation due to subsidence. According to the modeling, the peak pressure coincides with peak glaciations. In a previous study by Ostanin et al. [84], the fault system of the same dataset was studied in great detail. Also, several amplitude anomalies associated with the gas clouds have been identified based on the high RMS acoustic impedance contrast. Authors claim that all the anomalies are crosscut and follow the strike of faults or are bounded by the 1st order faults [84]. Based on the study of results in [79], we suppose that the conclusion that gas chimney shape follows the faults’ strike or bounded by the 1st order faults is far-fetched. Some of the anomalies are identified out of any 1st and 2nd order faults. The azimuth of the gas cloud does not coincide with the azimuth of the 1st order faults, also as the fault inclination has no similar geometry to the gas cloud. In addition, the mechanism of four and more cycles of fault drainage activation and deactivation by uniform ice load is justified but only assumed. Also, the fault properties depend on several factors, including the degree of disaggregation, dissolution, cementation and mineral precipitation, cataclasis, and stress state [85,86,87].

In the work of Mohammedyasin et al. [37], the plumbing system is interpreted from seismic not only by soft reflections but also by acoustic masking. The soft reflections are interpreted into two groups of anomalies: structurally unconformable to the background reflectors and structurally conformable (aligned in the same direction as the background reflectors). The last is stratigraphically controlled and has no spatial relationship with the deep-seated faults. The acoustic masking analysis helped to delineate three chimneys, as shown in Figure 2b. The author of the work emphasized that the geometry of chimneys signifies fluid leakage at different times for Christmas-Tree (Chimney 3) and regular and consistent leakage for cone-shaped (Chimney 1 & 2). Because of the large lateral extension of each chimney from 1 to 10 km, the authors suggest that more hydrocarbons were leaked from the reservoir through the chimneys. However, in conclusion, the authors give a primary role of pathways for gas leakage to deep-seated faults because most of the trapped gas is on the hanging wall along the section of the major faults [37]. This justification has a weak point since the trapped gas indicates that escaping the gas is difficult. In contrast, the evidence of a chimney with a uniformly distributed not trapped gas confirms the evidence of regular and consistent leakage of fluid in the study area.

In the present work, we consider that the chimney takes the primary role of the gas escape from the reservoir with the porosity wave mechanism. In contrast, the leakage by fault drainage, capillary failure, and molecular diffusion takes the secondary role. We agree with the previous studies that the glacial loading creates a triggering overpressure to leak the fluids. Therefore, it launches the propagation of porosity waves. There are examples where the chimneys are observed, but the glacial loading lacks [88,89,90,91,92,93,94]. Other overpressure-forming mechanisms such as the organic matter maturation, avalanche shelf sedimentation [93], or salt tectonic [94], could play a role in these regions.

To simulate localized highly-permeable channels in the SW Barents Sea, we use the reservoir-sized coupled hydro-mechanical model with the initial conditions coming from the basin modeling. The model suggests porosity waves in porous media with decompaction weakening as a primary mechanism for the spontaneous formation of chimneys. Simulation results reproduce natural observations of chimney-like focused fluid flow structures seen in seismic data. Besides, reservoir simulations show the presence of smaller channels that are possibly below the seismic resolution. The location of modeled chimneys is controlled by stratigraphic layers’ topology, thickness, and compaction length. Our results show that the chimney’s width and propagation speed strongly depend on the compaction length, $L_c$, given by Equation (11), and viscous compaction time, $t_c$, given by Equation (12), which are functions of material parameters of layers such as effective bulk and shear viscosity and permeability. Reducing permeability or increasing porous fluid shear viscosity by one order of magnitude will reduce the compaction length by factor 3 and increase the time of chimney formation by the same factor. At the same, reducing solid bulk viscosity by one order of magnitude will reduce both the compaction length and the time of chimney formation by the same factor 3. The essential feature of our work is a coupling of the reservoir- and basin-scale modeling and prediction of spontaneous formation of chimney structures without artificially pre-defined flow pathways.

6. Conclusions

In this paper, we combine basin and reservoir scale models to understand the nature of seismic chimney as focused fluid flow.

A reservoir-scale hydro-mechanical model, being based on geological settings produced by basin modeling, allows for the predicting of chimney-like focused fluid structures in the reservoir cap rock. The model involves fluid flow coupled with viscous deformation of permeable rock, with a decompaction weakening effect.

The results of basin modeling suggest that the lack of chimneys above the Tornerose field is due to insufficient pressure from immature source rocks, in addition to poor temperature conditions for secondary cracking.

The basin model of the Snøhvit field provides a geological environment favorable for chimney formation. The huge gas clouds above the Snøhvit formation are considered gas-filled chimneys, which are produced by porosity wave propagation. It is concluded that chimney formation is the primary reason for the assumed gas leakage from the reservoir.

The steepness of the reservoir-caprock interface, the thickness of layers, location of initial porosity anomalies are found to be a reason for the high-permeable chimneys and their localization. The chimney formation takes about 13 to 40 years to reach the seafloor, depending on the considered values of bulk viscosity of the layers, while the primary resistance occurs in the seal rocks.