Kinetic reaction modelling of thermogenic emissions

The aims are to determine time-series of carbon emissions flux for comparison with PETM climate records (previous work reports final total mass of emissions17,25), to cover the full observed range of sill thicknesses and emplacement depths (previous work covers thinner, shallower sills25), and to calculate methane generation directly (previous work uses vitrinite reflectance as a proxy25,50). We combined several established methods to achieve these aims. All calculations were done in 1 spatial dimension: depth in the case of a sill, or horizontal distance in the case of a dyke. 1D calculations were used because the aim was to produce a simple parameterisation flexible enough to cover the observed range of sill and host-rock characteristics yet fast enough to be run very many times for stochastic modelling. The dyke calculations were used to test the numerical thermal calculations against analytical solutions for dykes, and to test the reaction calculations by comparison with published thermal aureole measurements for dykes.

Thermal maturation calculations were begun many millions of years before igneous intrusion to simulate burial heating and consequent maturation, which has a significant effect on final aureole width after intrusion. Methane already generated by burial maturation at the time of igneous intrusion is not included in post-intrusion emissions. During this stage, temperature is given by

$$\begin{array}{*{20}{c}} {T = T_0 + G\left( {z + Bt} \right)} \end{array},$$ (4)

where T 0 is the seabed temperature, G is the geothermal gradient, B is the burial rate and z is depth. A geothermal gradient of G = 30 °C km–1, burial rate of B = 10–5 m yr–1 and seabed temperature of T 0 = 10 °C were used for all calculations (all notation given in Supplementary Table 2). This G value may underestimate the true value in the Vøring and Møre Basins51 but may overestimate the true value in the Rockall Basin, where extreme crustal extension has occurred.

Instantaneous sill intrusion was assumed33. Thermal evolution was then tracked by accounting for conductive cooling of the igneous sheet and host rock according to

$$\begin{array}{*{20}{c}} {\frac{{\partial T}}{{\partial t}} = \kappa \frac{{\partial ^2T}}{{\partial z^2}}} \end{array},$$ (5)

where κ is the thermal diffusivity. Before solidification, this equation was solved using the method given in §4.19 of ref. 27, which accounts for latent heat of crystallisation. After solidification, Eq. 4 was solved using an FTCS explicit finite difference scheme. A uniform distance grid was used, and the grid spacing was kept small enough to generate smooth emissions flux time-series: typically 1 m up to 250 years post intrusion, 5 m between 250 and 10,000 years after intrusion, and 20 m at greater times. Seabed temperature was held fixed. Constant heat flux was assumed at the boundary/ies far from the sill/dyke, and the distance grid was made long enough to minimise edge effects. Burial heating after intrusion was not considered because only times of <1 Myr after intrusion are of interest here, during which burial heating is negligible in comparison with the thermal effect of the igneous sheet. At times well after complete solidification, an analytical solution is available to check these calculations (ref. 27, §4.21).

Four thermal maturation reaction groups were tracked31,32. The first group, RG1, describes thermal maturation of labile (oil-prone) kerogen to produce mostly oil and a smaller amount of methane. The second group, RG2, describes thermal maturation of refractory (gas-prone) kerogen to produce methane. The third group, RG3, describes thermal maturation of vitrinite, leading to increasing reflectance. The fourth group, RG4, describes cracking of oil (produced by RG1) to produce methane. Each reaction group covers a natural range of kerogen composition, with a corresponding range of activation energies. Therefore, each group is described by three parameters: a frequency factor A E , a mean activation energy E a and a standard deviation of activation energies σ (Supplementary Table 2).

Our reaction modelling implements the kinetic models of refs. 31 and 32. We approximated each normal distribution of activation energies with an odd number N E of discrete reactions, each with a constant activation energy and frequency factor. This strategy is more straightforward to implement and quicker to execute than numerical integration because progress of each reaction at each time step can be calculated using the analytical expressions. The concentration of a product P (i.e. the progress of the reaction) for one of the discrete reactions in RG1 to RG3 is given by

$$\begin{array}{*{20}{c}} {\frac{{{\mathrm{{d}}}\left( {1 \, -\, P} \right)}}{{{\mathrm{{d}}}t}} = - \bar k\left( {1 - P} \right)} \end{array},$$ (6)

where \(\bar k\) is the average reaction rate. We assume \(\bar k\) to be constant within each short time step of the finite difference calculation; the harmonic mean of the reaction rates for the temperatures at the start and end of the time step was used. The reaction rate at any temperature is given by the Arrhenius equation

$$\begin{array}{*{20}{c}} {k = A_{\mathrm{E}}{\mathrm{{exp}}}\left( { - \frac{{\overline {E_{\mathrm{a}}} }}{{R_{\mathrm{g}}T}}} \right)} \end{array},$$ (7)

where A E is the frequency factor, \(\overline {E_{\mathrm{a}}}\) is the mean activation energy of the discrete reaction and R g is the gas constant. The reaction progress at the end of the time step is then given by

$$\begin{array}{*{20}{c}} {P = 1 - \left( {1 - P_0} \right){\mathrm{{exp}}}\left( { - \bar k\,t} \right)} \end{array},$$ (8)

where P 0 is the progress at the start of the time step. We assumed that the total amount of kerogen was initially evenly distributed over the discrete reaction in each group, so that overall progress of the reactions group at each time step is given by

$$\begin{array}{*{20}{c}} {P = \mathop {\sum }\limits_{i = 1}^{N_{\mathrm{E}}} \frac{{P_i}}{{N_{\mathrm{E}}}}} \end{array}.$$ (9)

The discrete oil cracking reactions in RG4 are governed by

$$\begin{array}{*{20}{c}} {\frac{{{\mathrm{{d}}}\left( {1 \, - \, P} \right)}}{{{\mathrm{{d}}}t}} = - \bar k\left( {1 - P} \right) - C} \end{array},$$ (10)

where C(t) is the total rate of oil generation from the RG1 reactions. We assume C to be constant over each short time step so that

$$\begin{array}{*{20}{c}} {P = 1 - \frac{C}{{\bar k}} - \left( {1 - P_0 - \frac{C}{{\bar k}}} \right){\mathrm{{exp}}}\left( { - \bar k\,t} \right)} \end{array}.$$ (11)

We also assume that the oil generated at each time step is distributed evenly over the discrete oil cracking reactions, so that the terms involving C become \(C/\bar kN_{\mathrm{E}}\). The individual oil cracking reactions were then summed as before to find the reaction progress for the group.

A discrete reaction within any reaction group covers a range of activation energies between E 0 and E 1 given by

$$\begin{array}{*{20}{c}} {\mathop {\smallint }

olimits_{E_0}^{E_1} {\mathrm{{\Pi}}}\left( E \right){\mathrm{{d}}}E = \frac{1}{{N_{\mathrm{E}}}}} \end{array},$$ (12)

where Π(E) is the normal probability density function. Substituting this function and integrating gives

$$\begin{array}{*{20}{c}} {E_1 = E_{\mathrm{a}} - \sigma \sqrt 2 \,{\mathrm{erf}}^{ - 1}\left[ {{\mathrm{erf}}\left( {\frac{{E_{\mathrm{a}} - E_0}}{{\sigma \sqrt 2 }}} \right) - \frac{2}{{N_{\mathrm{E}}}}} \right]} \end{array}.$$ (13)

Thus, the values E 0 and E 1 that bound each of the N chunks of the overall activation energy distribution can be found by stepping across the distribution. The mean activation energy for each discrete reaction is then given by

$$\begin{array}{*{20}{c}} {\overline {E_{\mathrm{a}}} = N_{\mathrm{E}}\mathop {\smallint }

olimits_{E_0}^{E_1} E\,{\mathrm{{\Pi}}}\left( E \right){\mathrm{{d}}}E} \end{array},$$ (14)

which evaluates to

$$\begin{array}{*{20}{c}} {\overline {E_{\mathrm{a}}} = N_{\mathrm{E}}\left[ {\frac{{ - \sigma }}{{\sqrt {2\pi } }}{\mathrm{exp}}\left( {\frac{{ - \left( {E_{\mathrm{a}} - E} \right)^2}}{{\sigma \sqrt 2 }}} \right) - \frac{{E_{\mathrm{a}}}}{2}{\mathrm{erf}}\left( {\frac{{E_{\mathrm{a}} - E}}{{\sigma \sqrt 2 }}} \right)} \right]\begin{array}{*{20}{c}} {E_1} \\ {E_0} \end{array}} \end{array}.$$ (15)

The method developed here is a simple way of achieving our aims of directly estimating time-series of methane emissions flux over the observed range of individual sill dimensions. Additional effects such as latent heat of mineral organic maturation25, mineral dehydration reactions25 and detailed thermal characteristics of the host rock50 are included in other schemes for estimating thermogenic methane generation. Our method could be augmented to include these effects, though we anticipate that two other effects are more significant and will require more urgent attention. First, there is uncertainty in the proportion of generated methane released from the solid Earth (Discussion). Secondly, our NAIP sill province simulations show that 25% of sill aureoles interact with each other (Fig. 7g, h), implying that a wider range initial conditions for thermal structure and organic reaction progress should be covered in future calculations.

Parameterisation of thermogenic emissions results

We sought an analytically solvable expression with minimum free parameters for use in stochastic modelling. We found that the cumulative mass of carbon m therm (Δt) generated by thermal maturation adjacent to an igneous sheet can be conveniently parameterised as

$$\begin{array}{*{20}{c}} {m_{{\mathrm{therm}}}\left( {{\mathrm{\Delta }}t} \right) = M_{{\mathrm{therm}}}\left( {1 - {\mathrm{exp}}\left[ { - \left( {\frac{{{\mathrm{\Delta }}t}}{{\tau _{{\mathrm{decay}}}}}} \right)^p} \right]} \right)} \end{array},$$ (16)

where Δt is the time since intrusion, M therm is the total mass of carbon released by thermal maturation when the sheet has fully cooled and τ decay is the characteristic decay time for gas release. The flux of thermogenic methane is then found by differentiating this expression to give

$$\begin{array}{*{20}{c}} {q_{{\mathrm{therm}}}\left( {{\mathrm{\Delta }}t} \right) = p\left( {\frac{{M_{{\mathrm{therm}}}}}{{\tau _{{\mathrm{decay}}}}}} \right)\left( {\frac{{{\mathrm{\Delta }}t}}{{\tau _{{\mathrm{decay}}}}}} \right)^{p - 1}{\mathrm{{exp}}}\left[ { - \frac{{{\mathrm{\Delta }}t}}{{\tau _{{\mathrm{decay}}}}}} \right]} \end{array}.$$ (17)

This parameterisation is empirical but it reflects the real processes operating within the thermal aureole of the sill quite closely. The cumulative mass of carbon generated must increase until it tends to a constant value when the sill has cooled sufficiently to inhibit further maturation reactions. The analytical solution for conductive heating of the host rock shows that the thermal front moves outwards as a power law with p = 0.5 (ref. 27). Initially, we expect the thermogenic reaction front to move out from the sill at a similar rate to the thermal front but with a p exponent slightly >0.5 because reactions can take place at lower temperatures when more time is available. For the same reason, the emissions decay time will be of the same order but slightly larger than the conductive cooling time given by

$$\begin{array}{*{20}{c}} {\tau _{{\mathrm{cond}}} = \frac{{S^2}}{\kappa }} \end{array},$$ (18)

where S is the sill thickness and κ is the thermal diffusivity of the host rock. At later times, equivalent to maturation at distances further from the sill, the maximum temperature experienced by the host rock becomes exponentially smaller27. Thus, the expressions for m therm and q therm behave as power laws when Δt < τ decay and taper asymptotically to M therm or zero when Δt > τ decay .

The singularity at q therm (Δt = 0) arises because of the assumption of instantaneous intrusion. This assumption is reasonable because the intrusion period is on the order of days whereas the gas generation occurs over decades to thousands of years. To avoid the singularity when plotting at time t « 1 yr and when quoting initial fluxes, we average q therm over a preceding time period \(\bar t\) using

$$\begin{array}{*{20}{c}} {\overline {q_{{\mathrm{therm}}}} \left( {{\mathrm{\Delta }}t,\bar t} \right) = \frac{{M_{{\mathrm{therm}}}}}{{\bar t}}\left( {{\mathrm{exp}}\left[ { - \left( {\frac{{\left[ {{\mathrm{\Delta }}t \, - \, \bar t} \right]}}{{\tau _{{\mathrm{decay}}}}}} \right)^p} \right] - {\mathrm{exp}}\left[ { - \left( {\frac{{{\mathrm{\Delta }}t}}{{\tau _{{\mathrm{decay}}}}}} \right)^p} \right]} \right)} \end{array}.$$ (19)

An advantage of the empirical form developed above is that the parameter set (M therm , τ decay , p) required to specify m therm and q therm is straightforwardly related to the set of parameters (A, K, S, Z, w CH4 , β) that describes the dimensions of the igneous sheet and the characteristics of the host rock. These are: A, the surface area of the sheet; Z, the depth of intrusion of the sheet below the contemporary seabed; S, the maximum thickness of the sheet; β, an exponent that describes how the sheet thins towards its edges; w CH4 , the weight% organic matter expelled as carbon after thermal maturation; and K, the proportion labile:(labile + refractory) kerogen, where labile kerogen is the term for organic matter that initially produces oil which then cracks to gas, and refractory kerogen is the term for organic matter that produces gas directly31.

The total carbon mass parameter is determined using

$$\begin{array}{*{20}{c}} {M_{{\mathrm{therm}}} = A\,S\,Y\,d_{{\mathrm{host}}}w_{{\mathrm{CH}}4}} \end{array},$$ (20)

where Y is the final width of the thermal maturation reaction aureole scaled in units of maximum sheet thickness S, and d host is the density of the host rock. Note that our definition of Y includes the thermogenic reaction aureoles either side of the sheet (but not the igneous sheet itself), whereas some other studies use the scaled thickness of the aureole on one side only25. Our definition implicitly accounts for significant differences in thickness of the aureoles above and below a sill (Fig. 2), which arise because the speed of the contact maturation reactions is enhanced by the increase in pre-intrusion burial maturation with depth, and retarded by the fixed sediment surface temperature. These effects give a fourfold variation in Y across the depth range of interest (Supplementary Fig. 1). The functions Y crack (S, Z) and Y ref (S, Z), the scaled reaction aureole widths for methane from oil cracking and refractory kerogen maturation, respectively, were determined from multiple runs of the coupled thermal-kinetic reaction model that spanned the range of observed S and Z values (e.g. Fig. 2). Y was then found by averaging using

$$\begin{array}{*{20}{c}} {Y = K\,Y_{{\mathrm{crack}}} + \left( {1 - K} \right)Y_{{\mathrm{ref}}}} \end{array}.$$ (21)

Sills that vent fluid from their aureoles are intruded into the shallowest part of the sedimentary column, which is most affected by compaction. Sill emplacement depth therefore affects carbon emissions flux via host-rock density d host . The standard relationship between density and porosity is

$$\begin{array}{*{20}{c}} {d_{{\mathrm{host}}} = \phi \,d_{{\mathrm{water}}} + \left( {1 - \phi } \right)d_{{\mathrm{grain}}}} \end{array},$$ (22)

where d water = 1030 kg m–3 and d grain = 2650 kg m–3 are the densities of pore water and solid rock grains, respectively. The standard exponential porosity–depth relationship is

$$\begin{array}{*{20}{c}} {\phi = \phi _0\,{\mathrm{{exp}}}\left[ { - \frac{Z}{\lambda }} \right]} \end{array},$$ (23)

where we use a depositional porosity of ϕ 0 = 0.61 and a compaction length scale of λ = 2 km, both appropriate for mudrocks. These expressions show that the density at the observed mean emplacement depth (1.2 km) is 2110 km m–3, which was used for the calculations in Figs. 2, 3 and Supplementary Fig. 1.

The carbon emissions parameters τ decay and p are simple to relate to the sill dimension parameters. The emissions decay time is expressed as

$$\begin{array}{*{20}{c}} {\tau ^ \ast = \frac{{\tau _{{\mathrm{decay}}}}}{{\tau _{{\mathrm{cond}}}}}} \end{array}.$$ (24)

The functions τ*(S, Z) and p(S, Z) were determined from multiple runs of the coupled thermal-kinetic reaction model that spanned the range of observed S and Z values (Supplementary Fig. 1). As anticipated, τ* values are slightly greater than τ cond and p values are slightly greater than 0.5.

To account for variation in thickness across each sill within the rapid 1D calculation framework, we assumed radial symmetry in thickness, subdivided each sill into five annuli of equal surface area, calculated M therm and τ decay using the mean thickness within each annulus, and summed fluxes for each annulus to find the flux from the whole sill.

Parameterisation of carbon emissions from magma degassing

The total mass of carbon dissolved in the magma is

$$\begin{array}{*{20}{c}} {M_{{\mathrm{magma}}} = A\,S\,d_{{\mathrm{magma}}}w_{{\mathrm{{CO2}}}}} \end{array},$$ (25)

where d magma = 2750 kg m–3 is the density of the magma and w CO2 = 0.5% × 27%, the former being the typical wt% CO 2 and the latter the amount of carbon in the CO 2 molecule. Carbon is strongly incompatible in all primary minerals formed during solidification of basaltic sills. Almost all carbon dioxide dissolved in the magma therefore exsolves shortly before solidification, mixes with the super-critical host-rock pore fluid, and escapes through hydrothermal vents along with the thermogenic methane. The proportion of solidified sill to magma is a roughly linear function of time, so the mass of magmatic carbon exsolved is

$$\begin{array}{*{20}{c}} {m_{{\mathrm{magma}}}\left( {\Delta t \, < \, \tau _{{\mathrm{solid}}}} \right) = M_{{\mathrm{magma}}}\;t/\tau _{{\mathrm{solid}}}} \end{array}$$ (26)

and the flux of magmatic carbon emissions is

$$\begin{array}{*{20}{c}} {q_{{\mathrm{magma}}}\left( {\Delta t \, < \, \tau _{{\mathrm{solid}}}} \right) = M_{{\mathrm{magma}}}/\tau _{{\mathrm{solid}}}} \end{array},$$ (27)

where τ solid is the time to full solidification, which we parameterised as

$$\begin{array}{*{20}{c}} {\tau _{{\mathrm{solid}}} = 0.0147 - 0.000141\;S + 0.00472\;S^2} \end{array},$$ (28)

which gives τ solid in years for S in metres. After the sill has fully solidified, we have

$$\begin{array}{*{20}{c}} {m_{{\mathrm{magma}}}\left( {\Delta t \ge \tau _{{\mathrm{solid}}}} \right) = M_{{\mathrm{magma}}}} \end{array},$$ (29)

$$\begin{array}{*{20}{c}} {q_{{\mathrm{magma}}}\left( {\Delta t \ge \tau _{{\mathrm{solid}}}} \right) = 0} \end{array}.$$ (30)

Sill and host-rock observation database

Our primary set of surface area observations comes from several 3D seismic reflection surveys close to the centre of the NAIP. These were checked with reference to 2D data from across the NAIP. Sill lengths measured on 2D surveys (e.g. Supplementary Fig. 2) were converted to areas assuming a circular planform. The sill area histogram for each 2D survey was corrected for the fact that the sill chord intersected by a randomly oriented seismic line rarely represents the true sill diameter. For both 3D and 2D data, we excluded surface areas < 10 km2 because hydrothermal vents are rarely observed associated with such small sills, so their methane likely migrates over a timeframe longer than that of interest here. Sill surface areas measured from satellite images of the Karoo sill province are shown for comparison in Fig. 4a but were not included in the surface area distribution used for stochastic modelling.

Depths of intrusion were obtained from the subset of sills for which the coeval seabed can be identified from associated vents and/or forced folds. Measured emplacement depths were backstripped to correct for post-emplacement compaction, using the same porosity model parameters as for the host-rock density calculation.

Maximum sill thickness was measured for the subset of sills that show seismic reflections from both top and base, on both 3D and 2D datasets. The maximum sill thickness used in the stochastic modelling was obtained by first applying the linear regression to the diameter equivalent to the stochastically chosen surface area (assuming a circular sill), and then perturbing this result by a thickness selected from a uniform distribution between ±100 m.

We measured sill thickness as a function of radius both directly, where top and bottom reflections were visible, and indirectly from forced folds in the contemporary seabed. Thickness profiles near sill tips can only be estimated from forced folds because the resolution of oil-industry seismic data prohibits reflections from the bottom of the sill when its thickness is smaller than about 70 m. Fold-derived thickness profiles were corrected for compaction owing to post-intrusion sediment deposition. Other factors such as flexural host-rock deformation, anomalous cementation, fold-crest erosion and interaction with adjacent forced folds could also possibly corrupt sill thickness profiles derived from folds. We minimised influence of these complicating factors by only selecting examples where the thickness profiles from decompacted folds and direct sill imaging agree within measurement uncertainty. Maximum thickness estimated from the forced fold was used to scale the thickness profiles for both the forced fold and the directly imaged sill. The distance between the fold-derived thickness maximum and the fold-derived edge of the sill is then used to scale distance. The two profiles from a sill centre to the tips on either side were scaled separately in terms of distance.

Stochastic carbon emissions modelling procedure

The stochastic modelling algorithm has nine steps. The first step is to determine τ repeat , the time elapsed since intrusion of the previous sill: either specify τ repeat a priori or allow τ repeat to vary as a function of mantle plume flux, sill intrusion density data and sill province area data. Secondly, define a new sill–vent system by drawing eight random numbers for use in steps 3–5; physical properties thus specified remain fixed for subsequent time steps. Thirdly, determine new sill dimensions: define sill surface area A, depth of intrusion Z, maximum thickness S and thickness profile exponent β using the cumulative probability distributions in Fig. 4. Step 4 is to determine host-rock characteristics for the new sill: define the weight fraction that converts to methane w CH4 and the kerogen composition parameter K using the cumulative probability distributions in Fig. 4. Step 5 is to determine thermogenic carbon generation parameters for the new sill: use the new thermogenic carbon emissions parameterisation to estimate the final mass of generated carbon M therm (A, K, S, Z, w CH4 , β), the emissions decay time τ decay (S, K, Z) and the emissions decay exponent p(S, K, Z) from the sill dimensions and host-rock characteristics obtained in steps 3 and 4. Step 6 is to determine new sill location, if required: determine polar coordinates based on the plume head model31 and known sill province geography, determine aureole overlaps, and correct M therm to avoid double-counting. Step 7 is to determine direct magma degassing carbon emissions parameters for the new sill: the final mass of degassed carbon M magma (A, S) and the emissions decay time τ solid (S, Z) are estimated from the sill dimensions obtained in step 3. Step 8 is to determine carbon emissions flux and cumulative emitted carbon from the entire sill province by summing thermogenic emissions (q, m)(M therm , τ decay , p, t) and magma degassing emissions (q, m)(M magma , τ solid , t) from all sill–vent systems that exist at this time step. Finally, repeat steps 1–8 as required: for fixed τ repeat , the total number of sills is specified a priori; for variable τ repeat , the calculation is continued until τ repeat goes to infinity, thereby determining the total number of sills.

Combining carbon emissions from the entire sill province

We calculated carbon emissions flux and cumulative emitted carbon from the entire sill province by summing thermogenic and magma degassing emissions from many sill–vent systems (step 8 of the stochastic procedure). Say that at time t after sill province initiation, n sills have been intruded. The initiation times for these sills have already been determined sill-by-sill and are represented as t 0 (i), i = 1, 2… n. Time relative to intrusion of each sill is denoted as

$$\begin{array}{*{20}{c}} {\Delta t\left( i \right) = t - t_0\left( i \right)} \end{array}.$$ (31)

Emissions flux from the province is

$$\begin{array}{*{20}{c}} {Q_{{\mathrm{prov}}}\left( t \right) = \mathop {\sum }\limits_{i = 1}^n \left( {q_{{\mathrm{therm}}}\left( {\Delta t\left( i \right)} \right) + q_{{\mathrm{magma}}}\left( {\Delta t\left( i \right)} \right)} \right)} \end{array}$$ (32)

and the cumulative mass is

$$\begin{array}{*{20}{c}} {M_{{\mathrm{prov}}}\left( t \right) = \mathop {\sum }\limits_{i = 1}^n \left( {m_{{\mathrm{therm}}}\left( {\Delta t\left( i \right)} \right) + m_{{\mathrm{magma}}}\left( {\Delta t\left( i \right)} \right)} \right)} \end{array}.$$ (33)

The carbon isotopic composition of the emissions is

$$\begin{array}{*{20}{c}} {\delta ^{13}C_{{\mathrm{prov}}}\left( t \right) = \frac{1}{{Q_{{\mathrm{prov}}}\left( t \right)}}\mathop {\sum }\limits_{i = 1}^n \left( {\delta ^{13}C_{{\mathrm{therm}}}\left( {\Delta t\left( i \right)} \right)q_{{\mathrm{therm}}}\left( {\Delta t\left( i \right)} \right) + \delta ^{13}C_{{\mathrm{magma}}}q_{{\mathrm{magma}}}\left( {\Delta t\left( i \right)} \right)} \right)} \end{array}.$$ (34)

Isotopic composition of carbon emissions

The isotopic composition of mantle-derived carbon is –7‰. The isotopic composition of thermogenic carbon depends on the composition and also the temperature of maturation of organic material in the host rock52. Both of these relationships are observed in the data compilation in Supplementary Fig. 4a. We used these data to define the relationships between δ13C and vitrinite reflectance, R o , a proxy for maturation temperature, for refractory (gas-prone) kerogen

$$\begin{array}{*{20}{c}} {\delta _{{\mathrm{ref}}}^{13}C = 15.4\log R_{\mathrm{o}} - 41.3} \end{array}$$ (35)

and labile (oil-prone) kerogen

$$\begin{array}{*{20}{c}} {\delta _{{\mathrm{lab}}}^{13}C = 22.0\log R_{\mathrm{o}} - 32.0} \end{array}.$$ (36)

We then used the coupled thermal-kinetic reaction model to determine the relationship between reaction temperature and vitrinite reflectance (Supplementary Fig. 4b). Combining the δ13C(R o ) and R o (T) relationships gives the carbon isotopic composition as a function of reaction temperature δ13C(T) (Supplementary Fig. 4c). The temperature of each reaction front can be tracked during numerical thermogenic reaction modelling, allowing us to estimate the carbon isotopic composition of thermogenic methane from both kerogen types over time for any sill (Supplementary Fig. 4d). These time-series show two phases. Between sill intrusion and full solidification, the carbon isotopic composition is relatively heavy and decreases little through time. At the time of full solidification, there is a kink in the composition curve, and thereafter the isotopic composition gets lighter over time. This behaviour arises because the temperature of the reaction front (where most of the methane flux is generated) decreases slowly until the sill solidifies fully and more rapidly thereafter until the background temperature is regained.

Estimating sill intrusion repeat time from mantle plume flux

There is as yet no effective method to determine the rate of LIP carbon emissions on a sub-thousand-year timeframe. One important difficulty is that it is impossible to measure the period τ repeat between intrusion of successive sills directly using radiometric or biostratigraphic dating. Here we derive an alternative constraint on the sill intrusion repeat time τ repeat based on the mantle convection process that generated the NAIP. If n sills are intruded in time t, then τ repeat and the corresponding frequency f are

$$\begin{array}{*{20}{c}} {f = \frac{1}{{\tau _{{\mathrm{repeat}}}}} = \frac{{{\mathrm{{d}}}n}}{{{\mathrm{{d}}}t}}} \end{array}.$$ (37)

To link τ repeat to mantle plume flux, first expand this expression in terms of the surface area of the sill province A LIP (t) to give

$$\begin{array}{*{20}{c}} {f = \frac{{{\mathrm{{d}}}n}}{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}\frac{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}{{{\mathrm{{d}}}t}}} \end{array}.$$ (38)

We define the sill area-density (i.e. the number of sills per unit area) as

$$\begin{array}{*{20}{c}} {\rho = \frac{{{\mathrm{{d}}}n}}{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}} \end{array}.$$ (39)

If expansion of the NAIP sill province directly reflects sub-plate expansion of the hot mantle that generated the sill magma, then

$$\begin{array}{*{20}{c}} {Q = \frac{{{\mathrm{{d}}}A_{{\mathrm{mantle}}}}}{{{\mathrm{{d}}}t}} = \frac{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}{{{\mathrm{{d}}}t}}} \end{array},$$ (40)

where Q is the mantle plume area flux, averaged vertically across the plume head of surface area A mantle . In this case, the expression for intrusion frequency becomes

$$\begin{array}{*{20}{c}} {f = \frac{1}{{\tau _{{\mathrm{repeat}}}}} = \rho Q} \end{array}.$$ (41)

More generally, expansion of the sill province is smaller than Q because thick lithosphere inhibits decompression melting, so that A LIP < A mantle . This effect can be incorporated by expanding in terms of plume head area to give

$$\begin{array}{*{20}{c}} {f = \frac{1}{{\tau _{{\mathrm{repeat}}}}} = \rho \frac{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}{{{\mathrm{{d}}}A_{{\mathrm{mantle}}}}}\frac{{{\mathrm{{d}}}A_{{\mathrm{mantle}}}}}{{{\mathrm{{d}}}t}}} \end{array}.$$ (42)

A model of the plume head is required to link dA mantle /dt to Q. We used a simple kinematic model that has been previously applied to model Cenozoic development of the Iceland Plume38, including estimating Paleocene–Eocene plume flux from dynamic support histories derived from sedimentary successions side of Scotland37,38,43, and estimating Neogene plume flux from the V-Shaped Ridges of oceanic crust south of Iceland34,35. In this model, plume mantle flows radially outwards from a point source between two rigid parallel plates that represent the top and base of the asthenosphere. The original model for a radially symmetrical plume head is

$$\begin{array}{*{20}{c}} {A_{{\mathrm{mantle}}} = \pi R^2 = 1.5Qt} \end{array},$$ (43)

where R is the radius of the plume head. We assume that the marker for the edge of the plume head is the location of maximum temperature anomaly, which is expected to be the location of peak melting, and to lie directly beneath the location of peak dynamic support of the overlying plate. The factor 1.5 occurs because this plume head marker travels at the maximum speed within the plume head channel, which is 1.5 times greater than the mean speed across the channel for a Poiseuille velocity profile38. If the plume head is elliptical, then the central term of Eq. 43 becomes πγR2, where γ is the aspect ratio (i.e. the ratio of short and long radii) and R is now the long radius of the ellipse. Differentiating gives

$$\begin{array}{*{20}{c}} {\frac{{{\mathrm{{d}}}A_{{\mathrm{mantle}}}}}{{{\mathrm{{d}}}t}} = 1.5Q} \end{array}.$$ (44)

If we define a function

$$\begin{array}{*{20}{c}} {\Gamma = \frac{3}{2}\frac{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}{{{\mathrm{{d}}}A_{{\mathrm{mantle}}}}} = \frac{3}{{4\pi \gamma R}}\frac{{{\mathrm{{d}}}A_{{\mathrm{LIP}}}}}{{{\mathrm{{d}}}R}}} \end{array}$$ (45)

to relate the areal expansion of the plume head and the sill province, then the expression for the sill intrusion time period becomes

$$\begin{array}{*{20}{c}} {f = \frac{1}{{\tau _{{\mathrm{repeat}}}}} = \rho \left( R \right)\Gamma \left( R \right)Q} \end{array}.$$ (46)

Relating sill distribution and carbon emissions

Although τ repeat , the time period between emplacement of successive sills, is a fundamental parameter that controls gas emissions flux Q prov , τ repeat is too small to measure directly using the geological record. It is therefore useful to derive an expression that relates Q prov to the spatial distribution of sills ρ, which can be measured directly. Results of stochastic modelling with constant sill recurrence time indicate an inverse relationship between τ repeat and Q prov at long-time steady state (Eq. 1). The derivation in the previous section relates τ repeat to ρ (Eq. 46). Combining these two relationships eliminates τ repeat :

$$\begin{array}{*{20}{c}} {\rho = \frac{{Q_{{\mathrm{prov}}}}}{{1.65\,{\mathrm{\Gamma }}\,{\mathrm{Q}}}}} \end{array}.$$ (47)

This expression can be used to investigate whether the ρ observations are compatible with Q prov estimates obtained independently from sedimentary records of environmental change. Figure 7c shows the ρ dataset overlain by a Gaussian curve that was used to calculate the NAIP emissions histories in Fig. 8. Figure 7c also shows an alternative ρ curve, calculated using Eq. 47 with Q prov from a sedimentary record5, the median value for Q, and Γ from Fig. 7e. The two ρ curves are quite different near the sill province centre, but they are both compatible with the data since ρ measurements near the province centre are sparse. Supplementary Fig. 5 shows NAIP emissions flux histories calculated using both ρ curves. Both a smooth rise to a lower peak flux and also a more irregular rise to a higher peak flux are probably compatible with available ρ data at present, and it should be possible to distinguish between these two scenarios in future with more ρ measurements near the province centre.

The NAIP Q prov model based on the alternative ρ model is often lower than the sediment-based Q prov curve from which the alternative ρ model was derived via Eq. 47 (Supplementary Fig. 5, purple versus blue). This discrepancy is expected for two reasons. First, Eq. 47 assumes long-time steady state via Eq. 1, so it should underestimate ρ when Q prov increases rapidly. Secondly, Eq. 47 does not account for overlap between sill aureoles, so it should underestimate ρ at later times. Nevertheless, Eq. 47 provides a useful new starting point for estimating the extent to which uncertainties in ρ and also Γ might explain differences between carbon source- and carbon sink-derived estimates of Q prov .

Determining sill locations

Emission histories can be calculated without determining sill locations, and this procedure was used in the constant-τ repeat calculations (Fig. 5). The constant-τ repeat simulations represent maximum possible long-time emissions because they implicitly assume that each new sill intrudes host rock capable of generating methane. There are two reasons for assigning sill locations. First, the resulting maps provide a more intuitive description of our proposed link between mantle plume head and sill intrusion processes (Fig. 8g–j; Supplementary Movie 1). Secondly, sill locations can be used to ensure that when a sill intrudes host rock that has already been thermally matured by a previous sill, the emissions are not counted twice in the combined emissions histories.

To define a sill location, two of the random numbers that define the sill in step 2 of the stochastic modelling algorithm are used to choose the sill location radius and the sill location azimuth. Since melt generation is proportional to mantle temperature, the sill location radius was chosen based on the location of the mantle thermal anomaly specified by the plume head model. Equation B14 of ref. 38 gives the temperature structure within the plume head channel as a function of plume head area, depth and time:

$$\begin{array}{*{20}{c}} {T\left( {A_{{\mathrm{mantle}}},z,t} \right) = \frac{1}{{Q\delta \sqrt {2\pi } }}{\mathrm{{exp}}}\left[ {\frac{{ - \left( {A_{{\mathrm{mantle}}}/Qf_{\mathrm{z}} \, - \, t} \right)^2}}{{2\delta ^2}}} \right]} \end{array},$$ (48)

where δ is the standard deviation of the pulse of anomalously hot mantle and f z is a function that describes vertical variation in velocity across the plume head. If time is scaled so that t = t′δ, plume head area is scaled so that A plume = A′ Qδ and f z is scaled by the thickness of the plume head channel, this equation becomes

$$\begin{array}{*{20}{c}} {T\left( {{A^{{\prime}}},{z^{{\prime}}},{t^{{\prime}}}} \right) = \frac{1}{{\left( {A_{{\mathrm{plume}}}/{A^{\prime}} } \right)\sqrt {2\pi } }}{\mathrm{{exp}}}\left[ {\frac{{ - \left( {{A^{\prime}} /{f^{\prime}} _{\mathrm{z}} \, - \, {t^{\prime}} } \right)^2}}{2}} \right]} \end{array},$$ (49)

where \({f^{{\prime}}}_{\mathrm{z}} = 1.5\left( {1 - {z^{{\prime}}}^2} \right)\) for a Poiseuille flow profile38. This equation can be integrated numerically over depth to give the mean temperature anomaly as function of radius and time (Supplementary Fig. 6a). If this mean temperature function is integrated over radius and suitably scaled, the result can be interpreted as a cumulative probability distribution for the sill location radius (Supplementary Fig. 6b). In practice, we used area rather than radius for the horizontal coordinate and subtracted the position of the plume head marker (taken as the location of maximum temperature anomaly) in order to form a cumulative probability function for the sill location area anomaly δA′ for each sill (Supplementary Fig. 6d). This function was used as a look-up table to determine the sill location radius from the appropriate random number for each new sill. To create Figs. 8g–j and Supplementary Movie 1, we used Q = 4 km2 yr–1 (Fig. 6) and δ = 40 kyr (ref. 38).

After assigning a sill radius, an azimuth was selected from the ranges that lay within the deep basin chain (Fig. 6a), assuming a uniform distribution. Geographical coordinates were calculated from the (radius, azimuth) pair assuming the same elliptical plume head location and geometry used to determine mantle plume flux.

Determining overlap between sill aureoles

No new thermogenic methane is generated where the aureole of a newly intruding sill overlaps an existing sill or its aureole. This effect was modelled by modifying the M therm values of each sill. Each new sill was compared with existing sills in turn. The distance d between the centres of each sill pair was determined from their radial coordinates. The sill footprints overlap partially if d < (R 1 + R 2 ) and fully if d < |R 1 − R 2 |, where the subscripts distinguish the two sills. For each partially overlapping sill pair, the area of the lens-shaped overlap is

$$\begin{array}{ccccc}\\ A_{\mathrm{O}} = & \hskip -8.3pcR_1^2{\mathrm{cos}}^{ - 1}\left( {\frac{{d^2 + R_1^2 - R_2^2}}{{2dR_1}}} \right) + R_1^2{\mathrm{cos}}^{ - 1}\left( {\frac{{d^2 + R_2^2 - R_1^2}}{{2dR_2}}} \right)\\ \\ & \hskip -0.5pc- 0.5\sqrt {\left( {R_1 + R_2 - d} \right)\left( {R_1 - R_2 + d} \right)\left( {R_2 - R_1 + d} \right)\left( {R_1 + R_2 + d} \right)} \\ \end{array}.$$ (50)

The depths of the top and base of each sill are Z ± 0.5 S, and the depths of the top and base of the aureole were approximated by Z ± 0.5 S (Y + 1). These depths were used to test for vertical overlap between aureoles analogous to the test for overlapping footprints. The overlapped aureole volume is A O Z O , where Z O is the vertical overlap. The overlapped volume as a proportion of the total aureole volume is

$$\begin{array}{*{20}{c}} {\varphi = \frac{{A_{\mathrm{O}}Z_{\mathrm{O}}}}{{ASY}}} \end{array}.$$ (51)

This fraction was used to update the total carbon mass parameter of each sill to \(\left( {1 - \varphi } \right)M_{{\mathrm{therm}}}\).

Two factors complicate this basic procedure. If a new sill overlaps multiple existing sills, it becomes awkward to determine φ because the multiple overlapping volumes may overlap each other. It is not worth solving this problem in full bearing in mind the various simplifying assumptions already made, including circular sills, symmetrical aureoles and uniform initial geotherm. Instead, we note the probability that a point lying within a large volume V is overlapped by at least one of k small independent blobs each of volume v is \(1 - \left( {1 - \frac{v}{V}} \right)^k\). Hence an estimate of the factor to correct M therm for k patches of aureole with fractional overlaps φ 1…k is

$$\begin{array}{*{20}{c}} {\varphi ^ \ast = 1 - \mathop {\prod }\limits_{i = 1}^k \left( {1 - \varphi _i} \right)} \end{array}.$$ (52)