Geodynamic modelling framework

All 2D and 3D geodynamic simulations have been conducted with the codes M2TRI and M3TET (refs 18, 41), respectively. These codes are written in MATLAB and solve for viscous incompressible flow, advection-diffusion of temperature, advection of composition and melting processes of mantle rocks that reach their pressure-, composition- and volatile-dependent solidus temperature.

Governing equations and numerical methods

Mantle flow is described as viscous creep of an incompressible fluid (‘Stokes flow’) using the Boussinesq approximation, that is, density variations are only considered in the buoyancy term. We solve the equations for conservation of mass

and momentum

using the viscous stress–strain rate relationship

See Supplementary Table 4 for a complete list of variables and parameters used in the calculations. We assume a Newtonian rheology with the dynamic viscosity η depending on temperature T (in units of Kelvin), pressure p and weight fraction χH 2 O of water in the mantle rocks:

The pre-exponential factor A controls the viscosity increase during dehydration:

where A max is the maximum increase in viscosity after complete dehydration. For numerical reasons we limit the viscosity variations to be between 3 × 1018 Pa s and 5 × 1022 Pa s. Density variations depend on temperature changes and the cumulative degree of melting F:

where α is the thermal expansion coefficient and β parameterizes the compositional buoyancy resulting from melting.

Energy conservation is formulated as

where ΔH=ΔS·T is the latent heat of fusion, and ΔS is the entropy of fusion.

We use the finite element method to solve for viscous flow and heat conduction. Advection of temperature and compositional fields are handled by the method of characteristics (also called semi-Lagrange method; for example, ref. 42). The matrix equation for the 2D Stokes flow problem is derived using the consistent penalty method combined with so-called Powell–Hestenes iterations (a formulation similar to that used in the MILAMIN code43). All matrix equations are solved using the Cholesky direct solver of the numerical library SuiteSparse44. In 3D, where the larger number of unknowns requires iterative solvers, we use a Schur-complement formulation45 to solve the velocity-pressure problem. This leads to an outer iteration to solve for pressure and an inner iteration updating the velocity field. For both iterations we use Conjugate Gradient (CG) solvers46, preconditioned by the pressure mass matrix (outer CG) and a single V-cycle of geometric multigrid (inner CG). M3TET is parallelized using MATLAB’s Parallel Computing Toolbox in combination with the MUTILS47 tools. For more details on the solution techniques see refs 18, 41.

The formulation to calculate melting rates in the mantle is based on ref. 48 and the reader is referred to this reference for the details of the melting model. We define solidus functions TS of the form

where is the solidus temperature at surface pressure, is the solidus dependence on pressure and is the solidus dependence on depletion (Supplementary Table 4). We further consider wet melting and have implemented the parameterization described in ref. 49 into our framework by modifying the solidus-depletion gradient depending on the rock’s current water content. Water is treated like an incompatible chemical element that preferentially enters the melt phase. We use a partition coefficient D H 2 O =0.01 (ref. 49) and assume fractional melting, where the water content XH 2 O of a mantle rock is calculated from its initial water content and cumulative degree of melting F (ref. 50)

The solidus-depletion gradient is further modified to mimic the drop in melt productivity once clinopyroxene is exhausted from the residue (for example, ref. 51). For simplicity, we assume that this is the case for a degree of melting F>0.16 and increase by a factor of 10 to reduce productivity (Supplementary Fig. 4).

2D mid-ocean ridge models

The numerical domain is 200 km deep and 400 km wide, and we assume symmetry along the ridge axis, that is, we model one half of the ridge. Temperature boundary conditions are 0 °C at the top and a potential mantle temperature T M =1,335 °C at the bottom inflow boundary. The symmetry axis (left side) and outflow boundary (right side) are thermally insulating. A constant half-spreading rate is prescribed along the top. Velocities along the inflow and outflow boundary are solved for in the first time step, with all buoyancy forces set to zero to avoid convective instabilities. From the second time step on velocities along the in- and outflow boundaries are fixed to these values and buoyancy forces are considered. We prefer this method for deriving the velocity boundary conditions over an analytical corner flow solution because the calculated in- and outflow velocities reflect the viscosity structure of the uppermost mantle. Calculations start with an initial temperature field calculated from the GDH1 plate model52 and we ran all calculations until a steady state in melt production was reached. Mantle composition and potential mantle temperature have been calibrated to produce a uniform crustal thickness of about 7 km at full spreading rates faster than 50 mm per year, in agreement with seismic observations53. We conducted a total of 360 2D numerical calculations in which we have varied (i) the half-spreading rate v HS between 1 and 100 mm per year, (ii) the initial water content of the upper mantle (50 p.p.m. w , 100 p.p.m. w and 200 p.p.m. w ; see ref. 54), and (iii) the maximum increase of the mantle rock viscosity during melting-related dehydration A max (factor increase of 5, 10, 50 and 100). The value of the latter is somewhat uncertain so that we explored a broad range. The consequences of different parameter choices for the predicted global magma and CO 2 mobilization rates are summarized in Supplementary Table 5. The global MOR model presented in the main text and in Fig. 4 is based on 2D simulations that assume 100 p.p.m. w water in the mantle source and a viscosity increase of factor 50 during melting-induced dehydration.

3D plume geodynamic models

Modelling the ascent and melting of a mantle plume beneath a moving lithospheric plate requires a 3D geodynamic model. We assume symmetry across a vertical plane parallel to the direction of plate motion and through the center of the plume, that is, we model one half of the plume. The numerical domain extends 400 km into the direction perpendicular to plate motion, 1,000 km into the direction of plate motion, and to 400 km depth. To have a precise control on the thickness of the overriding lithosphere we have excluded the rigid lithospheric plate from the model domain. Instead, the top of the domain corresponds to the base of the lithosphere and plume upwelling stops when this boundary is reached. Accordingly, we use a plate age-dependent heat flux boundary condition at the top to mimic the heat flow into the overlying lithosphere. A constant plate speed is prescribed along the top. A potential mantle temperature T M =1,335 °C is prescribed at the inflow boundary (left wall) upstream of the plume, all other side boundaries are insulating. The horizontal inflow velocities through the left wall are important because they create the ‘mantle wind’ through which the plume ascends. Here we initialize a viscosity-dependent vertical Couette flow profile with plate speed at the top and zero-velocity at 400 km depth. Free horizontal outflow is allowed through the right boundary downstream of the plume, while front and back boundaries are symmetric. Along the bottom we prescribe T M =1,335 °C except for a half-circular plume inflow region with 100 km radius, where we define a Gaussian-shaped temperature anomaly and a parabolic inflow profile. We control the plume strength by prescribing its excess temperature and buoyancy flux. Inflow velocities are then adjusted such that the desired buoyancy flux is exactly matched for the plume’s thermal structure. Depending on buoyancy and rheology the plume self-adjusts its radius well before reaching the melting region. All calculations ran until the plume head had passed the melting zone and a steady melt production established. We used the best-studied mantle plume/hotspot system Hawaii to calibrate the mantle composition for the mantle plume parameter study: For a plume buoyancy flux Q B =8,700 kg s−1 (ref. 55) and an excess temperature T exc =290 °C (ref. 56) we aimed to achieve an onset of melting at ∼170 km depth consistent with the seismically inferred melt-rich zone above this depth57. The total melt production below to the ∼91 Myr ago old lithosphere58 of ∼84 km thickness52 should be about 10 m3 s–1 (ref. 59). Meeting these constraints required a more fertile and wetter composition for mantle plume material compared to the ‘upper mantle’ composition . We conducted a 3D parameter study in which we systematically varied four key parameters related to the mantle plume and the plate tectonic setting of the associated hotspot: (i) plume buoyancy flux Q B (500, 1,000, 4,000, 6,000 and 10,000 kg s−1), (ii) plume excess temperature T exc (100, 200 and 300 °C), (iii) lithosphere thickness h L at the hotspot location (50, 70 and 90 km), and (iv) plate speed v L relative to the hotspot (10, 40 and 80 mm per year). The resulting 126 3D calculations span a 4D parameter space (Supplementary Fig. 5), within which the 43 oceanic hotspots (see Supplementary Table 2) have been linearly interpolated to derive their melt production and CO 2 release rates. The advantage of the parameter-space approach is that it allows us to tests how critically the predicted global magma fluxes and CO 2 mobilization rates depend on the four parameters that characterize each mantle plume/hotspot. The results of these sensitivity tests with the global plume model are summarized in Supplementary Table 6.

CO 2 mobilization and global mid-ocean ridge melting model

Assuming fractional melting50, the rate of CO 2 mobilization into the melt phase in each model calculation was calculated using an initial concentration of 140 p.p.m. w CO 2 (ref. 19) and a partition coefficient of D CO 2 =0.01. Note the latter is a very moderate estimate since CO 2 is likely to behave more incompatible than water, for which we use the same value. Each numerical calculation provides a rate of CO 2 mobilization into the melt phase per metre ridge axis for a particular spreading rate. The estimate for the total CO 2 mobilization along the global MOR system is derived using the global distribution of spreading rates22. Note that our global estimates (Fig. 4 and Supplementary Table 3) are in very good agreement with independent estimates for global magma production and CO 2 release at MORs19. Melting beneath MOR results primarily from the decompression of mantle rocks. The enhanced CO 2 mobilization into the melt phase at times of dropping sea level is calculated by superimposing the additional sea level-induced decompression onto the decompression of mantle rocks during upwelling. For a given amplitude Δh SL in sea level drop over a time span Δt SL and a sea water density ρ w =1,030 kg m−3 we calculate the rate of pressure change in the underlying mantle, which is then considered in the melting formulation.

CO 2 mobilization and global mantle plume melting model

For the global hotspot melting model we use a similar approach. The steady state solution of each 3D calculation is used to calculate the baseline melt production and CO 2 mobilization into the melt phase, while the superimposed decompression during a prescribed sea level drop gives the enhanced values during ice sheet growth.

In contrast to the submarine MOR, oceanic hotspot melting zones are located below large volcanic edifices that are partially subaerial. A falling sea level, however, will only change the hydrostatic pressure on the seafloor. To quantify the amplitude of the sea level-induced pressure signal below the island we have calculated the viscoelastic response of an oceanic lithosphere carrying an island using the free surface and viscoelastic capabilities of M3TET. We initialized cone-shaped islands of different sizes, applied a water depth-dependent normal force boundary condition at the top and calculated the island’s subsidence until an isostatic equilibrium established (cf. Fig. 3). Then we simulated the falling sea level around the island by reducing the normal forces and monitored the pressure change in the viscous asthenospheric mantle. Elastic parameter choices are summarized in Supplementary Table 4. The pressure drop increases from a minimum value central below the island to the full sea level-induced signal at some distance to the island (Fig. 3). We have used these distance-to-island dependent pressure signals to calculate the enhanced melt production and CO 2 mobilization rates in the plume melting region. For simplicity we have grouped all oceanic hotspots into four island-categories depending on their subaerial area (see Supplementary Table 2).

The derivation of global hotspot-related melt production and CO 2 mobilization rates is more complicated than for MOR because more parameters control the melting of a particular hotspot. The first parameter is the plume buoyancy flux, where we use the values given in ref. 55 and in ref. 60. The second parameter is the plume excess temperature, which is based on ref. 56. Unfortunately, not every plume with estimated buoyancy flux also has an estimate for excess temperature and vice versa. To fill these data gaps we calculate the missing value using a linear interpolation between Hawaiian-like values (8,700 kg s−1 and 290 °C) representing the ‘strong plume’ case and values that are typical for a weak plume (500 kg s−1 and 150 °C). In case that both excess temperature and buoyancy flux are unconstrained we assume a weak plume with 500 kg s−1 and 130 °C. The 3rd and 4th parameters are speed and age of the lithosphere at the hotspot location, respectively. We use the hotspot locations given in ref. 60 and the NUVEL-1A model61,62 to first identify the lithospheric plate hosting the hotspot. The plate motion relative to the hotspot is then calculated using the global plate motion model described in ref. 63. Lithosphere thickness is calculated from the age of the lithosphere at the hotspot location (using the data set in ref. 58) and then converted to plate thickness using the GDH1 model52 for the thermal evolution of oceanic lithosphere. We define the lithosphere-asthenosphere boundary at the depth where 90% of the basal temperature in the GDH1 model is reached. Following the argumentation in ref. 64 we assume a compositional lithosphere of 50 km thickness that forms during melting and dehydration of mantle rocks at mid-ocean ridges. Hence, we assume a minimum plate thickness of 50 km even if the GDH1 model predicts a thinner thermal lithosphere. Using these four parameters, each oceanic hotspot/plume (Supplementary Table 2) is located within the 4D parameter space spanned by the 3D calculations (Supplementary Fig. 5), and melt production and CO 2 mobilization are linearly interpolated. This approach had the advantage that we were able to test different estimates for plume buoyancy fluxes given, for example, in refs 55, 65 and to test how critically the predicted global fluxes depend on each of the four parameters. The results of these tests, summarized in Supplementary Table 6, show that the baseline fluxes are especially sensitive to variations in plume excess temperature and lithosphere thickness. The relative increase in global magma and CO 2 flux during sea level drop are, nonetheless, very robust and do not significantly change.

Carbon cycle simulations

The Box model of the Isotopic Carbon CYCLE–BICYCLE—is used to assess the impact of volcanic emissions on atmospheric CO 2 concentrations. BICYCLE has been widely used for a broad range of Pleistocene paleoclimate applications29. BICYCLE is a model of the atmosphere–ocean–terrestrial biosphere-part of the global carbon cycle, in which the response of the deep ocean sediments (carbonate compensation feedback) to any change in the carbon cycle, including an external CO 2 flux, is estimated with a pulse response function. It contains one homogeneous mixed atmosphere box, ten boxes in the ocean and seven boxes on land. The ocean is separated in surface, intermediate and deep waters in the Atlantic, Southern Ocean and Indo-Pacific. We used the BICYCLE model to assess the variations in atmospheric CO 2 resulting from injecting volcanic CO 2 into the model, either subaerial directly into the atmosphere or submarine into the water column. For these volcanic CO 2 release experiments anomalies in atmospheric CO 2 are calculated based on the model set-up as previously published29. Carbon cycle simulations are run for 7 kyr into quasi–equilibrium for climate conditions found 83 kyr ago leading to a simulated atmospheric CO 2 concentration of 227 p.p.m. v for a set-up of atmosphere–ocean including carbonate compensation feedback. The volcanic CO 2 release scenarios will be analysed as anomalies in this setup.

In detail, and based on the geodynamic modelling results, four carbon cycle scenarios have been simulated. In these simulations all boundary conditions are fixed at values for 83 kyr ago. In scenario S1 a sea level drop of 60 m over 5 kyr releases a total of 400 Gt CO 2 (188 Gt CO 2 from MOR+212 Gt CO 2 from hotspots) corresponding to an annual flux of 80 Tg CO 2 . In S2 a sea level drop of 60 m over 15 kyr releases a total of 446 Gt CO 2 (188 CO 2 from MOR+258 Gt CO 2 from hotspots) corresponding to an annual flux of 30 Tg CO 2 . In S3 a sea level drop of 80 m over 15 kyr releases a total of 601 Gt CO 2 (252 CO 2 from MOR+349 Gt CO 2 from hotspots) corresponding to an annual flux of 40 Tg CO 2 . In scenario S4 a sea level drop of 100 m over 10 kyr releases a total of 697 Gt CO 2 (314+383 Gt CO 2 from hotspots) corresponding to an annual flux of 70 Tg CO 2 .

Code availability

The M3TET and M2TRI codes are co-owned by Hamburg University, GEOMAR as well as Royal Holloway and the source code is not freely distributed. To help with all questions concerning reproducibility of the numerical results, further details on the used numerical implementations as well as data analysis procedures, additional data files and possibly precompiled sub-functions can be requested from the first author.

Data availability

All relevant geodynamic modelling data are available from the first author. The carbon cycle modelling data are available at Pangaea with the doi: 10.1594/PANGAEA.874815.