Evidence for massive emission of methane from a deep-water gas field during the Pliocene.

Geologic hydrocarbon seepage is considered to be the dominant natural source of atmospheric methane in terrestrial and shallow-water areas; in deep-water areas, in contrast, hydrocarbon seepage is expected to have no atmospheric impact because the gas is typically consumed throughout the water column. Here, we present evidence for a sudden expulsion of a reservoir-size quantity of methane from a deep-water seep during the Pliocene, resulting from natural reservoir overpressure. Combining three-dimensional seismic data, borehole data and fluid-flow modeling, we estimate that 18-27 of the 23-31 Tg of methane released at the seafloor could have reached the atmosphere over 39-241 days. This emission is ∼10% and ∼28% of present-day, annual natural and petroleum-industry methane emissions, respectively. While no such ultraseepage events have been documented in modern times and their frequency is unknown, seismic data suggest they were not rare in the past and may potentially occur at present in critically pressurized reservoirs. This neglected phenomenon can influence decadal changes in atmospheric methane.

The present‐day atmospheric methane (CH4) budget accounts for a variety of natural sources, including geologic processes, such as natural gas seepage from petroleum‐bearing sedimentary basins and fluid manifestations in geothermal areas (1, 2). Hydrocarbon seeps located onshore, shallow offshore, and coastal areas are the main geologic sources, while CH4 released by deep (≳300–400-m) ocean seeps typically does not reach the atmosphere as the majority of the released gas dissolves and is oxidized in the water column (3, 4). However, massive deep‐water seepage events could contribute a substantial amount of CH4 to the atmosphere, even from depths >1000 m (5–7). Acoustic and seismic imaging provides evidence for such events in the form of giant pockmarks on the seabed and seepage anomalies in subsurface sedimentary formations (8).

In addition to contributing to atmospheric CH4, massive deep‐water seepage events would also constitute an important source of temporal variability in geologic emission rates. Whereas nongeologic sources, such as wetlands and biomass burning, are known to exhibit important interannual or long‐term trends (9, 10), geologic sources have generally been assumed to be constant in time for purposes of CH4 budgeting (10, 11) and in investigations of preindustrial atmospheric CH4 isotope ratios (12). However, the rate of natural‐gas emissions from geologic sources is highly variable on multiyear to geologic timescales (2, 13). Massive episodic seepage events from terrestrial and shallow‐water sources have likely played a role in past climate changes (13–15).

Here, we present evidence for a sudden and massive deep‐water release of CH4 in the Faroe‐Shetland Basin (FSB; northeastern margin of the Atlantic Ocean) during the early Pliocene. We consider high‐resolution three-dimensional (3D) seismic data and borehole log and core data, where the latter were used to calibrate the geophysical observations and to constrain the reservoir properties. We combine these data with fluid‐flow modeling to estimate that 23–31 Tg of CH4 was expelled from a single subsurface reservoir over the course of 39–241 d, and that a substantial fraction of this CH4 would have reached the atmosphere.

Results

Geophysical Evidence for Gas Seepage.

Geophysical and well data confirm the present‐day occurrence of a hydrocarbon accumulation in a gas field within the FSB. This hydrocarbon accumulation, belonging to the Tobermory gas field (TGF), is hosted in a sand‐rich fan deposit known as the Strachan Fan (16) (Fig. 1 ). The reservoirs of this field have thermogenic gas with CH4 concentrations >99 vol% and a methane‐to‐ethane ratio of ∼125, and that originates from highly mature source rocks of the Upper Jurassic Kimmeridge Clay Formation (17; ). Gas generation occurred at about 175 °C (18). At this temperature, the isotopic signature of methane, for any type of kerogen, is typically enriched in 13C relative to the atmosphere, with δ13C values exceeding −43‰ (19, 20).

Fig. 1.

Location and direct and geophysical evidence of the Tobermory gas field (TGF) in the Faroe‐Shetland Basin (FSB). (A) The 3D seismic data cover the northern portion of the Strachan Fan. The TGF is one of many oil and gas fields located offshore of the United Kingdom. (B) Well 214/4–01 intersects the top Strachan Fan (TSF) at a depth of 2,660 m MD (measured depth) and encountered a gas–water contact (GWC) at 2,712 m MD. The Strachan Fan represents the reservoir interval (R.) and is overlain by a polygonally faulted clay‐rich seal and a biosiliceous mud‐rich overburden. The overburden, limited by the Upper‐Eocene Unconformity (UEU) and the Intra‐Neogene Unconformity (INU), hosts the seepage zone (S.Z.). Gas is encountered from the seal to the Paleocene lava flow (PLF), suggesting present‐day charge of hydrocarbons to the reservoir. The gas is methane dominated, although ethane is present. Interval velocity ranges from 1,790 to 2,512 m/s. These velocities were used to convert interpreted surfaces to depth on the time‐migrated seismic data. (C) The GWC, identified by the response of the rms amplitude map and related to the TGF, conforms to the structure of the TSF. The anomaly of the GWC is partially attenuated by the overlying S.Z. (Fig. 2).

Fig. 2.

Seepage zone and paleopockmarks above the TGF. (A) The seepage zone is located above the prominent GWC (see Fig. 1). The amplitude response of the individual amplitude anomalies (blue–red–blue [b–r–b]) is opposite with respect to the seabed (red–blue–red [r–b–r]) (Insets) and, thus, indicates an anomalous local decrease of acoustic impedance with depth. The top of the seepage zone is limited by the INU, implying that the expulsion of methane terminated at the end of the emplacement of this erosional event. (B) The INU exhibits eight prominent erosive depressions distributed above the GWC and interpreted as paleopockmarks. (C) Closeup of paleopockmark 4. The paleopockmarks are composed of circular subdepressions, which are interpreted as the vents or vent sites where gas was expelled. (D) Amplitude map along Horizon X (see Fig. 2) showing aligned amplitude anomalies indicating the presence of conduits in these specific regions.

Two geophysical observations provide evidence for a single vigorous release of CH4 from this reservoir during the Pliocene: 1) A region of acoustic amplification above the gas accumulation indicates the presence of compressible fluids (gas) and is interpreted as a seepage zone, and 2) eight irregular depressions in the paleoseabed above the seepage zone are interpreted as pockmarks formed by vigorous fluid venting.

The seepage zone is evidence that gas was vented from the reservoir. It is limited at the top by the Intra Neogene Unconformity (INU; Fig. 2), a regional hiatus formed by the action of the North Atlantic Deep Water current flowing southeastwards through the FSB and dated as being Early Pliocene in age (21). The seepage zone has a volume of 1.5 to 4.6 × 107 m3 (Fig. 2 and see ) and is consistent in appearance and structure with other seismic observations of gas‐expulsion phenomena above gas fields in other basins worldwide (22). The extent of the seepage zone coincides closely with the seismically identified extent of the prominent gas–water contact associated with the underlying TGF (Figs. 1 and 2). This implies that the individual amplitude anomalies composing the seepage zone indicate shallow occurrences of CH4 that migrated from the underlying reservoir.

The pockmarks are evidence that gas venting was vigorous. These eight irregular depressions are mapped at the INU directly above the seepage zone (Fig. 2). These depressions range from 700 m to 2.4 km across and are characterized by an erosive base and an onlap fill (Fig. 2). Whereas the INU formed by deep‐water erosion, these depressions have closed perimeters and are, therefore, unlikely to have formed by the action of bottom currents. Instead, the depressions are interpreted as paleopockmark craters based on their similarity to pockmarks described elsewhere (8). The erosive character of these depressions (Fig. 2) suggests that fluid vented through the seabed at a sufficiently large velocity to mobilize and excavate the seafloor sediments to a depth of 40–50 m (Fig. 2). This characteristic, typically observed at the top of blowout pipes (23) and hydrothermal vents (24), indicates a single occurrence of rapid and sustained fluid expulsion from the subsurface rather than slow seepage (3). The crestal position of these large depressions relative to the gas field further supports their interpretation as paleopockmarks, suggesting that they were the main venting sites for the underlying seepage zone (Fig. 2) and that high‐pressure fluids (CH4 and water) migrated vertically from the underlying TGF into the ocean.

Lastly, the occurrence of anomalies only up to the paleoseabed (the INU) and the occurrence of pockmarks only at the paleoseabed are evidence that gas was expelled in a single event rather than via gradual or intermittent seepage.

Estimate of the Amount of CH4 Released into the Ocean.

We assess the impact of this massive release by estimating: 1) the mass of CH4 stored in the TGF prior to the release, 2) the mass of CH4 remaining in the reservoir, the seepage zone, and the conduits after the release, 3) the mass of CH4 released into the ocean, and 4) the mass of CH4 emitted to the atmosphere. To incorporate the uncertainty and potential spatial variability associated with these estimates, we assign a statistical distribution to each physical input quantity (see and ), and we propagate these distributions throughout our analysis via a Monte Carlo framework. In the main text, we report input quantities as the mean ± the SD of their assigned distribution, and we report the resulting estimates as ranges from the first quartile value to the third quartile value.

We estimate the mass of CH4 stored in the TGF prior to the release based on its structure, a four‐way‐dip anticlinal trap (Fig. 1). This trap developed progressively through the Neogene due to in‐plane compression related to ridge push from the northeastern Atlantic spreading axis (17). We reconstruct the trap morphology at the time of the release using the fossilized opal amorphous to cristobalite/tridimite transition (opal A–CT) present in the study area (25). The opal A–CT is a diagenetic boundary associated with the dissolution and reprecipitation of biosiliceous sediments and is observed on seismic data as a strong‐amplitude seabed‐simulating reflection with the same polarity as the seabed (26) (Fig. 2). In the FSB, the opal A–CT was active during the Pliocene (25), so the trap morphology prior to the release can be reconstructed by determining the deformation of the opal A‐CT with respect to its active position. This reconstruction shows that the anticlinal trap was already in place prior to the release but structurally shallower than its present‐day configuration by 9 ± 5.8 m and with a gross capacity of 7.7–8.6 × 108 m3 (see ); note that this constitutes a relatively small gas field.

Direct calibration with borehole 214/4–01 confirms that the present-day GWC is located at the spill point (i.e., the trap is currently full). The paleolocation of the gas–water contact is difficult to determine, but the regional source rock has been producing hydrocarbons since the late Cretaceous (); we, therefore, assume that the trap was also full at the time of the expulsion. We also assume that the porosity and gas saturation in the reservoir prior to the expulsion were uniform and similar to the present‐day values of 0.34 and 0.93, respectively, as determined during the drilling of well 214/4–01 (Fig. 1 and ), which has an off‐axis intersection with the TGF. We, therefore, estimate that, prior to the release, the trap contained 2.1–2.4 × 108 m3 of CH4 or 40–46 Tg of CH4 at paleoreservoir pressure and temperature (P–T) conditions (see ).

We next estimate the mass of CH4 remaining in the reservoir, the seepage zone, and the conduits after the release. We estimate the former quantity by assuming that all mobile gas is expelled from the reservoir such that only residual (trapped) gas remains. The mass of CH4 remaining in the reservoir after the release is then a fixed fraction of the prerelease mass. We estimate a residual gas saturation of 0.16 ± 0.023, a reasonable but conservative range given that values of <0.12 have been measured in the laboratory (27, 28). The mass of CH4 stored in the seepage zone is less well constrained due to uncertainty in the associated gas saturation (29). We, therefore, use a wider range of gas saturations, 0.20 ± 0.05, acknowledging that seismic amplitude is insensitive to variations in gas saturation (30). For the conduits, we calculate the gross volume based on their geometry (see ) and otherwise applying the same methodology as for the seepage zone. We, thus, estimate that 6–7 Tg of CH4 remain in the reservoir and 5–11 Tg remain in the seepage zone and the conduits, suggesting that the balance, 23–31 Tg, was released into the ocean (see ).

Estimate of the Duration of Venting and the Amount of CH4 Released into the Atmosphere.

We now estimate the duration of this emission and the mass of CH4 that would have reached the atmosphere. The duration is important because CH4 has an atmospheric lifetime of ∼10 y (2, 31) and a 20-y global warming potential 86 times that of CO2. To estimate these quantities, we first consider the flux of CH4 from the reservoir to the paleoseabed (the INU).

The existence of the gas reservoir suggests that the overburden, a faulted sequence rich in both clay and biosiliceous mud (Figs. 1 and 2), acted as an effective seal until and following the release. The release was most likely driven by regional overpressure due to tectonic activity, a common geological process in which regions of the subsurface are compressed by the motion of tectonic plates. The widespread presence of synchronous sand injectites between the Strachan Fan and the paleoseabed suggests that this overpressure occurred relatively suddenly and over a large region (). These sand injectites were emplaced in a single short‐lived event immediately after the formation of the INU (32). Injectites form when the pore pressure in an interval rich in unconsolidated sand exceeds the total vertical lithostatic stress at that depth, causing hydraulic failure of the immediate overburden and, then, liquefaction and injection of sand (33). Injectites can be dated via the forced folding of the overburden that occurs when the sand slurry intrudes in the form of shallow sills (34). In the FSB, the forced folding due to shallow injectites and the paleopockmarks above the seepage zone occur along the same seismic horizon, calibrated as the INU (). This supports the argument that the formation of the injectites and the seepage from the gas reservoir were approximately synchronous () and, importantly, that both events were sourced from the sand‐rich Strachan Fan (Fig. 1 and ). Note that the gas in the reservoir would have been a small part of this regional fluid expulsion.

We, therefore, assume that the regional pressure that drove both the formation of sand injectites and the release of CH4 was equal to the total vertical lithostatic stress of 25–27 MPa, which corresponds to a fluid overpressure of 7–9 MPa (see ). Note that the total vertical lithostatic stress is a lower bound since higher pressures could have occurred prior to seal breach (34 and see ). Tectonic activity is the most plausible mechanism for the sudden generation of a regional overpressure of this magnitude; other common causes of overpressure, such as fluctuations in sea level or changes in sedimentation rate, are thought to be too gradual and too modest (35). For example, an overpressure of this magnitude would have required a sudden drop in sea level of >700 m, which is unlikely. Tectonic compression is further evidenced by the large number of extensive fold structures related to tectonic activity in the Neogene (17, 36, 37) and as observed in many other basins (38, 39). Note that the gas would not have invaded the permeable seal before this pressure was reached because of the large capillary entry pressure (>100 MPa based on composition, 40).

Hydraulic failure of the overburden would most likely have occurred by exploiting preexisting weaknesses, such as faults and polygonal faults (41). Once opened, these conduits would have provided a high‐conductivity pathway for the release of CH4 into the ocean (42, Fig. 3). The presence of aligned sand extrudites and amplitude anomalies in the shallow section of the seepage zone also supports the idea of venting through linear features, such as faults (Fig. 2 and ). As a result, we conceptualize venting as the flow of a gas–sand suspension through fault-like conduits that open rapidly at the beginning of the expulsion and remain open for the duration of gas expulsion before closing as the overpressure eventually relaxes (Fig. 3). We then use a flow model based on the Darcy–Forchheimer equation to estimate that the expulsion of CH4 into the ocean was completed in a period of 39–241 d (see ).

Fig. 3.

Conceptual model of the pre‐, syn‐, and postseepage events. (A) Preseepage condition with gas emplaced in the Strachan Fan deposit. No bypass was present at the time, and CH4 accumulated up to the spill point. The opal A‐CT transition was active and seabed simulating (25). (B) Overpressure, probably triggered by tectonic compression (see the main text), caused the opening of conduits by exploiting preexisting weaknesses in the seal (Inset). This process allowed rapid CH4 expulsion to the seabed (INU) and the formation of both pockmarks and seepage zone. (C) Recharge of the TGF. No further seepage phenomena occurred after the INU. UEU = Upper‐Eocene Unconformity.

The release took place at an estimated water depth of 1,090 ± 86 m (see ). Gas consumption by the ocean is thought to be effective at preventing CH4 from reaching the atmosphere (3, 5, 6), but the rate of consumption depends strongly on the size of the gas bubbles and on the local thermodynamic conditions. Typical underwater seeps produce relatively small gas bubbles (diameter <∼50 mm) (12). In the present context, we estimate the bubble size from the modeled rate of expulsion using an empirical relation, resulting in diameters of 48–98 mm (see ). We model the effect of expansion, dissolution, and gas exchange for bubbles of this size (see ) to find that 18–27 Tg (79–93%) of the CH4 released into the ocean would have reached the atmosphere (Fig. 4). This represents 44–62% of the CH4 stored in the TGF (Fig. 4). With a volumetric methane/ethane ratio of ∼125 (17), this emission would have included about 0.15–0.22 Tg y−1 of ethane, a photochemical pollutant and ozone precursor. The tightness of these results (recall that each range is from the first to the third quartile from a Monte Carlo analysis) suggests that our estimates are robust to substantial uncertainty in the input parameters within the constraints of our modeling assumptions.

Fig. 4.

Results and impact of the emission to the atmosphere. (A) Boxplot showing median (red dash), interquartile range (IQR, blue box), and lower and upper 1.5 × IQR whisker (black line) of the mass of CH4 from the simulated geologic model. The amount of methane emitted into the atmosphere (A) from the TGF is 44–62% of the total initially present in the reservoir (R). The balance is either dissolved into the ocean (O) or remains as residual in the reservoir (Resid.), the conduits (C), or the seepage zone (SZ). (B) The emission of CH4 from TGF would account for ∼10% of the present‐day natural emissions (top down estimate in ref. 43). The emission from TGF represents about four times the net interannual variation of CH4 from all sources (43).

Impact of the Gas Seepage on Atmospheric CH4.

The abrupt release of 18–27 Tg CH4 is an extraordinarily large emission compared to ordinary rates of individual seepage sites, which are typically on the order of Mg or a few Tg CH4 per year (1, 2). We propose the term “ultraseepage” to define such reservoir‐size gas expulsions. The emission corresponds to about 10% of the current annual emission from all natural sources (Fig. 4) and to about 28% of current global annual CH4 emissions from oil and gas operations (43). Its climate impact is equivalent to >1 mo of total current US anthropogenic CO2 emissions (43) for a 100-y atmospheric impact time horizon, or >4 mo for a 20-y time horizon. Using a box model of the global atmosphere, this emission would have resulted in a nearly instantaneous increase of ∼9 ppb CH4 in the atmosphere globally (see and ). Since it is not a sustained CH4 source, most of this increase would have been removed over the following 10–20 y (44). This period is short compared to the time resolution of preindustrial ice core CH4 measurements. An event of this magnitude today, however, would be readily detectable with global monitoring stations, which have observed interannual variations over the past three decades of about 5–13 ppb (45). The ethane release (0.15–0.22 Tg y−1) is about 10% of present global emissions from all geosources (46, 47).

Discussion

As explained above, this gas expulsion was driven by higher-than‐lithostatic pore‐fluid overpressure in the main reservoir of the TGF, most likely due to tectonic activity. Overpressures of this magnitude are not uncommon, as evidenced by the globally diverse occurrence of reservoir‐scale injectites (32). In addition, a number of different mechanisms commonly lead to overpressures of lower magnitude (35); for example, critical fluid-overpressure conditions can be triggered by seismicity, as observed for mud volcanoes (48). These events could still be sufficient to cause similar atmospheric impacts in shallower settings than the one studied here (35). Seismic data from other basins show several other examples of potentially rapid and massive gas release due to catastrophic seal failure, and similar or even much larger CH4 release events could also occur in modern times (49, 50). The emission studied here is, therefore, not an isolated and exceptional event; rather it represents a previously neglected type of seepage that can have a substantial global environmental–atmospheric impact (11).

Our results suggest that further work is needed in a number of areas. Although we have allowed for substantial uncertainty in our estimates based on fluid-flow modeling, it is clear that the pressure-driven opening of conduits, the high-speed flow of gas and sand through these conduits and into the ocean, and the evolution of bubble plumes from deep sources are topics that require substantial further study to better constrain estimates like the one presented here. Our results also highlight the importance of identifying existing reservoir‐size accumulations of light hydrocarbons that are under critical overpressure and seal-failure conditions. Studying and potentially monitoring these accumulations as is currently performed for on-shore sources, such as mud volcanoes (48), will provide important insight into the frequency of these ultraseepage events, which is central to better assessing their potential occurrence and impact.

Materials and Methods

Geophysical and Well Data.

The 3D seismic data were acquired and processed by Petroleum Geo-Service in 2000. The data were processed with a standard processing sequence for marine seismic data (51) and finalized using a time‐Kirchhoff migration, zero phase, and American polarity (an increase in acoustic impedance with depth is associated with a positive reflection amplitude, RC + peak). The in‐line and x‐line spacing is 25 m and the vertical resolution, based on 1/4 of the dominant wavelength (52), is 7 m for the interval of interest (0–3,000 ms two-way travel time [TWT]). The seismic data were interpreted using Schlumberger’s Petrel software.

Well 214/4–1 was completed in June 1999 by Mobil North Sea Ltd. The well reached a total depth of 4,110 m (true vertical depth). The well data used in this study comprise density and γ‐ray logs completed after the perforation of the borehole, cuttings derived from the drilling of the geologic formations, and cored rock samples from the reservoir interval (core data, ).

We determine the present-day depths of the seabed, the INU, the UEU, and the reservoir from our geophysical data calibrated against well 214/4–1 (Fig. 1). We determine the paleodepths of these features by estimating and subtracting the subsidence of the basin , which we vary as part of our Monte Carlo analysis (). We assume throughout that the water pressure is hydrostatic and that the temperature is ∼0 °C at the seabed/paleoseabed and increases linearly with depth with a constant geothermal gradient ().

Statistical Methods.

A Monte Carlo framework was developed in order to calculate the mass of CH4 in the reservoir (before and after the emission), the conduit zone, the seepage zone, the ocean (see below, “Bubble-size modeling”), and the atmosphere (see below, “Atmospheric impact”). To obtain a statistically stable result, a total of 106 samples were used for each parameter ().

Model.

The Monte Carlo framework is applied to the geologic model derived from the interpretation of the 3D seismic and well data. The conceptual model is as follows:

A certain amount of gas was stored in the TGF prior to the emission ().

A triggering event resulted in a sudden and substantial overpressure (above lithostatic), producing fault opening, sand remobilization, and fluid expulsion (Fig. 3).

The fluid expulsion depleted part of the reservoir, produced a seepage zone above the reservoir, and emitted a certain amount of gas into the ocean at a certain rate.

The gas in the ocean was transported to the atmosphere via bubbles.

The TGF was recharged after the emission.

The presence of gas in the reservoir prior to the emission is proven by the formation of the seepage zone itself (Fig. 2 and ). The overpressure event triggering the emission is confirmed by the formation of sand injectites (). The emission to the ocean is calculated from the mass balance of the residual gas remaining in the reservoir, in the conduit zone, and in the seepage zone. The emission to the atmosphere is based on bubble modeling constrained by the gas flux estimated from the overpressure conditions (above lithostatic).

Volume and Mass of CH4.

The gross volume of the reservoir was obtained by calculating the volume enclosed between the top Strachan Fan surface and the isodepth surface intersecting the spill point of the top Strachan Fan surface. This operation was completed using Schlumberger’s Petrel software by applying a standard workflow, which includes interpretation of the seismic horizons, gridding (convergent interpolation), and calculation of the volume enclosed between two surfaces.

The net volume of gas in the reservoir was calculated as ϕ, where is the net to gross associated with the lithological proportion of sand over the total volume of the rock interval, is the porosity, and is the gas saturation. These parameters are chosen from distribution functions defined using real data (). The mass of CH4 is obtained by multiplying by the density of CH4 at paleo P-T conditions.

The gross volume of the seepage zone is calculated by counting the number of seismic voxels with amplitude above the amplitude threshold that separates anomalous from background amplitudes (). To capture the uncertainty behind the determination of anomalous versus background amplitude values, an amplitude cutoff is chosen for each simulation from a normal distribution as part of our Monte Carlo framework (). The voxel volume is given by the bin size 25 × 25 m, multiplied by the time‐to‐depth conversion of a time sample (4 ms) as function of a chosen interval velocity—a parameter chosen from a distribution function (). The net volume and the mass of CH4 in the seepage zone is calculated using the same approach used for the reservoir. The parameter used to constrain the properties of the biosiliceous mudstone hosting the seepage zone, namely, net to gross and porosity, are based on core data from equivalent geologic formations retrieved at the Ocean Drilling Program (ODP) site 643 (53). In order to capture the uncertainty about the gas saturation in the seepage zone, a distribution function is defined based on previously modeled gas saturations in the leakage zone (29; see ). The mass of CH4 is obtained by the product of the net volume of CH4 in the seepage zone with the density of CH4 at paleo P-T conditions.

We assume that the conduits close after the overpressure has dissipated and, therefore, that a negligible amount of CH4 is trapped and stored in the conduit region.

Conduit Properties.

We observe a cluster of eight pockmarks at the paleoseabed (the INU), each of which contains a series of vent sites (46 vent sites in total). For example, Pockmark 4 and the nine associated vent sites are shown in Fig. 2. For each pockmark, we assume that all of the associated vent sites are connected to the reservoir by the same conduit comprising several fault-like segments, each with its own aperture and across-slope length (33 segments in total across all conduits; ). We assume that each of the individual vents is circular with diameter . We assume that the across-slope length of each conduit segment is equal to the linear distance between neighboring vent sites as measured from 3D seismic data (). We take the total rectangular inlet area at the top of the reservoir across all conduits to be . We take the total outlet area at the paleoseabed across all conduits to be , where is the total number of vent sites across all conduits. The actual transition between each rectangular inlet and the associated circular outlets is likely to be complex; for simplicity, we assume that the total cross-sectional area of the conduit zone tapers linearly from to . We adopt a distribution of the dip angle of the conduits (relative to horizontal) that is consistent with the dip of the widespread regional polygonal-fault system (). The vertical extent of the conduits is the vertical distance between the top of the reservoir and the INU. We generate a distribution of the along-slope height of the conduits by dividing the vertical extent by .

We assume that, in response to the overpressure, the conduits would have opened relatively suddenly to their full aperture at the beginning of the expulsion. To produce a distribution of values for our Monte Carlos analysis, we first estimate a set of aperture values from the set of across-slope segment lengths using an empirical relationship between and for normal faults derived from observations at outcrops, (54). We, then, use the statistical properties of this set (min, max, mean, and SD) to create a truncated normal distribution for , the mean and SD of which are about 0.205 and 0.065 m, respectively. For each iteration of our Monte Carlo analysis, we use a single value of for all 33 of the conduit segments and all 46 vent sites. We assume that the aperture remains fixed for the duration of gas expulsion before closing as the overpressure eventually relaxes.

Overpressure.

The formation of the sand injectites and the forced folding requires a pore‐fluid pressure that is greater than the total lithostatic vertical stress at the depth where the injectites are present (32). We calculate the lithostatic stress at a paleodepth equal to the vertical center of the gas column in the reservoir by measuring the thickness of the overburden and by using a distribution for the vertical density based on logs from well 214/4–01, leading to a value of 25–27 MPa. The pore-fluid overpressure prior to the emission, calculated by subtracting the hydrostatic pressure from the vertical lithostatic stress, was 7–9 MPa. As a lower bound, we assume that the reservoir pressure during CH4 release was equal to this lithostatic stress prior to seal breach (34). We neglect the variation in this overpressure due to the changing gas column, which introduces an error of a few percent into our estimate of the gas flux but allows for a much simpler calculation.

Flux along the Conduits.

We assume that a mixture of gas and sand is expelled through the conduits. We model this flow using the Darcy–Forchheimer equation,

where is the pressure gradient along the conduit, is the volume flow rate through the conduit, is the cross-sectional area of the conduit, is the pressure difference between the reservoir (, assumed to be lithostatic as discussed above) and the seabed (, assumed to be hydrostatic), is the density of the gas–sand mixture, is the viscosity of the gas–sand mixture, is the Darcy permeability of the conduit, and is the Forchheimer (inertial) permeability of the conduit.

The flow properties and depend on the conduit aperture , whereas the fluid properties and depend on the volume fraction of sand in the gas–sand suspension . The latter is unknown, so we allow it to vary from 0 to the maximum value of (55). Lower values may be more likely since our evidence suggests that the conduits have closed completely, but this full range provides a more conservative estimate since sand increases both the effective viscosity and the effective density of the suspension (see below).

For and , which are relevant to slow (viscous/Darcy) flow, we use

where is the viscosity of CH4. We calculate a single value of for each iteration of our Monte Carlo method based on paleo P-T conditions at the vertical center of the conduit using the NIST Chemistry WebBook (56). This expression for is the well-known one for the permeability of a channel with aperture . This expression for is due to Eilers (57) for a fluid–solid suspension. For and , which are relevant to fast (inertial/Forchheimer) flow, we use

where is the density of CH4 and is the density of sand mineral (quartz). For each iteration of our Monte Carlo method, we calculate values of at several key depths (e.g., within the reservoir, middle of the conduit, at the seabed) based on paleo P-T conditions using the NIST Chemistry WebBook (56). This expression for is from ref. 58. This expression for is a basic volumetric average. We use this same mixture density in the gravity term.

To derive an expression for the venting rate, we integrate Eq. from the top of the reservoir to the paleoseabed. We take the volume flow rate and the fluid properties to be uniform within the conduits, and we assume that the total cross-sectional area varies linearly in from at the top of the reservoir to at the paleoseabed. The total flow rate through the conduits is then given implicitly by

where . Equation is readily solved via the quadratic equation and reduces to the expected result for Darcy flow as . The volume rate of gas expulsion is .

Bubble-Size Modeling.

Given that a volume flow rate of gas is vented from venting sites, we next estimate the size of the gas bubbles that would have been generated at the point of venting. The high-speed venting of gas–sand mixtures from natural conduits into deep water is poorly understood, with most experimental and field observations limited to low-flux gas seeps. We estimate the bubble size using an empirical correlation between the flow rate per venting site and the bubble diameter d (59),

where is a dimensionless empirical constant. Previous experimental work suggests that (59). We use a smaller value of to account for the breakup of bubbles due to turbulence at large venting rates and the presence of sand. For each iteration of our Monte Carlo method, we assume that all of the expelled bubbles are the same size; variation in bubble size is, therefore, captured across the full set our Monte Carlo results but not within each individual iteration.

Once the bubbles are formed, they will rise through the water column. As they rise, they evolve due to the exchange of gases with the surrounding water and the decrease in pressure with depth. Existing models for bubble evolution are complex and heavily parameterized (e.g., 3, 4). We develop a simple model by combining essential ingredients from a few existing models.

The terminal rise velocity of an individual bubble is typically parameterized as a function of its size. We do so using the work of Davies and Taylor (equation 2.4 of ref. 60). We relate the size of the bubble to the local pressure using the ideal gas law with a compressibility factor for pure methane at 0 °C (equations 12 and 13 of ref. 61).

In our setting, the gas-hydrate stability zone (GHSZ) extends from below the paleoseabed up to about 285 m below mean sea level. As a result, bubbles in the water column spend most of their rising time within the GHSZ. We, therefore, assume that a hydrate shell grows instantaneously on the outside of the bubbles at the sediment–water interface and that this shell disappears instantaneously when the bubbles reach the upper edge of the GHSZ.

We assume that bubbles contain CH4, N2, and O2, with an initial composition that is pure CH4. We model gas exchange with the ocean using the standard model (e.g., equation 1 of ref. 62) with associated mass-transfer coefficients (equations 4–6 of ref. 62). We reduce the mass-transfer coefficient by 70% in the GHSZ to account for the hydrate shell (section 3.2 of ref. 63). We calculate the maximum (saturated) concentration of N2 and O2 using Henry’s law with appropriate Henry’s constants (3). The solubility of CH4 is a strong function of temperature and pressure. Above the GHSZ, we calculate the maximum concentration of CH4 by fitting a cubic polynomial to experimental data (Table 4 of ref. 64). Below the GHSZ, we assume that the maximum concentration of CH4 is roughly independent of depth and given by its value at the top of the GHSZ (e.g., section 5.2 and Fig. 4 of ref. 65). We integrate the resulting system of coupled ordinary differential equations in MATLAB with the built-in solver ODE45.

Note that we have neglected all bubble–bubble interactions, both for rising and for mass transfer. We have also assumed a one-way coupling between the bubbles and the ocean by assuming constant background concentrations of dissolved gas and a constant size-dependent drag. Most of these effects, such as the net increase in the background concentration of CH4 in the water column and the enhanced rise velocity of bubble plumes relative to individual bubbles, would have increased the amount of CH4 that reached the atmosphere.

Atmospheric Impact.

The atmospheric impact of the emission () was computed using a one‐box model of the global atmosphere in steady state such that

where is the global average CH4 mixing ratio at year , is global annual CH4 emissions, and is the global average atmospheric CH4 lifetime (9 y). The model initialization of baseline global total CH4 emissions of 200 Tg ⋅ y−1 and 65- ppb CH4 is consistent with previous analyses (11, 12).