Consequences of glacial cycles for magmatism and carbon transport at mid-ocean ridges

TL;DR

This study models how glacial cycles influence magmatism and carbon transport at mid-ocean ridges, revealing that sea-level changes cause significant, lagged fluctuations in volcanic carbon emissions, with implications for climate feedback mechanisms.

Contribution

The paper introduces a new model for mid-ocean ridge response to sea-level changes, improving upon previous assumptions and incorporating thermodynamic effects of mantle carbon.

Findings

Magma and carbon flux fluctuate by up to 20% and 10%, respectively, over glacial cycles.

Peak emissions lag sea-level peaks by less than 20 kyr, with lag times sensitive to melt segregation rates.

Fluctuations are primarily controlled by changes in melting rate, not melt transport times.

Abstract

Magmatism and volcanism transfer carbon from the solid Earth into the climate system. This transfer may be modulated by the glacial/interglacial cycling of water between oceans and continental ice sheets, which alters the surface loading of the solid Earth. The consequent volcanic-carbon fluctuations have been proposed as a pacing mechanism for Pleistocene glacial cycles. This mechanism is dependant on the amplitude and lag of the mid-ocean ridge response to sea-level changes. Here we develop and analyse a new model for that response, eliminating some questionable assumptions made in previous work. Our model calculates the carbon flux, accounting for the thermodynamic effect of mantle carbon: reduction of the solidus temperature and a deeper onset of melting. We analyse models forced by idealised, periodic sea level and conclude that fluctuations in melting rate are the prime control on magma and carbon flux. We also discuss a model forced by a reconstruction of eustatic sea level over the past 800 kyr. It indicates that peak-to-trough variations of magma and carbon flux are up to about 20% and 10% of the mean flux, respectively. Peaks in mid-ocean ridge emissions lag peaks in sea-level forcing by less than about 20 kyr and the lag could well be shorter. The amplitude and lag are sensitive to the rate of melt segregation. The lag is much shorter than the time it takes for melt to travel vertically across the melting region.

Tables

Table 1. Table 1: Physical parameters and corresponding dimensionless model parameters. Values as used in calculations unless otherwise stated. See Appendix  A for justification.

Parameters	Value	Unit	Description
${\rho }_w$	1000	kg/m3	Sea-water density
$\rho$	3300	kg/m3	Mantle density
$W_0$	2	cm/yr	Mantle upwelling velocity
$H$	130	km	Height of melting column
$H_{\text{dry}}$	65	km	Height of dry melting column
$\Delta S$	0.1	km	Peak-to-trough amplitude of sea-level fluctuation
$F_{\text{max}}$	0.2		Maximum degree of melting
$D_c$	${10}^{-4}$		Partition coefficient of carbon
$n$	2		Exponent in permeability-porosity relationship
$𝒬$	${10}^5$		Liquid flux scale
${\Gamma }^{\ast }$	$0.4$		Melting rate scale
$ℳ$	$0.2$		Effect of carbon on the mantle solidus scale
${\varphi }_{\max }$	$0.0014$		Maximum steady-state porosity at top of column
$w_0$	$2.8$	m/yr	Maximum melt velocity at top of column
$t_{\mathrm{dry}}$	$46$	kyr	Melt transport time across dry melting region
$\alpha$	30	∘	Dip of decompaction channel

Equations

$$\varphi w+(1-\varphi )W=1,$$

$$\Gamma =\frac{\partial \varphi }{\partial t}-\frac{\partial }{\partial z}\left[(1-\varphi )W\right].$$

$$\varphi (w-W)W_0$$

$$\Rightarrow 1-W$$

$$\mathcal{Q}\equiv \frac{\Delta \rho gk}{\mu W_0}$$

$$T=T_0+{\gamma }^{-1}(P-P_0)+M{c}_0(1-c_s),$$

$$P=P_0-\rho gHz+{\rho }_{w}gH(S-S_0).$$

$$z=\delta (t)\equiv \frac{{\rho }_w}{\rho }(S-S_0).$$

$$\frac{\mathrm{D}T}{\mathrm{D}t}=-\frac{L}{c_p}\Gamma ,$$

$$\Gamma =\underset{{\Gamma }_0}{\underbrace{{\Gamma }^{\ast }}}\underset{{\Gamma }_P}{\underbrace{-{\Gamma }^{\ast }\frac{d\delta }{dt}}}\underset{{\Gamma }_c}{\underbrace{+\mathcal{M}\frac{\mathrm{D}{c}_s}{\mathrm{D}t}}},$$

$${\Gamma }^{\ast }\equiv \frac{\rho gH{c}_p}{\gamma L},\mathcal{M}\equiv \frac{c_{p}M{c}_0}{L}.$$

$${\Gamma }_P\equiv -{\Gamma }^{\ast }\frac{d\delta }{dt}=-{\Gamma }^{\ast }\frac{d}{dt}\left[\frac{{\rho }_w}{\rho }\left(S-S_0\right)\right].$$

$$\mathcal{M}={\Gamma }^{\ast }\left(1-\frac{H_{\mathrm{dry}}}{H}\right).$$

$$F_{\max }={\Gamma }^{\ast }-\mathcal{M}.$$

$$\frac{\partial }{\partial t}\left[\varphi c_l+(1-\varphi )c_s\right]+\frac{\partial }{\partial z}\left[\varphi w{c}_l+(1-\varphi )W{c}_s\right]=0,$$

$$c_s=D_c{c}_l,$$

$$\varphi =0,c_s=1(\text{at }z=\delta ).$$

$$\frac{\partial \varphi }{\partial t}+\frac{\partial Q}{\partial z}={\Gamma }^{\ast }-{\Gamma }^{\ast }\frac{d\delta }{dt}+\frac{\mathrm{D}c}{\mathrm{D}t},$$

$$D\frac{\mathrm{D}c}{\mathrm{D}t}+\frac{\partial }{\partial t}\left(c\varphi \right)+\frac{\partial }{\partial z}\left(cQ\right)=0,$$

$$c\equiv \mathcal{M}{c}_s,D\equiv D_c/(1-D_c).$$

$$Q(\varphi )=\mathcal{Q}{\varphi }^n(1-\varphi {)}^2+\varphi .$$

$$Q_c\equiv \varphi w{c}_l=\frac{Qc}{\mathcal{M}{D}_c}.$$

$${\varphi }_{\max }\sim {\left(\frac{F_{\max }}{\mathcal{Q}}\right)}^{1/n}={\left(\frac{F_{\max }\mu W_0}{k\Delta \rho g}\right)}^{1/n},$$

$$w_0\sim W_0\frac{F_{\max }}{{\varphi }_{\max }}.$$

$$t_{\mathrm{dry}}\sim n\frac{H_{\mathrm{dry}}}{w_0},$$

$$\delta ={\delta }_0{\mathrm{e}}^{i\omega t},$$

$$S=S_0+\frac{\Delta S}{H}{\mathrm{e}}^{i\omega t},$$

$$\varphi =\overline{\varphi }(z)+{\varphi }^{\prime }(z,t),c=\overline{c}(z)+c^{\prime }(z,t),$$

$$Q_c^{\prime }=\frac{Q^{\prime }\overline{c}+\overline{Q}{c}^{\prime }}{\mathcal{M}{D}_c}.$$

$$t_p^{\ast }=C\frac{1}{n(n-1)}{t}_{\mathrm{dry}}=C\frac{1}{n-1}\frac{H_{\mathrm{dry}}}{w_0},$$

$$A_Q\sim \frac{\Delta S}{W_0}\frac{{\rho }_w}{\rho }\frac{2\pi }{t_p},$$

Full text

mid-ocean ridges, glacial cycles, magmatism, carbon flux

1 Introduction

Terrestrial climate has changed dramatically between glacial and interglacial periods during the Pleistocene epoch. Continental ice sheets grow during glacial periods, causing a drop of up to 130 m in eustatic sea level; this decrease is recovered during interglacials (Lisiecki and Raymo, 2005; Bintanja et al., 2005; Waelbroeck et al., 2002; Siddall et al., 2010). The shift of mass-loading between continents and oceans affects subaerial and submarine volcanism (Jull and McKenzie, 1996) and its consequent carbon transfer from the solid Earth into the atmosphere/ocean (Huybers and Langmuir, 2009). During the Pleistocene, the climate system has varied on time-scales associated with Milankovitch orbital periods (Hays et al., 1976), indicating that glacial cycles are externally forced by variations in insolation (and its distribution over latitude). Huybers and Langmuir (2009) argued that glacial–volcanic coupling creates an internal amplification of climate variations. They later hypothesised (Huybers and Langmuir, 2017) that this feedback explains why glacial cycles of the last $\sim$700 kyr have a period of 100 kyr whereas the dominant Milankovitch forcing has a periodicity $\sim$41 kyr (Abe-Ouchi et al., 2013). However, the climate-feedback hypothesis depends on the amplitude of the glacial/interglacial fluctuation in carbon emissions and its lag with respect to changes in sea level. The present manuscript aims to develop a rigorous theory to predict the amplitude and lag.

Eighty percent of global volcanism occurs at mid-ocean ridges (MOR), where tectonic-plate divergence induces upwelling of the underlying mantle. Here, magma is produced by decompression melting, a quasi-isentropic process (McKenzie and Bickle, 1988; Langmuir et al., 1992). Up to about 20% partial melting occurs in the upper $\sim$100 km of the mantle, within a zone that extends $\sim$100 km on each side of the plate boundary (Katz, 2008; Keller et al., 2017). The melt segregates under its buoyancy, which supplies magma to the ridge axis and forms the oceanic crust. Because melting is driven by pressure change, and because variations of sea level affect the static pressure below, it is plausible that glacial cycles modulate magmatic production (Huybers and Langmuir, 2009; Lund and Asimow, 2011). Indeed a simple estimate and more detailed calculations by Crowley et al. (2015) indicate that crustal thickness could change by $\sim$10% (see also Tolstoy, 2015). Crowley et al. (2015) argued that such fluctuations provide a mechanism for the formation of abyssal hills with Milankovitch periodicity. This idea is controversial; the standard theory holds that seafloor topography is tectonically controlled (e.g., Olive et al., 2015). Lund and Asimow (2011) hypothesized that sea-level variations impacted the hydrothermal activity and the geochemistry of seawater, which has found some support in sedimentological records (Lund et al., 2016; Costa et al., 2017).

The mantle contains about 104 times more carbon than the atmosphere and ocean combined (Sleep and Zahnle, 2001; Dasgupta and Hirschmann, 2010; Hirschmann, 2018). This carbon is transferred from the mantle into the ocean-atmosphere system by volcanism and returns to the mantle in subduction zones (Zahnle and Sleep, 2002; Kelemen and Manning, 2015). At mid-ocean ridges, Hirschmann (2018) estimated that 120$\pm$26 Mt CO2/yr is extracted and emitted from the solid Earth; the majority of studies cited by Hirschmann (2018) have estimates that differ by within a factor of $\sim$3. The carbon flux could be modulated by variations in MOR magmatism during glacial cycles. Burley and Katz (2015) hypothesised that the key coupling mechanism is the pressure-driven variation in the depth of the onset of silicate melting. A drop in sea level reduces the static pressure and hence deepens the onset of first silicate melting. The downward motion of this boundary enhances the flux of carbon into the melting region. The opposite occurs when sea level rises. Burley and Katz (2015) predicted that 100 m changes in sea level with periodicity in the range 20–100 kyr would produce fluctuations up to 10% in the emission rate of CO2. Their hypothesised mechanism creates a lag of 50–80 kyr between a peak in the forcing and a peak CO2 emission rate. This delay arises from the time required for carbon-enriched (or depleted) melts to travel from the base of the silicate melting region to the surface. Huybers and Langmuir (2017) found that lags of 10–50 kyr are conducive to a negative feedback that would pace glacial cycles at a frequency of 1/100 kyr*-1* during the Pleistocene epoch.

Burley and Katz (2015) invoked significant assumptions and approximations in the development of their model. Two assumptions are particularly relevant here. First, they neglected the fluctuations of melting rate, porosity, and melt-transport speed that arise from sea-level variation (Lund and Asimow, 2011; Crowley et al., 2015). Second, they neglected the effect of carbon on mantle melting, which is to drastically lower the solidus temperature at constant pressure or, equivalently, to increase the pressure of first melting at constant entropy (e.g., Gaetani and Grove, 1998; Dasgupta and Hirschmann, 2006). Indeed, many studies have shown that even low concentrations of volatiles ($\sim$100 ppm) can induce the formation of low-degree melts at depths far greater than that of the anhydrous solidus temperature (e.g., Dasgupta et al., 2013; Dasgupta, 2018).

The present study aims to develop a more robust mathematical theory by removing simplifying assumptions made by previous work (Crowley et al., 2015; Burley and Katz, 2015). Our models are more complex than previous work in terms of the thermodynamic effect of CO2 and consistently modelling all the consequences of pressure fluctuations. In particular, our theory accounts for both the pressure-induced variations in the onset of (volatile-enriched) melting and also the pressure-induced variations in the melting rate. In this context, we show that the melting-rate variations created by oscillating sea level cause melt-flux variations, and that these are the primary cause of variations in CO2 emissions. Variations in the carbon concentration are secondary and variations in the onset depth of melting are inconsequential. As in previous work, we quantify the model sensitivity to sea-level change in terms of the admittance, which is defined as the amplitude ratio of response to forcing, as a function of frequency. We obtain a similar admittance of carbon flux as did Burley and Katz (2015), but we find that carbon emissions lag the causative changes in sea level by less than 20 kyr, assuming that melt ascent at mid-ocean ridges is no slower than 1 m/yr. Faster melt extraction corresponds to shorter lags. The lag is much shorter than the melt travel time because extra melts are generated throughout the column.

The manuscript is structured as follows. In Section 2, we describe the physical model, the governing equations and the mathematical strategy used to analyse them. In Section 3, we describe the results. We present the steady and time-dependent model predictions of porosity, carbon concentration, melt flux and carbon flux. We determine the admittance and lag of the fluxes as a function of the period of the sea-level cycle. In Section 4, we compare our results to previous models (Crowley et al., 2015; Burley and Katz, 2015) so as to isolate and discuss the differences. Finally, we consider a calculation of melt and carbon flux variations arising from a reconstruction of sea level over the past 800 kyr. Appendices provide details of the derivations and analyses that support our findings.

2 Methods

Here, a physical model is developed in mathematical terms to quantify the effect of sea-level variations on the melt and carbon flux to the ocean–atmosphere system. We generalize a standard, steady-state, one-dimensional melting-column model of the upwelling mantle (Ribe, 1985). The crucial generalization is to account for the time-dependent melting caused by a time-dependent sea level. We also calculate the evolution of the carbon concentration and account for its thermodynamic effect (Dasgupta and Hirschmann, 2006). Figure 1 is a schematic diagram of the melting column, depicting the effect of variable sea level and carbon on the depth at which the mantle crosses the solidus temperature. A pseudo-2d model is constructed from a series of columns.

2.1 Physical model

The physical model consists of a mechanical model of two-phase flow, a thermo-petrological model of melting, and a chemical model of carbon transport.

2.1.1 Mechanical model

Deep beneath a mid-ocean ridge, the mantle upwells at a rate $W_{0}$. At some depth (pressure) during this rise, its temperature reaches and exceeds the solidus temperature; at this point, it begins partial melting. The melt produced has a lower density than the residual solid and so rises buoyantly. This difference is small relative to the mean density, and hence we make a Boussinesq approximation and neglect density differences between phases, except when calculating the buoyancy. The volume fraction (or equivalently mass fraction) occupied by the melt is $\phi$. To quantify the key physical controls on the system, we derive a set of non-dimensional governing equations. In these equations, we non-dimensionalize lengths with the height of the melting column $H$, velocities with $W_{0}$, and time with $H/W_{0}$. However, we keep thermodynamic variables (temperature $T$ and pressure $P$) as dimensional quantities.

Bulk mass conservation (averaged across the two phases) requires that the bulk, one-dimensional flux is constant and equal to the mantle upwelling rate. In non-dimensional units, bulk mass conservation is

[TABLE]

where $w$ is the liquid velocity and $W$ is the solid velocity. This expression allows us to determine the liquid flux $Q\equiv\phi w$ in terms of the solid flux.

The volumetric melting rate, hereafter referred to as the melting rate, can be determined by considering mass conservation in the solid phase and is given by

[TABLE]

Melt segregates from the solid under gravity $g$ because it is less dense by an amount $\Delta\rho$. The melt has a viscosity $\mu$ and, at the scale of mantle grains, it inhabits an interconnected network of pores with permeability $k\phi^{n}$ (e.g., Miller et al., 2014; Rudge, 2018). The prefactor $k$ is a constant. The exponent $n$ determines the sensitivity of permeability to porosity and has been estimated to be $2\leq n\leq 3$. Then the melt segregation velocity is

[TABLE]

where the implication follows from equation (1). The parameter

[TABLE]

is the ratio of the rate of buoyancy-driven magma segregation to the rate of mantle upwelling. Thus, given the porosity, the solid flux can be determined from equation (3) and then the liquid flux from equation (1).

In this formulation, we have neglected the isotropic stress associated with compaction. Compaction stress typically varies on length scales much shorter than $H$ and plays a role only in narrow boundary layers (Ribe, 1985). However, compaction stresses can give rise to transient features called compaction waves or magmons (Scott and Stevenson, 1984; Richter and McKenzie, 1984). These could potentially interact with the time-dependence caused by sea-level variation, modifying the rate of chemical transport (Jordan et al., 2018).

2.1.2 Thermo-petrological model

The thermo-petrological model is used in concert with energy conservation to determine the melting rate. A simple parameterisation (Rees Jones et al., 2018; Bo et al., 2018) of solidus temperature that increases with pressure and decreases with carbon content is given by the linear relationship

[TABLE]

where $T$ is the solidus temperature, $P$ is total pressure, $\gamma$ is the slope of the carbon-free solidus, $M$ is the dependence of the solidus temperature on the carbon concentration $c_{0}$ of the unmelted mantle, and $c_{s}$ is the carbon concentration of the solid phase after it has been scaled by $c_{0}$. The last definition means that the $c_{s}=1$ before the onset of melting. $T_{0}$ and $P_{0}$ are the temperature and pressure at the onset of melting.

We introduce a dimensionless, time-dependent sea level $S(t)$ with time mean $S_{0}$; we introduce an independent vertical coordinate $z$. The origin $z=0$ is taken to be the bottom of the melting column, the depth at which upwelling mantle achieves the solidus temperature in the absence of sea-level variation ($S=S_{0}$). All lengths are scaled by $H$.

The total pressure is affected by sea level. Neglecting dynamic (compaction) pressure, the total pressure within the melting column is

[TABLE]

Deviations of sea level from $S_{0}$ lead to variations in the depth of the onset of mantle melting at

[TABLE]

Energy conservation is given by a temperature equation that accounts for advective transport and latent heat release (Ribe, 1985)

[TABLE]

where $L$ is the latent heat and $c_{p}$ is the specific heat capacity. Note that $\mathrm{D}/\mathrm{D}t\equiv\partial/\partial t+\partial/\partial z$ is a Lagrangian derivative accounting for bulk advection, i.e., advection by both the solid and liquid phases. We do not consider the effect of compressibility, which is relatively small. We substitute equations (5) and (6) into (8) to obtain the melting rate

[TABLE]

where the dimensionless parameters are defined by

[TABLE]

$\Gamma^{*}$ is the isentropic, dry melt productivity and $\mathcal{M}$ is the significance of carbon for melting. The latter is the ratio of the sensible heat associated with the solidus-depression caused by carbon to the latent heat. Dry melting corresponds to $\mathcal{M}=0$; carbonated melting to $\mathcal{M}>0$. Equation (9) shows that there are contributions to the melting rate from standard decompression melting $\Gamma_{0}$, from the pressure change associated with sea-level variation $\Gamma_{P}$, and from the effect of carbon on melting $\Gamma_{c}$. Sea-level variation has a direct effect on melting since

[TABLE]

Sea-level variation also has an indirect effect on the melting rate by changing the concentration of carbon, a mechanism captured in the term $\Gamma_{c}$.

The behaviour of the melting rate has a crucial transition between low-degree, carbonated melting and dry melting. This occurs when $z=\mathcal{M}/\Gamma^{*}$. Indeed, we can make an alternative interpretation of $\mathcal{M}$ in terms of the total depth $H$ of the melting region (the carbonated solidus depth) relative to the dry solidus depth $H_{\mathrm{dry}}$. In particular

[TABLE]

Note that $H_{\mathrm{dry}}={\gamma\Delta T}/{\rho g}$, where $\Delta T$ is the difference between the mantle potential temperature and the solidus temperature at the top of the melting column (Ribe, 1985). The maximum degree of melting $F_{\mathrm{max}}$ is equal to the liquid flux at the top of the column (Ribe, 1985). This satisfies, to an excellent approximation at small porosity,

[TABLE]

Thus given estimates for $H_{\mathrm{dry}}$, $H$ and $F_{\mathrm{max}}$, the dimensionless parameters $\Gamma^{*}$ and $\mathcal{M}$ can be uniquely determined (see Appendix A for details).

2.1.3 Chemical model

The final part of the model captures the transport of carbon. Diffusion and dispersion of carbon are slow compared to advection, so conservation of carbon in the two phases is given by

[TABLE]

where $c_{l}$ is the dimensionless concentration of carbon in the melt (again scaled by $c_{0}$).

We assume that carbon behaves as an incompatible element (e.g. Rosenthal et al., 2015) that partitions preferentially into the melt according to

[TABLE]

where $D_{c}$ is a partition coefficient for carbon. Equation (15) allows us to eliminate $c_{l}$ in equation (14) and obtain a single equation for carbon transport.

2.1.4 Boundary conditions

At the depth of the onset of melting $z=\delta(t)$, there is no porosity and the carbon concentration is that of the far-field, upwelling mantle. Thus appropriate boundary conditions are

[TABLE]

2.1.5 Governing equations

To synthesise the model components above, we combine the mechanical, thermal and chemical equations to obtain a coupled system. We eliminate the liquid and solid fluxes using equations (1) and (3) and the melt rate using equations (2) and (9). Then the system of equations for the evolution of porosity and carbon concentration is

[TABLE]

where we simplified the notation by defining a scaled carbon concentration and effective partition coefficient, respectively,

[TABLE]

The liquid flux $Q$ depends on porosity according to

[TABLE]

We also evaluate the carbon flux in the melt, $Q_{c}$, which is defined by

[TABLE]

2.2 Decomposition into steady and fluctuating components

The melting column has a steady-state, mean behaviour in the absence of sea-level fluctuations. Figure 2 shows a typical steady-state solution of the governing equations. Table 1 gives the standard set of parameters we use in calculations unless otherwise stated. The main sensitivity is to the speed of melt flow and the mantle upwelling rate (which depends on the spreading rate). Both of these effects are contained within the parameter $\mathcal{Q}$. Our standard choice of $\mathcal{Q}$ results in melt flow that is comparable in magnitude to previous studies (Crowley et al., 2015; Burley and Katz, 2015) but slower than has been inferred for melt velocity based on the Iceland post-glacial melt pulse (Maclennan et al., 2002; Swindles et al., 2017; Eksinchol et al., 2019)). We explore model sensitivity to the choice of $\mathcal{Q}$ in figure 5 and the Supplementary Information (section S1).

The melting column (Figure 2) can be divided into three regions: ‘wet’, ‘transitional’ and ‘dry’, as labelled on the figure. In the wet region, occupying roughly the bottom half of the column, a small amount of carbon-rich melts are generated. The porosity and hence the liquid flux remain very small. There is almost no melt segregation. In the dry region, occupying roughly the top half of the column, the degree of melting increases. The liquid flux increases and the solid flux decreases. The transport of carbon becomes dominated by the liquid phase. Around the depth of the carbon-free solidus, there is a transitional region separating the wet and dry regions. We discuss the steady-state behaviour and choice of parameters further in A.

Two important scales emerge from the steady solution: the maximum porosity at the top of the column

[TABLE]

and the liquid velocity at the top of the column, which, in dimensional units, is given by

[TABLE]

Thus, if the melt flux parameter $\mathcal{Q}$ is large, melt extraction is efficient and so the maximum porosity is small and the melt velocity is large. The dimensional melt transport time $t_{\mathrm{dry}}$ across the dry-melting region is approximately given by

[TABLE]

where the factor of $n$ arises from the $z$-dependence of the velocity.

The variations in sea level associated with glacial cycles are on the scale of 100 m, while the depth of the melting region is tens of kilometres. Furthermore, the density ratio $\rho_{w}/\rho\approx 0.3$. These considerations imply that the variations $\delta$ in the depth of first partial melts associated with sea-level change are small, according to equation (7). Therefore, we assume that the time-dependent fluctuations of all quantities are relatively small and linearize the governing equations about the steady state. We further decompose the time-dependent fluctuations into harmonics of dimensionless frequency $\omega$. Thus

[TABLE]

where $\delta_{0}\ll 1$ is the maximum fluctuation. This is appropriate for either periodic forcing of frequency $\omega$ or by taking the Fourier transform of a record of sea level (in the latter case, $\delta_{0}$ is a function of $\omega$). The dimensionless sea level satisfies

[TABLE]

where $\Delta S$ is the maximum dimensional sea-level fluctuation. If we define $\Delta\delta\equiv 2\delta_{0}$, then $\Delta S=H\Delta\delta(\rho/\rho_{w})$.

We decompose the porosity and bulk carbon concentration into steady parts (denoted with an overline) and time-dependent parts (denoted with a prime),

[TABLE]

where $\phi^{\prime}(z,t)=\delta_{0}\mathrm{e}^{i\omega t}\hat{\phi}(z)$ and $c^{\prime}(z,t)=\delta_{0}\mathrm{e}^{i\omega t}\hat{c}(z)$. We derive equations governing the fluctuations of porosity and carbon concentration in Appendix B. These are used to deduce the contribution of sea-level variations to the melt and carbon fluxes. Equivalent notation is used for the decomposition of the fluxes into steady and fluctuating parts. The fluctuating part of the carbon flux has contributions from both the fluctuation of melt flux and carbon concentration:

[TABLE]

Software to reproduce the calculations and the results shown in the figures in the remainder of this paper is available (Cerpa et al., 2019).

2.3 Pseudo-two-dimensional model of melt focusing

The melting region beneath a mid-ocean ridge is not columnar; rather it is a volume that encloses upwelling, melting mantle (Forsyth et al., 1998). In a vertical plane normal to the ridge axis, the shape of the melting region can be approximated as triangular (Langmuir et al., 1992). Magma produced off-axis is focussed along a decompaction channel at the base of the lithosphere toward the ridge axis (Sparks and Parmentier, 1991).

Figure 1 illustrates our pseudo-two-dimensional model (see Appendix D for full details). We assume that the melting region comprises an array of independent columns that deliver magma into a decompaction channel (Sparks and Parmentier, 1991). The decompaction channel transports both the mean flux and variations. Following previous work, we consider two simple assumptions for this transport: instantaneous (Burley and Katz, 2015) and finite-rate (Crowley et al., 2015). Magma and carbon flux variations are delivered to the ridge according to an integral over columns from the axis out to some maximum focusing distance (equation (56)). In the case of finite-rate focusing, the transport time causes a phase-delay $\tau$ that increases with distance to the ridge $x$.

3 Results

Variation in sea level causes variation in pressure and hence variation in (i) the onset depth of melting and (ii) the melting rate throughout the column. We refer to (i) as the ‘basal-flux mechanism’ and to (ii) as the ‘internal-melting mechanism’. These mechanisms, in turn, drive variation in the melt and carbon flux.

3.1 Example of fluctuations due to sea-level changes

Figure 3 shows the response of porosity, carbon concentration, melt flux and carbon flux to a sea-level cycle that has a peak-to-trough magnitude of 100 m and a period of 100 kyr. While these numbers are chosen for illustration, they roughly correspond to the sea-level variation experienced in the late Pleistocene (past 800 kyr).

We first discuss the coupled evolution of the porosity and carbon concentration, since these are the primary fields. We discuss the evolution working from the bottom of the melting column to the top, following the direction of the steady-state liquid and solid flow.

Figure 3a shows the porosity fluctuation. In the wet-melting region, near the bottom of the melting column, the basal-flux mechanism is significant. Here the steady-state porosity increases with height (figure 2a) and hence, in order to maintain zero porosity at the onset of melting $z=\delta(t)$, a positive sea-level fluctuation leads to a negative porosity fluctuation. Conversely, the steady-state carbon concentration decreases with height (figure 2b). So a positive sea-level fluctuation leads to a positive carbon fluctuation (figure 3b). Hence the porosity is in antiphase with sea level but the carbon concentration is in phase. The magnitude of the porosity variation is relatively small because the steady-state porosity is also very small. By contrast, the magnitude of the carbon variation is more significant because of the sharper variation in the steady carbon concentration. The internal-melting mechanism causes little change to $\phi$ and $c_{l}$ because there is no melt segregation in the wet-melting region. Thus there is an almost perfect balance between the melting-rate fluctuation caused by decompression melting and that caused by the carbon concentration fluctuation. Carbon is buffering the system to counteract the extra internal decompression melting caused by sea-level variation. In terms of the melting rate given by equation (9), $\Gamma_{P}\sim-\Gamma_{c}$.

In the transitional region, around the depth of the dry solidus, the fluctuations of porosity and carbon concentration change significantly. It is in this zone that the melt segregation becomes important, since this is where the steady porosity and melt flux start to increase significantly. Melt segregation breaks the buffering capacity of carbon described previously $(\Gamma_{p}\nsim-\Gamma_{c})$ because $\Gamma_{c}$ is reduced in magnitude. Indeed, segregation leads to a decrease in the magnitude of the carbon fluctuation. There is also a shift to a slightly positive phase (N.B., throughout the manuscript, a ‘positive phase’ is achieved when the peak of a fluctuation shifts to earlier times, as in the transition region in figure 3b). Furthermore, the porosity fluctuation increases with depth and undergoes a phase shift, which we explore below.

In the dry-melting region, near the top of the column, the dominant contributions to porosity evolution come from internal decompression melting $\Gamma_{p}$ and fluctuations in the upward transport of melt. The resultant phase observed is intermediate between the phase associated with the basal flux ($-S$) and that associated with internal melting ($-\dot{S}$, where a dot represents a time derivative). Concurrently, the carbon fluctuation continues to decrease in magnitude, in the manner outlined above.

The melt flux fluctuation (figure 3c) has the same phase as the porosity fluctuation. The amplitude of the flux increases more rapidly with height than porosity because the steady-state flux also increases. The carbon flux fluctuation (figure 3d) has contributions from both the melt flux and carbon concentration fluctuations, which must be weighted by the steady state as given by equation (27). The former contribution is larger, so the carbon flux largely follows the melt flux fluctuation. The carbon concentration fluctuation has the opposite phase so, in part, offsets the melt flux fluctuation and slightly advances the phase, so that the peak in carbon flux occurs slightly earlier than the peak in melt flux.

3.2 Effect of the period of sea-level fluctuations

These flux variations are sensitive to the period of forcing. Figure 4 shows the melt and carbon flux fluctuations at three different periods, reflecting the dominant periods of sea-level variations in the late Pleistocene. We focus our discussion on the behaviour at the top of the melting column. At very short periods (e.g., 23 kyr), the melt and carbon flux fluctuations are in antiphase with sea level. With increasing forcing period, there is a switch in behaviour. For a sufficiently long period (e.g., 100 kyr), the system behaves as described in section 3.1, with fluxes proportional to the rate of decrease in sea level. This results in a phase advance of about a quarter of a period as the forcing period increases.

The behaviour of the system with a short forcing period (high frequency) can be approximated using an analytical method. A detailed calculation is presented in Appendix C. Here, we describe the main physical ideas and insights of this analysis. At high forcing frequency, melt segregation during one cycle is minimal and the porosity and carbon concentration respond almost instantaneously throughout the melting column to the pressure change induced by the change in sea level.

This behaviour changes when the forcing period is proportional to the melt travel time across the dry melting region $t_{\mathrm{dry}}$ given by equation (23), and occurs at a dimensional critical period $t_{p}^{*}$ given by

[TABLE]

where $C\approx 1$ is a dimensionless prefactor (roughly corresponding to a 20% phase shift). For the parameters given in table 1, the critical period is about 23 kyr, consistent with figure 5b described below. To a good approximation, the critical period does not depend on the solidus-depressing thermodynamic effect of carbon. Similar arguments can be applied to the evolution of carbon, in which case dilution by melt transport determines the critical period given by equation (52b) in Appendix C.

3.3 Admittance and lag

Figure 5 summarizes the effect of the period of sea-level forcing in terms of the admittance and lag of melt and carbon fluxes. Here the admittance (panel a) is defined as the peak-to-trough magnitude of the time-dependent part of the flux at the top of the melting column normalized by the steady-state flux. This quantity (sometimes called ‘relative admittance’) is proportional to the amplitude of the sea-level cycle. Therefore, we report admittance as a percentage per 100 m peak-to-trough sea-level fluctuation. We define the lag (panel b) as the difference between the time of the peak flux and the time of peak rate of decrease in sea level $-\dot{S}$. This choice of baseline, while somewhat arbitrary, is motivated by the fact that it corresponds to a maximum in $\Gamma_{P}$, i.e., to the peak rate of generation of extra partial melts due to sea level.

The calculated lag at small forcing periods is approximately proportional to $t_{p}/4$ because the fluxes are antiphase with sea level, rather than in phase with the rate of decrease in sea level. At periods longer than the critical period of equation (28), the phase approaches that corresponding to the rate of decrease in sea level. Thus, although the calculated lag of melt flux continues to increase with increasing period (reflecting the finite melt transport time), it does so at a much smaller rate.

At long periods, the admittance decreases with period because longer period corresponds to slower sea-level change. The melt flux admittance $A_{Q}$ can be estimated (Lund and Asimow, 2011) by comparing the effect of sea-level variation on pressure to the static pressure and hence the relative melting rate ($\Gamma_{P}/\Gamma_{0}$). This gives

[TABLE]

which is a good estimate at periods longer than about 100 kyr [see Supplementary Information, section S2].

Conversely, at short periods, the admittance tends toward a constant that is independent of period. The magnitude of this constant can be estimated analytically, as done in Appendix C where we derive the following approximations for admittance of porosity $A_{\phi}$, melt flux $A_{Q}$ and carbon flux $A_{Q_{c}}$,

[TABLE]

where $\Delta S$ is the dimensional magnitude of the sea-level fluctuation and $w_{0}$ is the dimensional melt velocity at the top of the column. The physical meaning of these expressions can be interpreted as follows, making the approximation of negligible melt segregation on the timescale of one period of sea-level variation. The fluctuating part of the porosity is equal to the pressure change associated with sea level multiplied by the productivity. The steady-state melt flux (which is equal to the porosity multiplied by the melt velocity) is equal to the pressure change across the dry melting region multiplied by the productivity and the mantle upwelling rate. This allows us to estimate the steady-state part of the porosity. The combination of the estimates for the steady and fluctuating parts of the porosity shows that the admittance is equal to the ratio of melt to mantle velocity multiplied by the ratio of pressure change associated with sea level to that across the dry melting region, giving equation (30a). The admittance of melt flux is a factor of $n$ larger than that of porosity because melt flux increases with porosity as a power law with exponent $n$, so $A_{Q}=nA_{\phi}$, giving equation (30b). The fluctuation of carbon concentration is the opposite to that of porosity; this is required to keep the bulk concentration constant, so $c^{\prime}\overline{\phi}=-\overline{c}\phi^{\prime}$. Physically, an increase in porosity dilutes the carbon concentration in the melt. Finally, by combining this with the carbon flux fluctuation given by equation (27), we obtain equation (30c).

Figure 5(c,d) shows that this theory can be applied to obtain simple estimates of the maximum admittance and forcing period at which that maximum occurs. The maximum admittance of melt flux occurs at a forcing period $t_{p}^{\mathrm{max}}$ that is proportional to the melt transport time across the dry-melting region $t_{\mathrm{dry}}$, as previously suggested by Crowley et al. (2015). In particular, we write

[TABLE]

where $C_{Q}=0.67$ is a prefactor chosen so that the analytical estimations fit our model calculations at the reference value of $\mathcal{Q}$. This period is intermediate between the long- and short-period limits. The magnitude of the maximum admittance of melt flux is a factor of about $B_{Q}=1.12$, i.e., 12% greater than the approximate formula (30b) applicable for small forcing period. Physically, this maximum occurs when the porosity fluctuations caused by melting are positively reinforced by fluctuations of the upward transport of melt (these contributions reinforce each other at intermediate $t_{p}$).

The same approach can be applied to carbon fluxes. We use the same notation for the prefactors, except replacing subscript Q (melt flux) with subscript ${}_{Q_{c}}$ (carbon flux) and using equation (30c), which is the estimate of carbon flux admittance at small forcing period.

These results are sensitive to the fluid dynamical properties of the system. Figure 5(a,b) shows the effect of reducing melt flux parameter $\mathcal{Q}$ by a factor of 4. This could correspond, for example, to a reduction of permeability, an increase in melt viscosity or an increase in mantle upwelling rate. This reduces liquid velocities and hence increases the melt transport time. The admittance is reduced. At small periods, $A_{Q}$ is reduced by a factor of 2 ($=4^{1/n}$, since $n=2$, see equation (48)). At large periods, the effect is much more modest. Indeed, reducing the permeability can slightly increase the carbon flux fluctuation at sufficiently long periods. The lag is typically (but not always) increased by increasing the melt transport time (smaller $\mathcal{Q}$) and the critical period $t_{p}^{*}$ increases, consistent with equation (28). The behaviour of the carbon flux is complicated because it is affected by two contributions: one from the melt flux fluctuation and the other from the carbon concentration fluctuation, as shown in equation (50). The latter tends to be in phase with $S$ and so peaks earlier in the cycle, which leads to the negative lags with respect to $-\dot{S}$ at high forcing periods (figure 5b). Carbon concentration fluctuations are also responsible for the non-monotonic sensitivity of admittance to permeability.

4 Discussion

This study builds on previous work by incorporating fluctuations in the melting rate throughout the melting column (‘internal melting’) and also by considering the thermodynamic effect of carbon on melting. Burley and Katz (2015) considered only fluctuations introduced by variation in the melting-onset depth with sea level (‘basal flux’). Crowley et al. (2015) considered only internal melting and calculated only melt fluxes. They did not calculate the concentration of carbon, nor consider its thermodynamic effect.

The simplifying assumptions made in previous studies can be tested within our framework. First, the thermodynamic effect of carbon can be assessed using a ‘dry-melting model’ ($\mathcal{M}=0$). Second, the importance of internal melting can be assessed by excluding it from the equations governing the fluctuations (in addition to $\mathcal{M}=0$, we force $\Gamma_{P}=0$; see Appendix B). We label this a ‘basal-flux model’ because it includes only this mechanism.

Figure 6 shows that there are minimal differences between our dry- and wet-melting models near the surface. This is because carbon has only a minor effect on melting outside the wet-melting region. Fluctuations in melt flux and the melt-flux admittance (blue line figure 6d) behave similarly in both cases. In particular, $A_{Q}$ approaches a constant admittance at small forcing periods. This finding contrasts with Crowley et al. (2015), so we test whether this difference arises from our 1d simplification by creating pseudo-2d models using the same methodology (see section 2.3 and Appendix D). The admittance of the pseudo-2d model with instantaneous focusing (light blue, figure 6d) is similar to that of the 1d model. In contrast, when assuming a finite-focusing rate (dark blue, figure 6d, we predict sharp decrease in the melt-flux admittance at small period. The latter case is consistent with the model of Crowley et al. (2015). However, at forcing periods of Milankovitch cycles, there are minimal differences between 2d and 1d models. In addition to 2d effects, storage in crustal magma chambers and turbulent mixing of carbon from the MOR to the surface ocean and atmosphere would attenuate admittance at timescales less than about 1 kyr.

Figure 7 shows our basal-flux model, which we now compare to figures 6a and 7c in Burley and Katz (2015). In both cases, the admittance decreases with forcing period and the lag is insensitive to forcing period. Using the same parameters as Burley and Katz (2015) (thin line), we find a good agreement between our models. Thus we regard our basal-flux model as a one-dimensional representation of Burley and Katz (2015). In our model with reference parameters (thick line; larger permeability than in Burley and Katz (2015)), the admittance is increased and the lag is diminished. At a forcing period of 100 kyr, the admittance is two times higher and the lag is one third lower than in the models with values of Burley and Katz (2015). Next, we compare the basal-flux model, and thus Burley and Katz (2015), with the full wet-melting model (figures 4 and 5). The basal-flux model gives a slightly greater admittance of carbon flux at the dominant forcing periods and a different admittance structure at small periods. There are also profound differences in the timing of the surface fluxes because, in the basal-flux model, fluctuations created at the base must travel to the surface, which occurs approximately over the melt travel time. By contrast, accounting for internal melting means that much of the effect of sea-level fluctuation is generated closer to the surface (this is true independent of volatile content). Thus our more general models that incorporate both basal-flux and internal-melting mechanisms differ markedly from Burley and Katz (2015) in terms of the lag predicted.

Finally, we assess the implications of our calculations for carbon fluxes from mid-ocean ridges over the Pleistocene. Figure 8 shows our predictions based on the reconstructed sea-level variation over the past 800 kyr. Our calculations indicate that sea-level fluctuations have driven substantial variation in melt and carbon fluxes from the mid-ocean ridge system. The melt and carbon fluxes depart from the mean, steady-state values with a total range of about 20% and 13% respectively, and are therefore potentially significant contributors to variation in crustal thickness and variation in magmatic carbon fluxes to the ocean/atmosphere. Figure 9a,c shows that these estimates are particularly sensitive to the melt transport (reported as a function of steady-state melt velocity at the top of column) but not strongly influenced by 2d effects.

The melt and carbon fluxes in figure 8 are most closely related to the reconstructed rate of decrease in sea level $-\dot{S}$. Variations of this rate directly force variations in the rate of decompression melting $\Gamma_{P}$, emphasising the significance of the internal-melting mechanism. The negative flux excursions are more pronounced than the positive ones. This reflects the fact that glacial cycles are marked by gradual decreases in sea level (slow glaciation) and sharp increases (rapid collapse of ice sheets). The carbon flux variation (fig. 8d) is dominantly caused by the variation in melt flux (fig. 8c), rather than by the variation in carbon concentration. However, the carbon concentration plays a mitigating role, such that the carbon flux variation is about half the melt flux variation. Thus estimates of variation in melt flux from an observed variation in crustal thickness could be used to estimate variation in mid-ocean ridge carbon emissions. Figure 8e,f shows that the peaks of melt and carbon flux lag the forcing by about 5 kyr. Figure 9b,d shows that this estimates is sensitive to the melt transport and also influenced by 2d effects, particularly if lateral melt focussing to the ridge axis is slow.

Carbon is a significant greenhouse gas and so any variation in magmatic carbon flux can potentially act as a feedback on glacial cycles. Our model, driven by Pleistocene fluctuations of sea level (figure 8), indicates that rapid deglaciations were followed by a significant decrease in the carbon flux from mid-ocean ridges, which in turn would have reduced atmospheric carbon, a potentially significant negative climate feedback. Because this feedback is nonlinear, the periods of forcing and response can be different. Indeed, Huybers and Langmuir (2017) argued that a $\sim$120 kyr glacial cycle could arise from the 41-kyr orbital forcing when the lag of magmatic carbon emissions is 10–50 kyr. For Burley and Katz (2015), variation in MOR emissions is driven by changes in the basal carbon flux, which is caused by variation in the depth of first melting. This produces lags of about 50–80 kyr, controlled by the melt migration rate across the whole melting column. According to Huybers and Langmuir (2017), such lags would promote longer glacial cycles than those of the Pleistocene epoch. In the present work, we account for internal melting throughout the column. This leads to a shorter lag of about 5 kyr with reference parameters (figure 8f), which is much less than the melt travel time. Figure 9d shows that the lag could be even shorter if the melt velocity is faster than our preferred value, or up to about 20 kyr if the maximum melt velocity is slower (1 m/yr) and/or lateral melt focusing is slow enough. Hence, at least for the upper range of plausible lags that we predict, the mechanism proposed by Huybers and Langmuir (2017) is viable. It would be interesting to revisit these feedbacks with our revised model of carbon fluxes from mid-ocean ridges.

Our prediction of lag for mid-ocean-ridge emissions mainly depends on the segregation rate. Assuming melt ascent by diffuse porous flow, microstructural measurement of the permeability of rocks (e.g., Miller et al., 2014) suggests speeds of the order of 1 m/yr. However, observations of U-series disequilibria are consistent with melt ascent speeds of several tens of metres per year (Rubin and Macdougall, 1988; Stracke et al., 2006). And the magmatic response in Iceland to the last deglaciation indicates rates of 50 m/yr or higher (Maclennan et al., 2002; Eksinchol et al., 2019). Our lag predictions also depend on the rate of mantle upwelling. Current full-spreading rates range from roughly 1 to 15 cm/yr globally (Bown and White, 1994), inducing mantle upwelling rates in the range 0.5–10 cm/yr. At our reference permeability scale, this range of mantle upwelling rates would produce lags within the range given in Fig. 9. Overall, our models show that the lag is always less than about 20% of the melt travel time.

The theory developed here and in previous work calls out for a test by comparison with observations. Unfortunately, time-series of lava compositions with appropriate duration and resolution are unavailable (although see Ferguson et al. (2017)). Moreover, it is unlikely that carbon dioxide would leave any observable, temporal signal after degassing. However, records of hydrothermal elemental fluxes in sediments provide a proxy for temporal variations in hydrothermal activity and ultimately in magmatic budgets over the last glacial cycles (Lund and Asimow, 2011). The hydrothermal proxies in sediments can be accurately dated and the timing of peak hydrothermal activity can be compared to predictions from our theory. Lund et al. (2016) reported time-series of Fe and Mn fluxes in the sediments of the Southern East Pacific Rise at 11∘S over the last 200 ka. They showed that the peaks in these fluxes lag the two previous maxima in the rate of sea-level decrease by about 15 kyr. Middleton et al. (2016) found that sediments from 26∘N on the Mid-Atlantic Ridge document an increase in elemental fluxes (Fe, Cu) in hydrothermal systems concomitant with the most rapid sea-level decrease leading to the last glacial maximum. Relatively short lags between peaks in rate of sea-level change and peaks in melt flux are consistent with our models forced with the sea-level reconstruction of the Middle and Upper Pleistocene. Furthermore, consistent with the control by parameter $\mathcal{Q}$ on the lag of melt flux (figure 5, see also figure 9), the higher lags observed at the Southern East Pacific Rise might be due to the greater half-spreading rate there, compared to that of the Mid-Atlantic Ridge.

We found differences between one- and two-dimensional models only at small forcing frequency. However, our pseudo-2d models are limited in that they do not consider potential complexities of flow focusing to the axis or lateral variations within the melting region. The latter could be caused by changes in mantle upwelling rate or melt-localisation instabilities (e.g., Kelemen et al., 1995; Keller et al., 2017; Rees Jones and Katz, 2018). To improve our understanding of the relationships between glacial cycles and mid-ocean ridge magmatism, next steps should include forcing two-dimensional models of MOR magmatism with local reconstructions of sea-level variations. They should also include comparisons with time-series of hydrothermal activity, trace elements in basaltic glass chips, and oceanic crustal thickness.

Acknowledgements

The authors thank P. Asimow and P. Huybers for insightful reviews and J.F. Rudge for comments on an early version. This research received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement number 279925 and under Horizon 2020 research and innovation programme grant agreement number 772255. N.G.C. acknowledges support from the University of Montpellier and public funding through ANR under the “Investissements d*′*avenir” programme with the reference ANR-16-IDEX-0006. D.R.J. acknowledges research funding through the NERC Consortium grant NE/M000427/1, NERC Standard grant NE/I026995/1 and the Leverhulme Trust. We thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme Melt in the Mantle that was supported by EPSRC Grant Number EP/K032208/1. We thank the Deep Carbon Observatory of the Sloan Foundation.

Appendix A Steady state for fixed sea level, parameter estimation

The steady-state porosity and carbon concentration are obtained by taking the steady state of equation (17), integrating once with respect to $z$, and applying boundary conditions ([math]). We find

[TABLE]

where $\overline{Q}=Q(\overline{\phi})=\mathcal{Q}\overline{\phi}^{n}(1-\overline{\phi})^{2}+\overline{\phi}$. These algebraic equations can be readily solved; an example solution is shown in figure 2.

The steady-state melting model can be analyzed to understand its basic behaviour. In turn, this analysis can be used to estimate the dimensionless parameters of the system in terms of (relatively) easy-to-measure quantities. The porosity in the melting column is controlled by the balance of melt production and melt extraction. As melt is generated, the porosity increases but the carbon concentration decreases because carbon partitions as an incompatible element into the melt. We can estimate the porosity $\phi_{\mathrm{max}}$ and carbon concentration $c_{\mathrm{min}}$ at the top of the column $z=1$ by assuming that the porosity is small and carbon is very incompatible. We define precise conditions necessary for these limits below. Under these approximations

[TABLE]

The corresponding concentration of the solid and liquid phases, respectively are

[TABLE]

We then determine appropriate conditions under which the approximations hold by estimating a posteriori the magnitude of each of the neglected terms. In particular, we require

[TABLE]

The first expression ensures that the term $+\overline{\phi}$ can be neglected in the flux $\overline{Q}$. Then $\overline{\phi}\ll\Gamma^{*}-\mathcal{M}<1$. The second expression ensures that $\overline{c}\ll\Gamma^{*}-\mathcal{M}$. Physically, these conditions mean that the melt extraction is fast enough to maintain a small porosity and the compatibility is low enough that carbon does not suppress the solidus much near the top of the column.

We now relate our governing parameters to other quantities of interest, which may be easier to measure or estimate in practice. First, in the limit of small porosity, the maximum flux is equivalent to the maximum degree of melting $F_{\mathrm{max}}$ (Ribe, 1985). In particular,

[TABLE]

So, using equation (33a), the maximum porosity satisfies

[TABLE]

Recall that $\mathcal{M}$ is related to $\Gamma^{*}$ through equation (12), so

[TABLE]

For example, we take the following broadly accepted values from the literature: $F_{\mathrm{max}}\approx 0.2$, the carbonated solidus is at $H\approx 130$ km, and the dry solidus is at $H_{\mathrm{dry}}\approx 65$ km. Then $\Gamma^{*}=0.4$ and $\mathcal{M}=0.2$ (e.g., Klein and Langmuir, 1987; Forsyth et al., 1998).

Second, we estimate the liquid velocity at the top of the column

[TABLE]

Given an estimate of the melt velocity (say from observations of uranium-series disequilibrium (Rubin and Macdougall, 1988; Stracke et al., 2006) or the Iceland post-glacial melt pulse (Maclennan et al., 2002; Swindles et al., 2017; Eksinchol et al., 2019)), we can estimate

[TABLE]

For example, if the melt velocity is 140 times faster than the mantle upwelling rate, $w_{0}/W_{0}=140$. If also $n=2$, which is appropriate given that the porosity is small (Rudge, 2018), then $\mathcal{Q}=10^{5}$. As expected, $\mathcal{Q}$ is relatively large since melt extraction is fast and porosities remain small. Also, $D\sim 10^{-4}$ for carbon (Rosenthal et al., 2015) so all the constraints required in equation ([math]) are satisfied.

Despite the consensus on the physical quantities above (similar parameters are used in previous studies (Crowley et al., 2015; Burley and Katz, 2015)) there is uncertainty arising both from the indirect nature of the constraints and from geographical variation. This is particularly true of $\mathcal{Q}$, which depends on permeability, spreading rate and melt viscosity. We explore model sensitivity to the choice of $\mathcal{Q}$ in figures 5 and 9 and the Supplementary Information (section S1).

Appendix B Effect of fluctations in sea level

The equations governing the time-dependent fluctuations are obtained by linearizing equation (17) about the steady part $(\overline{\phi},\,\overline{c})$ of the solution. In particular, we neglect terms that contain $\delta_{0}^{2}$ or higher powers. We collect terms proportional to $\delta_{0}$ and find

[TABLE]

Note that $\hat{Q}=\hat{\phi}\,\tfrac{d\overline{Q}}{d\overline{\phi}}$. Appropriate boundary conditions are found by linearizing equation ([math]) to obtain

[TABLE]

Hence equation (41) is a pair of coupled ordinary differential equations that can be solved for $\hat{Q}$ (and hence $\hat{\phi}$) and $\hat{c}$.

We also develop a ‘basal-flux’ model akin to that of Burley and Katz (2015). In this model, we neglect the contributions from internal melting. In particular, we neglect the term $-i\omega\Gamma^{*}$ in equation (41a) and $i\omega\overline{c}\Gamma^{*}$ in equation (41b). Results using this model are presented in figure 7 and discussed in section 4.

Appendix C Approximate solutions valid in the limit of large forcing frequency

When the period of sea-level fluctuation is short, the frequency $\omega$ is large. In this limit, we derive approximate solutions of the equations given in Appendix B. These approximate solutions are useful in that they allow us to identify the physical mechanisms of importance in this regime (as described in the main text), as well as simple estimates of quantities of interest, such as the melt and carbon flux.

The key methodological idea is that the imaginary part of the fluctuating quantities is very much smaller than the real part, by a factor of $\omega^{-1}\ll 1$. Mathematically, this can be seen by inspection of equation (41) (we discuss the physical meaning in section 3.2). When $\omega$ is very large, the collective group of terms multiplied by it must be very small. This then shows that $\hat{\phi}$ and $\hat{c}$ must be (approximately) equal to some real function of $\Gamma^{*}$, $D$, $\overline{\phi}$ and $\overline{c}$, all of which are real. So the fluctuations are approximately equal to their real parts. Then it can be seen that the imaginary parts are a factor $\omega^{-1}\ll 1$ smaller. With this in mind, we write:

[TABLE]

In this expansion, a subscript $r$ denotes the real part of a quantity and a subscript $i$ denotes the imaginary part scaled with the frequency. Thus, for example, both $\hat{Q}_{r}$ and $\hat{Q}_{i}$ are real quantities. We substitute this decomposition into equation (41), take the real and imaginary parts, and collect powers of $\omega$.

The leading order balances, in which we collect terms proportional to $\omega$, gives

[TABLE]

These expressions can be straightforwardly rearranged to give explicit solutions for the real part of all the fluctuations in terms of the mean state. In particular

[TABLE]

Note that $\hat{Q}_{r}=\hat{\phi}_{r}\,\tfrac{d\overline{Q}}{d\overline{\phi}}$.

Then, the next order balances, in which we collect terms proportional to $\omega^{0}$, gives

[TABLE]

As before, we can rearrange these expressions to give explicit expressions for $\hat{\phi}_{i}$ and $\hat{c}_{i}$, since the terms involving real parts of the fluctuations are known from equation (45). We can also obtain the imaginary part of the melt flux fluctuation by noting that $\hat{Q}_{i}=\hat{\phi}_{i}\,\tfrac{d\overline{Q}}{d\overline{\phi}}$.

Our approximate expressions for the real part of all the fluctuating quantities allow us to estimate the admittance (in terms of melt and carbon fluxes), since the full fluctuation is dominated by the real part. The admittance of melt flux (in units of per 100 m peak-to-trough sea-level fluctuation) is

[TABLE]

where all quantities are evaluated at the top of the melting column (and we convert to a percentage when plotting). In A, we discussed how the steady-state variables can be simplified at the top of the melting column. If we make the approximation that $D\ll\overline{\phi}\ll 1$ and $\overline{c}\ll\overline{\phi}$, then

[TABLE]

where $\Delta S$ is the dimensional magnitude of the sea-level fluctuation and $w_{0}$ is the dimensional melt velocity at the top of the column. Thus the admittance is proportional to the pressure fluctuation induced by sea level relative to static pressure over the dry melting region (this quantity is typically very small). However, it is also multiplied by $w_{0}/W_{0}$, the ratio of melt velocity at the top of the column to the mantle upwelling velocity, which is typically fairly large. Thus the admittance can be significant. The admittance of porosity is smaller than that of melt flux by a factor of $\tfrac{1}{n}$, since we lose the factor coming from $\tfrac{d\overline{Q}}{d\overline{\phi}}$.

The admittance of the carbon flux $A_{Q_{c}}$ can be estimated

[TABLE]

where the first term comes from the fluctuation in melt flux and the seccond term comes from the fluctuation in carbon concentration. By making the same approximations to the steady state used to derive (48), we find that the admittance of carbon flux is smaller than that of melt flux by a factor $\tfrac{n-1}{n}$. In particular,

[TABLE]

The admittance of carbon flux is less than that of the melt flux because the variation in carbon concentration partially compensates the variation in melt flux. While these predictions only apply to the short period regime, they are consistent with the wider pattern observed in figure 8 for the sea-level record over the past 800 kyr. Admittance of the porosity is half that of the melt flux which is about double that of the carbon flux, consistent with the above theory when $n=2$.

Moreover, by calculating the imaginary part of the fluctuations, we can also estimate the phase of the fluctuation. Finally, by comparing the magnitude of the real and imaginary parts, we can estimate the critical frequency below which this large $\omega$ regime no longer applies. Physically, this allows us to estimate when the fluxes are proportional to sea level and when they are proportional to the rate of change of sea level. Both of these depend only on the ratio of the scaled imaginary part of the fluctuation to the real part. Approximate formulae for these are

[TABLE]

Note that equation (51a) also applies to melt flux, since $\hat{Q}_{i}/\hat{Q}_{r}=\hat{\phi}_{i}/\hat{\phi}_{r}$ . In deriving equation (51b), we assumed that $D\ll\overline{\phi}\Gamma^{*2}/\mathcal{M}F_{\mathrm{max}}$ . Finally, we can use equation (51) to infer the critical period above which the phase shift is significant. We find a dimensional critical period for the melt flux and carbon fluxes, respectively,

[TABLE]

where $C\approx 1$ is a dimensionless prefactor (which can be chosen to match a specific phase shift). For the parameters given in table 1, the critical period for the melt flux is about 23 kyr, consistent with figure 5(b). Note that ${H_{\mathrm{dry}}}/{w_{0}}$ is a measure of the transit time of melt across the dry melting region. The transit time based on the true $z$-dependent velocity is (approximately, when $\overline{\phi}\ll 1$) a factor of $n$ larger than the time based on the maximum melt velocity $w_{0}$. If the forcing period is much less than the melt transport time, there is little melt segregation over the forcing cycle and instead the porosity fluctuation and carbon concentration respond near instantaneously to the change in sea level.

Appendix D Lateral melt focusing in a pseudo-two-dimensional melting region

Following Langmuir et al. (1992), we consider a triangular melting region with a base at $z=0$ and apex at the ridge axis, $z=1$. As in the main text, all lengths are non-dimensionalized by $H$ and fluxes by $W_{0}$. The melting region is capped by a decompaction channel with dip $\alpha$ (Sparks and Parmentier, 1991), sketched in figure 1. At each distance $x$ from the axis, a melting column spans $0\leq z\leq z_{c}(x)$, where $z_{c}(x)\equiv 1-x\tan\alpha$ is the depth to the decompaction channel as a function of distance $x$ off-axis. All melting columns (except the one at $x=0$) empty into the decompaction channel, which focuses magma laterally up-dip to the ridge axis. Only magma that enters the decompaction channel within some maximum focusing distance $x_{f}$ actually arrives at the ridge axis (e.g., Katz, 2008; Hebert and Montési, 2010). Following Burley and Katz (2015), we choose $x_{f}$ to give a mean crustal thickness of 7 km at the ridge axis.

To keep things simple (and consistent with the triangular geometry), we assume a uniform upwelling rate at the bottom of the melting region for all $x$. Hence all columns are identical except for their height $z_{c}$. The rate $Q_{\text{mor}}$ at which magma arrives at the ridge axis is

[TABLE]

where ${\tau}(x)$ is the travel time required for melts that enter the decompaction channel at a distance $x$ to arrive at the ridge axis. The factor of two comes from the contributions from both sides of the axis, assuming mirror symmetry about $x=0$. Note that $Q_{\text{mor}}$ is a flux multiplied by a length, so the dimensional equivalent needs to be multiplied by a factor of $W_{0}H$. We use the same method for the carbon flux.

Then, writing $dx=-d{z}/\tan\alpha$, we express the integral as

[TABLE]

Note that $\tau(z)=\tau(z_{c}(x))$, where $z_{c}(x)$ is the top of the melting column at position $x$. It is understood that regardless of the argument ($x$ or $z_{c})$, $\tau$ represents the transit time of magma laterally along the sloping decompaction channel.

The choice of ${\tau}$ closes the model. We consider two cases. First, we take $\tau=0$ to represent instantaneous focussing as in Burley and Katz (2015). Second, we follow Crowley et al. (2015) in using the steady-state column transport time to define the focusing transport time. In particular, we assume that the focusing time of melt that enters the decompaction channel is the same as the time that would be required for melt to continue up the melting column to the top. Thus

[TABLE]

For either case, we use the decomposition ([math]) to separate the MOR magma delivery rate into steady and fluctuating parts as

[TABLE]

where $z_{f}\equiv z_{c}(x_{f})=1-x_{f}\tan\alpha$.

Note that the dimensional mean crustal thickness is given by

[TABLE]

where $\rho_{l}=2800$ kg m*-3* and $\rho_{c}=2900$ kg m*-3* are magma and crustal density, respectively.

The admittance and lag of the pseudo-2d models is computed identically to that of the column model, except in using $Q_{\text{mor}}$ from equation (56). Carbon fluxes are treated in the same way.