Variations in mid-ocean ridge CO2 emissions driven by glacial cycles

The geological record shows links between glacial cycles and volcanic productivity, both subaerially and at mid-ocean ridges. Sea-level-driven pressure changes could also affect chemical properties of mid-ocean ridge volcanism. We consider how changing sea-level could alter the CO2 emissions rate from mid-ocean ridges, on both the segment and global scale. We develop a simplified transport model for a highly incompatible element through a homogenous mantle; variations in the melt concentration the emission rate of the element are created by changes in the depth of first silicate melting. The model predicts an average global mid-ocean ridge CO2 emissions-rate of 53 Mt/yr, in line with other estimates. We show that falling sea level would cause an increase in ridge CO2 emissions with a lag of about 100 kyrs after the causative sea level change. The lag and amplitude of the response are sensitive to mantle permeability and plate spreading rate. For a reconstructed sea-level time series of the past million years, we predict variations of up to 12% (7 Mt/yr) in global mid-ocean ridge CO2 emissions. The magnitude and timing of the predicted variations in CO2 emissions suggests a potential role for ridge carbon emissions in glacial cycles.


Introduction
Glacial cycles transfer ∼ 5 × 10 19 kg of water between the oceans and ice sheets [Tushingham and Peltier, 1991], leading to accumulation and ablation of kilometres of ice on the continents and sea-level change of ∼ 100 m. In Iceland, for example, the pressure change associated with melting of the ice sheet since the last glacial maximum had well-documented consequences for the volcanic activity [Sigvaldason et al., 1992, Jull and McKenzie, 1996, Maclennan et al., 2002 and lava geochemistry [Maclennan et al., 2002]. More broadly, continental volcanism in both the northern and southern hemispheres shows increased activity associated with the last deglaciation [Gardeweg et al., 1998, Jellinek et al., 2004, Nowell et al., 2006, Huybers and Langmuir, 2009]. Huybers and Langmuir [2009] and Lund and Asimow [2011] hypothesised that the pressure variations caused by changing sea level during glacial cycles would affect mid-ocean ridge (MOR) volcanism. Crowley et al. [2015] documented variations of bathymetry near the Australian-Antarctic ridge that are possible evidence of such glacial effects.
A simple argument shows that crustal thickness and bathymetric variations resulting from sea-level variation should be expected. The melting rate of a parcel of mantle beneath a MOR is proportional to its depressurisation rate. As the parcel upwells, it depressurises due to the decreasing height of rock above it. The rate of change of pressure due to upwelling is the mantle density times the upwelling rate (∼ 3 cm/yr). Sea-level variation can modify this depressurisation rate; the pressure change due to varying sea level is the product of water density and the rate of change of sea level (up to 100 m in 10 kyr or 1 cm/yr). Mantle density is about three times greater than water density, so sea-level changes can modify melting rates by up to ±10%. Crowley et al. [2015] apply a paleo-sea-level reconstruction to a simulation of MOR melting and melt transport, predicting variations in crustal thickness and bathymetry consistent with observations. Given this evidence for glacial cycles affecting MOR melt production, it is reasonable to consider if, as in Iceland [Maclennan et al., 2002], the chemistry of the lavas is also affected. To investigate this we develop a model of the transport of a highly incompatible element from the asthenosphere through the melting region to the MOR. Highly incompatible elements partition strongly into the melt, rather than remaining in the residual solid. Approximating this as complete incompatibility creates useful simplications in modelling. For instance, a completely incompatible element's path through the melting region is entirely determined by the motion of the melt, without any need to consider that element in the solid. Furthermore, to leading order, for small perturbations to the melting rate (such as those caused by sea-level change), the mass-flow rate of the element through the MOR is constant.
In a simple model of melting beneath a mid-ocean ridge, a parcel of mantle upwells adiabatically beneath the ridge axis, cooling slightly due to its expansion. The pressure-dependent solidus temperature of that parcel decreases as it ascends; at some depth (or, equivalently, pressure), the temperature of the parcel is equal to its solidus temperature. This depth, thought to be around 60 km, is called the depth of first silicate melting. With further upwelling, the parcel's temperature exceeds the solidus, and it partially melts. As soon as the first increment of melt is produced, 100% of the completely incompatible element that was present in the solid mantle is transferred to the melt. Because the mantle is permeable and the melt is less dense than the residue, the melt ascends faster than the solid, segregating from its source. Melt segregation and transport of incompatible elements thus begins at a depth that is pressuredependent. Variations in sea level will cause this depth of first silicate melting (and initiation of melt segregation) to rise and fall. The rate at which mantle crosses this boundary and delivers its content of incompatible elements to the melting region is the mantle upwelling rate minus the rate of upward motion of the boundary. So in this model, with mantle upwelling rate assumed to be constant, variations in the depth of first silicate melting cause variations in the flux of an incompatible element. As sea level falls, the depth of first melting increases, upwelling mantle crosses into the melting region faster, and the flux of incompatible element increases; the reverse is true for sea-level rise. Any perturbation to the melting rate within the melting region does not alter the mass of the incompatible element in the melt, it only dilutes (or concentrates) the element. Variations in melt-transport rate associated with melting perturbations are a secondary effect and are not considered here.
A more nuanced view of melting may disagree with this simple story in some of the details, especially with the inclusion of volatile elements that are present in small concentrations in the mantle. Experimental evidence suggests that CO 2 -rich melt forms at ∼ 250 km depth  and has a low viscosity that rises sharply with silica content at shallower depth [Kono et al., 2014]. If such carbon-rich melts can segregate from the solid mantle it would complicate the role of the transition to silicate melting at around 60 km. However, it remains an open question whether oxygen fugacity allows such melts to form and, if they do form, whether such tiny melt fractions can segregate from the solid mantle.  suggests carbonatite melt fractions reach ∼ 0.03 wt% deep in the mantle below ridges, which is at the lowest limit of carbonatite melt interconnectivity of 0.03-0.07 wt% in ∼ 0.05 mm olivine grains [Minarik and Watson, 1995]. The additional presence of water might increase the melt fraction to 0.06-0.1 wt%  by 150 km depth, but the threshold for interconnectivity of such melts is not known. In our calculations, we assume that these melt fractions do not segregate from the solid mantle until the onset of volatile-free peridotite melting at ∼ 60 km raises the melt fraction, creating an interconnected, permeable network of pores.
Among the highly incompatible elements. we focus on carbon, despite its active role in the thermodynamics of melt production, because variations in CO 2 emissions from the solid Earth are potentially important to understanding past variation of the climate. The solid Earth contains 10 10 -10 11 Mt carbon [Dasgupta and Hirschmann, 2010], orders of magnitude more than the atmosphere (0.6 × 10 6 MtC [Solomon et al., 2007]) and the oceans (4 × 10 7 MtC [Solomon et al., 2007]). Solid-Earth carbon emissions from MORs are estimated as ∼ 25 MtC/yr [Coltice et al., 2004, Cartigny et al., 2008, Marty and Tolstikhin, 1998] and from arc volcanoes ∼ 20 MtC/yr [Coltice et al., 2004]. As the largest reservoir, the solid Earth's carbon budget is known to control atmospheric carbon on multi-Myr timescales; geological ages show a correlation between volcanic activity and atmospheric CO 2 concentration [Budyko et al., 1987]. Hence there is evidence for both MOR volcanism being affected by glacial cycles, and for the effect of volcanic CO 2 emissions on atmospheric concentrations. Whilst we focus on CO 2 in our model, the same theory applies equally to other highly incompatible elements such as U, Th, Nb, Ba, and Rb.
The model is developed under the guiding principle that it should be simple enough that the connections between the assumptions and the outputs are readily traceable. The full model is assembled from independent, decoupled parts that capture the key physics with minimal complexity. Mantle flow is modelled by the passively-driven corner-flow solution [Batchelor, 1967, Spiegelman andMcKenzie, 1987]; lithospheric temperature structure and thickness is computed with a half-space cooling model [Turcotte and Schubert, 2002]. To quantify melt generation and transport we use one-dimensional compaction columns [Ribe, 1985, Hewitt, 2010 that are based on conservation of energy, mass, and momentum at steady state in a homogenous mantle. The outline of the melting region is given by a parameterised solidus [Katz, 2008]. A focusing width is applied such that melt focused to the ridge produces a maximum crustal thickness of 7 km.
A detailed discussion of the assumptions made in deriving the model is presented below. To summarise the results, the model predicts that a section of MOR spreading at 3 cm/yr half-rate will see a change in the rate of efflux of highly incompatible elements (e.g., CO 2 ) of ∼ 8% for a linear sea-level change of 100 m in 10 kyrs. For reconstructed sea level data and the present distribution of plate spreading rates, the model predicts global MOR emissions to vary from the mean by up to ±6%. These results are sensitive to the permeability of the mantle, which is a primary control on the rate of melt transport.
Section 2 details the model used to predict CO 2 emissions for a section of mid-ocean ridge. The behaviour of the model is demonstrated for simple scenarios of sea-level variation in sections 3.1-3.3, and then the model is applied to the global MOR system under a reconstructed sea-level history in sections 3.4 and 3.5. The results are discussed in section 4 and the key conclusions stated in section 5.

The Model
Our goal is to develop a method to compute the CO 2 emission rate E CO 2 (mass per unit time per unit length of ridge) from a segment of ridge. To achieve this we require a model of CO 2 flux into the melting region and also of its transport to the ridge. We approximate the behaviour of CO 2 as perfectly incompatible, and hence there is no exchange of CO 2 between phases during melt transport. The rate of ridge emission of CO 2 can then be quantified by integrating the mass flux into the base of the melting region f CO2 (mass/area/time) and using the travel time from the base of the melting region to the ridge. This is formulated as where x is the horizontal distance from the ridge axis, x f is the maximum distance over which melt is focused to the ridge axis, U 0 is the half-spreading rate, and t is time. A parcel of melt arriving at the ridge axis at time t was produced by mantle that crossed into the melting region at time t s . The travel time of melt from the base of the melting region to the ridge is represented as τ , which varies with lateral distance x from the ridge axis. Hence the source time is t s (x) = t − τ (x). The factor of two in eqn.
(1) arises from the symmetry of the melting region across the ridge axis. A sketch of half the melting region is shown in figure 1. Equation (1) requires an expression for f CO2 . This is a product of the rate at which mantle material crosses the depth of first silicate melting z m and the CO 2 concentration in that material. For generality, we consider an expression that allows for volatile-enriched insipient partial melting beneath the depth of first silicate melting, though we will later exclude this scenario from consideration. Thus, f CO 2 is written as where W m (x, U 0 ) is the upwelling rate of the mantle, w m is the upwelling rate of insipient melt, C CO 2 is a mass concentration of CO 2 , and φ is the volume fraction of melt; all of these are evaluated at depth z m and distance x from the ridge axis. The first term in parentheses on the right of the equation is the rate at which mantle crosses into the melting region. The concentrations and melt fraction may be considered steady-state, constant values as long as w m ≥ W m ≥ max (ż m (t)). This ensures that the rate at which material crosses the depth of first silicate melting is always positive or zero. Even the fastest sea-level changes on record, meltwater pulses during the last deglaciation, satisfy these conditions for most MORs: 20 m sea level change in 500 years givesż m = 1.3 cm/yr. However, meltwater pulse events are not resolved in the reconstructed sea-level series that we consider in this paper [Siddall et al., 2010], so the conditions w m ≥ W m ≥ max (ż m (t)) are true with only occasional exceptions for the slowest spreading ridges. With these conditions satisfied and assuming that either φ(z m ) = 0 or that w m = W m , equation (2) can be simplified to Here we should interpret C CO 2 as the mass concentration of CO 2 in the solid mantle plus comoving insipient melt, if it is present. This model could be modified to accomodate more complicated situations, but this is not given further consideration at present. The flux of CO 2 in equation (3) depends on the solid mantle upwelling rate at the depth of first silicate melting. Approximating this as passive (plate-driven) flow of isoviscous rock, mantle upwelling is given by the corner flow solution [Batchelor, 1967, Spiegelman andMcKenzie, 1987]. The vertical component of this solution, evaluated at z = z m , can be written as where the lithosphere is represented as a wedge with angle α c to the horizontal. We follow Spiegelman and McKenzie [1987] in computing the wedge angle to approximately match the plate thickness at a specified distance from the axis. We use α c = tan −1 (z l /x w ), such that the wedge intersects the upper boundary of the melting region z l at the maximum width of melting region x w (see appendix A.3 for details). This definition of the wedge angle ensures that the lithospheric wedge does not overlap with the melting region, and that the upwelling rate is small but non-zero at the extreme width of the melting region (0 ≤ W m U 0 ). An expression for the depth of first melting z m is needed in equations (3) and (4). Taking decompression as the only influence on local mantle temperature prior to melting, we model the depth of first melting as the intersection between an adiabatic temperature profile and the solidus temperature profile. We approximate both profiles as linear with respect to depth (details in appendix A.1) to get where T is the mantle potential temperature, T S ref is the solidus temperature at reference mantle composition and surface pressure, γ is the Clausius-Clapeyron slope for the mantle, ρ is the mantle density, g is gravitational acceleration, α is the coefficient of thermal expansion, c is the specific heat capacity, S is the sea-level deviation from a long-term mean, and ρ w is the density of water. The first term in equation (5) is the dry peridotite melting depth of ∼ 60 km, and the second term is the shift in melting depth due to sea-level. The only time dependent variable in equation (5) is the sea level S, so differentiating gives These equations state that silicate melting begins at a fixed pressure, but the depth corresponding to this pressure varies as sea level rises and falls.
The source time t s used in equations (1) and (3) depends on the travel time of melt to the ridge from any point x along the base of the melting region. For simplicity, we consider melt flow as following a vertical path from the base to the top of the melting region, then following a high porosity channel along the impermeable top of the melting region to the ridge axis, as illustrated by streamlines τ 1 and τ 2 in figure 1. This is a reasonable approximation of numerically modelled streamlines for homogenous mantle [e.g., Katz, 2008], where buoyancy forces drive vertical fluid flow in the majority of the melting region, with the compaction pressure only becoming large enough to deflect melt flow from the vertical within a few kilometres of melt-impermeable boundaries [Sparks and Parmentier, 1991]. We consider flow along the high porosity channel as instantaneous, motivated by the high flow rates expected there, as compared to vertical flow rates in the rest of the melting region [Katz, 2008]. To compute the travel time from the base of the melting region to the base of the lithosphere, we use a 1D compaction column from Hewitt [2010]. This model assumes Darcy flow and thermodynamic equilibrium for a two-component, homogenous mantle with a constant Clausius-Clapeyron slope. Following the reduced model of Crowley et al. [2015], we assume that small variation in melting rates due to sea-level change do not significantly affect melt velocity or travel time. Furthermore, in computing τ we take z m as constant, because changes in z m due to sea level are only tens of metres, changing τ by 1%. This gives a travel time where η f is the mantle melt viscosity, K 0 is the permeability of the mantle at reference porosity φ 0 , ∆ρ is the density difference between the solid and melt, n is porosity exponent in the permeability relation (n ≈ 2.7, Miller et al. [2014]), z l (x, U 0 ) is the depth of the upper boundary of the melting region, and Π is the adiabatic melt productivity (kg of melt produced per m 3 per m upwelling). This productivity is given by Hewitt [2010] as a ratio of thermodynamic parameters. To close the model we need an expression for the upper boundary of the melting region z l (x). This boundary is located where the local temperature equals the solidus temperature, which is controlled by the thickness of the conductively cooled boundary layer that forms the lithosphere. Crowley et al. [2015] also included the effects of adiabatic decompression and latent heat removal, but this leads to an expression for x l (z) that cannot be inverted. However, as shown in Appendix A.2, the depth of this melting boundary is approximated by an isotherm of the half-space cooling model, which can be expressed as z l (x). This isotherm is hotter than the low-pressure solidus, but the temperature difference compensates for the change in solidus temperature with respect to depth. We use where T l is the temperature of the upper boundary of melting region, assumed to be constant. T 0 is the temperature of the ocean floor, and κ is the thermal diffusivity of the mantle. Note that although we have neglected adiabatic decompression and latent heat of melting in equation (8), these are accounted for in the melting calculations.

Results
Below we demonstrate the behaviour of the model by a series of examples, then consider global CO 2 emissions for reconstructed sea level. We begin by modelling a unit length of mid-ocean ridge in the absence of sea-level change. This defines the baseline state of the model that we compare against when applying sea-level forcing. The first dynamic example is a single, linear sea-level change, which illustrates key characteristics of the model and emphasises the importance of the melt travel time. The emissions curve for this example approximates the Green's function for the model: the response to a step-change in sea level. We next consider simple, periodic sea level curves, demonstrating the system's response to an oscillatory forcing as a function of the frequency of that forcing. We then compute the model's response to a reconstructed Pleistocene sea-level record. Building on these calculations for a single section of mid-ocean ridge, we calculate predictions for a global model, created by summing sections with appropriate spreading rates over the full MOR system. We demonstrate this composite model by considering the global emissions response to a single sea-level change and then to the reconstructed Pleistocene sea-level curve.

Constant Sea Level & Baseline Emissions
Baseline emissions are defined as the steady-state emission rate of CO 2 from the MOR for constant sea level. This represents the background state, which is disturbed after sea-level change. Figure 2 shows baseline emissions as a function of plate spreading rate and demonstrates that, in the absence of changing sea level, CO 2 emissions per metre of ridge are approximately proportional to the half-spreading rate U 0 . A faster spreading rate drives faster mantle upwelling, bringing more CO 2 into the melting region. Faster spreading rates also increase the width of the melting region, leading to more CO 2 (and more melt) being focused to the ridge axis. However, not all melt produced is focused to the MOR; at some lateral distance, melts are frozen back into the lithosphere, rather than travelling to the ridge axis [e.g., Katz, 2008]. We incorporate this detail by enforcing a maximum focusing width x f such that melt focused to the MOR will produce crust ≤ 7 km thick in steady state, consistent with observations [White et al., 1992, Bown andWhite, 1994]. Crustal thickness is calculated as the volume flow rate of melt to the ridge axis times the density ratio of basaltic melt to oceanic crust, divided by the half-spreading rate. At half-rates U 0 ≥ 1 cm/yr, the focusing width is smaller than the width of the melting region. With the imposition of this limit on melt focusing, there is a slight change in the slope of the baseline emission curve at U 0 = 1 cm/yr (Fig. 2).
For the range of U 0 on Earth of up to 8 cm/yr, figure 2(a) shows baseline emissions of up to 2300 kg m −1 yr −1 . This is calculated using 0.38 kg CO 2 per m 3 of upwelling mantle (125 ppm CO 2 by weight), based on an average of MOR source mantle of 50-200 ppm CO 2 by weight [Dasgupta, 2013, Cartigny et al., 2008, Saal et al., 2002, Marty and Tolstikhin, 1998].
Subsequent sections present scenarios of changing sea-level. In these sections, plots depict emissions in terms of percentage difference from baseline emissions for the appropriate halfspreading rate; these percentage values apply to all highly incompatible elements. Percentage results can be converted to CO 2 mass via figure 2. The deviations from the baseline emission rate are a consequence of non-zeroż m in equation (3). By definition, total emissions are equal to this deviation summed with baseline emissions.

Single Sudden Change in Sea Level
Imposition of a sudden sea-level change exposes the behaviour of the model. We use a steep, linear ramp, which gives a box function inṠ. Figure 3 shows the predicted MOR CO 2 emission rate E CO 2 resulting from this sea level forcing. CO 2 emissions remain at the baseline level for ∼90 kyrs after the change in sea level, then there is a sharp rise in E CO 2 , representing an 8% increase in emission rate. This delayed response is due to the travel time of CO 2 from the base of the melting region to the ridge. After the 8% peak, E CO 2 falls sharply to about 3%, followed x w x f focusing width x f and full width x w of the melting region; all shown for varying half-spreading rate. Focusing width is equal to the width of the melting region when crustal thickness is less than 7 km (U 0 ≤ 1 cm/yr), and is otherwise limited such that crustal thickness does not exceed 7 km. The switch between these behaviours is marked by the grey dotted line at 1 cm/yr. E CO 2 is computed with CO 2 concentration in the mantle of 0.38 kgm −3 (125 ppm CO 2 by weight). by a slow decay until 130 kyrs, after which there is a linear drop back to the baseline level over the duration of the box-pulse inṠ. The origin of this emissions pattern in figure 3(c) is explored in figure 4. Figure 4 demonstrates how the shape of the emission response curve in figure 3(c) is a consequence of the travel time τ and its variation with respect to distance x from the ridge axis. This travel time, shown in figure 4(c), ranges between τ min for melts originating in the mantle beneath the ridge axis and τ max for distal melts that are focused laterally. The mixture of melts that arrives at the ridge at time t contains CO 2 that was transported in melts that initiated along the base of the melting regime. These deep melts formed and began to segregate from their source at times in the past t − τ (x). The CO 2 content of the segregated melt is different to the baseline case according to theṠ value at t − τ (x). Therefore, we can calculate the deviation from the baseline by considering the rate of sea-level change (alternatively,ż m ) acting on the melting region when each element of melt first segregated. This is represented in NegativeṠ is plotted so that peaks inṠ and consequent peaks in E CO 2 point in the same direction.
figure 4(c), and the integral of this plot with respect to x is directly proportional to E CO 2 in figure 3. The E CO 2 response is clarified by again considering a sharp, linear drop of sea level that occurs over the interval 0 ≤ t ≤ 10 kyrs. The drop in sea level causes the depth of the base of the melting region to increase, importing additional CO 2 into the melting regime. In panels (a-i), (a-ii), and (a-iii) of fig. 4 we see the box pulse receding into the past as time progresses over t 1 < t 2 < t 3 (figure 3 shows the emission-rate response at these times). The start and end times of the pulse are projected onto τ (x) in panel (b). At t 1 , the projection lines do not intersect τ (x), meaning that the CO 2 perturbation generated byṠ has not yet reached the ridge axis. Therefore the emissions curve in figure 3(c) remains at the baseline level. At t 2 , the projection lines span τ min ; the shallow slope of τ (x) near τ min means that CO 2 from a broad (30 km) region affected by theṠ pulse is arriving at the ridge axis at t 2 . This causes the spike in emissions shown in figure 3(c). At t 3 , the interval of sea-level change has receded far into the past. Distal melts from the base of the melting region at 50 to 58 km off-axis are arriving at the ridge axis. The narrowness of this band translates to a reduced (but non-zero) value of E CO 2 at t 3 in fig. 3(c). As the sea-level drop recedes further into the past, the emission rate drops to zero because the very distal melts (from x > x f with τ > τ max ) are not focused to the ridge axis.
In the limit of vanishing duration of sea-level change, the pulse ofṠ approximates to a Dirac delta function. A sum of shifted, scaled delta functions ofṠ can be used to approximate any sea level signal. The associated E CO 2 response of that signal can be represented as an appropriately weighted sum of the E CO 2 response to a delta function. Hence the E CO 2 response shown in figure 3(c) is an approximation of the Green's function.

Oscillating Sea Level
We now consider oscillatory sea level and discuss the concepts of lag and admittance. Figure 5 shows a pair of oscillating sea level scenarios and their predicted E CO 2 variation. The left column (i) shows a time series of alternating box-pulses inṠ; the right column (ii) shows a sinusoidal sea-level variation. Figure 5(c-i) has an oscillating series of peaks in E CO 2 resulting from a series ofṠ box-pulses in Figure 5(b-i). ThisṠ series is equivalent to summing single box-pulses from figure 3(b), with suitable offset and amplitude. Similarly, the CO 2 emissions can be represented as a sum of offset emissions spikes from single, linear changes in sea level. The train of emissions peaks in figure 5(c-i) and (c-ii) stabilises in amplitude after t ≈ τ max ; this transient represents the spin-up time of the model, associated with the tail of excess emissions shown in fig. 3. Figure 5(ii) shows sinusoidal sea level and provides the context to define the lag metric. Lag L is the time between a peak inṠ and the corresponding peak in emissions. Because the time interval around τ min kyrs before t has the largest influence on E CO 2 at t, the E CO 2 signal should lagṠ by about τ min . However, the lag will not be exactly τ min , as the influence ofṠ on E CO 2 is felt up to τ max years after the change in sea level. Thus the exact value of lag will be slighly Considering the left column, each box-pulse inṠ produces an emissions peak/trough as in figure 3(c), with interference where they overlap, giving the net effect seen in (c-i). The stready upward shift in the E CO 2 -peaks until ∼ 150 kyr is due to this overlap. If the first box-pulse inṠ had been in the opposite direction, E CO 2 peaks in (c) would be mirrored across E CO 2 = 0. greater than τ min , and we expect this difference to depend on the period of sea-level oscillation relative to τ max . In particular, when the period approaches and exceeds 2(τ max − τ min ), the lag becomes equal to the mean melt-transport time τ mean . Figure 6 shows lag, τ min , τ mean , and τ max for varying half-spreading rate, permeability constant, and sinusoidal sea-level period. We note that τ min ≤ L ≤ τ mean and therefore we assume L ≈ τ min with a small, systematic error. Sinusoidal variation of sea level also provides a context in which compute admittance. Admittance is the ratio of the response amplitude to the forcing amplitude as a function of the sinusoidal forcing period. We define two versions of admittance: absolute admittance, with units of kilograms CO 2 per metre of ridge per 100 m of sea-level change per year, and relative admittance, with units of percentage change from baseline emission rate per 100 m of sea-level change. The latter is absolute admittance divided by the baseline emission rate. Figure 7 shows absolute and relative admittance and how they vary with changing sea-level period, halfspreading rate, and permeability constant. We discuss both the trends and the oscillations of these curves, starting with absolute admittance. Absolute admittance (panels (a) and (b) of figure 7) depends on the period of sea-level oscillation, the permeability, and to a lesser extent, on the half-spreading rate. We consider these in turn. Shorter periods of sea-level variation at constant amplitude give larger values of |Ṡ| and |ż m |, and hence increase the temporal variation of f CO 2 . This causes increased deviation of ridge emissions from baseline. Increased mantle permeability and spreading rate both reduce the melt travel-time from the base of the melting region. In the melt travel-time calculation (eqn. (7)), the permeability appears as K 0 , while the spreading rate is directly proportional to the mantle upwelling rate term, W m . A reduced melt travel-time overall implies a smaller difference between τ max and τ min , and therefore a focusing of CO 2 from the base of the melting region to the ridge axis over a shorter interval in time. The uncertainty in mantle permeability translates to a broader spread of the absolute admittance curves than does the range of halfspreading rates considered.
Relative admittance is equal to the absolute admittance normalised by the baseline emissions rate. The baseline depends on half-spreading rate but not on permeability ( fig. 2(a)). We therefore see a difference between absolute and relative admittance in figure 7(b) and (d), respectively. For slow-spreading ridges, which have a low baseline emissions rate, the normalised variance (and hence the relative admittance) is larger.
The oscillations superimposed on the primary admittance trend are not physically significant, but are readily explained. They arise from variation in the number of sea-level half-cycles that fit into the time interval from t − τ max to t − τ min ; this is the time interval over which Ṡ (t − τ ) can contribute CO 2 emissions at time t. For oscillatory sea level, each positive or negative peak inṠ has an opposing emissions effect relative to the prior negative or positive peak. If there is an unmatched peak affecting the bottom of the melting region, the amplitude of E CO 2 variations is larger and the admittance is higher.
Broadly, the patterns of admittance and lag imply that the dominant sea-level changes, in terms of CO 2 emission variation, will be those with large, short period changes in sea level. The emissions variation associated with such changes will lag the forcing by approximately τ min . The modelled magnitude and lag of E CO 2 changes are affected by both K 0 and U 0 .

Reconstructed Pleistocene Sea Level
The simple scenarios of sea-level variation presented above give insight into the behaviour of the model, but are not necessarily representative of the variations that have occured naturally, over the past two million years. We move, therefore, to a model forced by the time-series of reconstructed global sea level from Siddall et al. [2010], shown in figure 8(a). Other reconstructions exist, but the differences between them are small enough that we follow Crowley et al. [2015] and consider only this one. Siddall et al. [2010] record data every 3 kyrs and, based on their reconstruction, the highest rates of sea-level change ( fig. 8(b)) meet the condition max (ż m (t)) < W m required for validity of equation (3). The result of applying reconstructed sea level to ridge E CO 2 is shown in figure 8(c). Figure 8(c) shows up to ±8 % range in E CO 2 for moderate half-spreading rate and permeability. The E CO 2 curve is, qualitatively, an offset version ofṠ with small variations smoothed outas expected from E CO 2 being approximated by convolvingṠ with the emissions response in figure 3(c). Within this framework, we now consider how to apply the model to global MOR emissions.

Global Mid-Ocean Ridges
The global MOR system is composed of ridge segments spreading at different rates, ranging from the ultra-slow Gakkel ridge to the fast East Pacific Rise. The baseline emissions depend on the half-spreading rate, as does the character of the emissions response to sea-level change. The global response to sea-level change should therefore be computed as the segment-scale response, integrated over the global MOR system, where G CO 2 is global MOR emissions of CO 2 in kg-yr −1 , l is arc length along the ridge, and L MOR is the total length of the MOR system. This integral can be approximated by discretising the half-spreading rate into bins U 0i and summing the local response in each bin. A weighting is applied to each entry in the sum to account for the total length of segments with half-spreading rates in that bin. The sum is written as where N is the total number of spreading rate intervals to sum over and L i (U 0i ) is the total length of MOR in a particular spreading-rate bin. The local emission rate E CO 2 in each bin is computed by adopting the average half-spreading rate of the bin and assuming that sea-level change is eustatic -the same for all segments globally. Gale et al. [2013] provide a catalogue of segment lengths and spreading rates for the global MOR system; the total ridge length is 61000 km with a mean half-spreading rate of 2.5 cm/yr. A histogram of these data is plotted in figure 9(a). Figure 9(b) shows, for each spreading-rate bin, the rate of baseline emissions per metre of ridge. Panel (c) then shows the product of ridge length and emissions rate per metre, giving the total emissions rate associated with each spreading rate bin. This is summed in eqn. (10) to give the total baseline global response. The global baseline emission rate thus predicted is 53 Mt CO 2 per year (12 Mt carbon per year). This can be compared to other estimates of 91 ± 45 MtCO 2 /yr [Coltice et al., 2004, Cartigny et al., 2008, Marty and Tolstikhin, 1998]. Our model assumes a sub-ridge mantle CO 2 concentration of 0.38 kg per m 3 ; the global baseline emission rate is directly proportional to this assumed concentration.
Before applying the reconstructed sea-level forcing to the weighted global emissions sum in eqn. (10), we consider the simple sea level forcing that was used to probe the behaviour of E CO 2 in fig. 3. This linear ramp in sea level is applied to global emissions in figure 10(a) for a range of mantle permeabilities. Global emissions in figure 10(a) are, unsurprisingly, more complex than the E CO 2 equivalent, as they consist of a summation of N E CO 2 peaks from figure 3, weighted according to the ridge length in each bin and offset due to the variation in τ with U 0i . Compared with figure 3, G CO 2 has a smaller percentage difference in the rate of CO 2 emissions and is spread out over a longer time-period.
Finally, we apply the reconstructed sea level time series to the global model. Figure 10 shows G CO 2 for reconstructed sea level, calculated and plotted for a set of three values of permeability constant K 0 in panels (d-i), (d-ii), and (d-iii). These curves demonstrate a reduction in G CO 2 range from 6.6 MtCO 2 /yr for K 0 = 10 −11 m 2 to 1.2 MtCO 2 /yr for K 0 = 10 −13 m 2 . The reduction in G CO 2 range occurs because K 0 affects the range of τ min that is implicit in the global sum in eqn. (10). A large value of K 0 gives higher overall permeability, shorter melttransport times, and a global range in τ min that is smaller. Therefore the emissions response to a box-pulse inṠ ( Fig. 10(a)) is temporally concentrated and attains a higher peak value. Hence, larger K 0 causes greater amplitude of variation in G CO 2 for reconstructed sea level.

Discussion
The model prediction for global MOR CO 2 emissions at constant sea level in figure 9 is 53 Mt CO 2 per year (14 Mt carbon per year), a value that is within error of estimates from analyses of mid-ocean ridge basalts of 91 ± 45 MtCO 2 /yr [Coltice et al., 2004, Cartigny et al., 2008, Marty and Tolstikhin, 1998]. However, note that our prediction is not independent of MORB studies because our choice of C CO 2 is based on observed MORB chemistry. The ranges in global CO 2 emissions under reconstructed variation of past sea level (Fig. 10) are 1 MtCO 2 /yr, 3 MtCO 2 /yr, and 7 MtCO 2 /yr for mantle permeabilities at 1% porosity of 10 −13 m 2 , 10 −12 m 2 , and 10 −11 m 2 respectively. These idealised predictions assume that 100% of CO 2 transported to the ridge axis is degassed into the oceans, rather than retained in the crust or mantle. This may, in fact, be rather accurate; Cartigny et al. [2008] estimated that over 80% of CO 2 in primitive MORB is degassed near the ridge axis. It has previously been speculated that variations in ridge CO 2 emissions would have an effect on global climate Langmuir, 2009, Lund andAsimow, 2011]; it is beyond the scope of the present manuscript to  fig. 3(a), (b)) (a) CO 2 emissions from the global MOR system; and for reconstructed sea level (b) sea-level change, (c) the negative rate of sea-level change, and (di-iii) CO 2 emissions from the global MOR system. The sea level series has been applied further back in time than shown, such that the left-most point in G CO 2 is correctly initialised by more than τ max kyrs of prior sea-level change. For (d-i,ii,iii) the lags in G CO 2 are, respectively, ∼ 60 kyrs, ∼ 120 kyrs, and ∼ 250 kyrs.
investigate this question, though it is a target for future work. Uncertainty in mantle permeability translates to uncertainty in both the amplitude of variations and the lag of global emissions from MORs. There are various experimental constraints on mantle permeability [e.g. Miller et al., 2014, Connolly et al., 2009; these tend to agree on the scaling with porosity, but disagree on the magnitude of K 0 . Furthermore, permeability is sensitive to grain size, a parameter that is poorly known for the mantle beneath MORs [although see Turner et al., 2015]. Our chosen range of K 0 is intended to accomodate these uncertainties, as well as represent an effective permeability for melt transport that may be channelised into high-porosity, high-permeability dunite channels [e.g. . Our K 0 range of 10 −13 -10 −11 m 2 encompasses a change in the amplitude of G CO 2 variation by a factor of 5, a difference in lag of 200 kyrs, and qualitative difference in the time-series of G CO 2 (Fig. 10). Therefore K 0 represents a leading source of uncertainty in the model. Uranium series disequilibria may provide an independent constraint on magma travel time from the base of the melting region [e.g., , although interpreting the various species in the decay chain is fraught with complexity. Preservation of 230 Th disequilibrium (half-life of 70 kyrs) suggests a permeability of K 0 ≥ 10 −12 m 2 .
Our model is based on the assumption that melt travels vertically from the depth of first melting to the top of the melting region and is then focused along a sloping, high porosity decompaction channel to the ridge axis [Sparks and Parmentier, 1991]. The travel time of the vertical flow is modelled by a 1D compaction column; we assume that transport in the decompaction channel is instantaneous. The systematic error introduced by the latter assumption is zero on the ridge axis (x = 0), and increases with distance x from the ridge axis. This means τ plotted in figure 4(b) is more accurate at small x, but increasingly underestimates τ for larger x. Therefore, assuming τ min is accurately modelled, τ max is probably too small, such that E CO 2 in figure 3 should have a longer tail on the right of the graph. However, long tails have little effect on the resultant E CO 2 pattern for complex sea-level changes or on the G CO 2 pattern for reconstructed sea level. Therefore the overall effect of including a finite travel time along the high porosity channel would be to make a small adjustment to the E CO 2 response. This suggests that assuming instantaneous travel time along the channel has little effect on the results of the model.
Another assumption made is that travel time is constant with respect to time, despite changes in melting rate (and thus, porosity) caused by changing sea-level. This follows the approach in the reduced model of Crowley et al. [2015], where the perturbations in porosity were taken as negligible disturbances to the travel time in a steady-state compaction column solution from Hewitt [2010].
A more significant assumption underpinning the model is that carbon dioxide behaves as a perfectly incompatible tracer -meaning that carbon concentration does not affect the mantle's physical or thermodynamic properties. However, carbon is not a trace element. In contrast, experiments by Kono et al. [2014] document the very low viscosity of insipient, carbon rich melts present at small melt fraction below the base of the silicate melting region. The experiments also show that viscosity rises sharply as the carbon is diluted by silicate melting. It would be challening to capture this variability in models, especially since the wetting properties (and hence the mobility) of carbon-rich melts are poorly constrained.
A more significant concern, however, is that treating volatiles as trace elements neglects their thermodynamic effect on melting. Small mantle concentrations of carbon affect the depth at which melting begins, though the melt fractions produced by this incipient melting are probably less than a few tenths of a percent . Our model assumes that these melts do not segregate until the onset of silicate melting. At such small porosity, it is unclear whether these carbonated melts can percolate. However, water-induced melting at the wet peridotite solidus of ∼90-120 km Langmuir, 2003, Dasgupta et al., 2013] increases the melt fraction. Again, the threshold of interconnectivity for such melts is not known, so it is possible that such deep, hydrous melts do not segregate, or do so very slowly. If the 230 Th disequilibrium observed in young MOR lavas originates with melt segregation in the presence of garnet, it would support the hypothesis of efficient segregation of hydrous melts. Overall, our model depends only on the presence of a pressure-dependent boundary that separates non-segregating melts, below, from segregating melts, above. The sharpness that is required of this boundary is unclear. Finally, our model assumes a chemically and thermally homogenous mantle, which is certainly not true of the natural system [e.g., Dalton et al., 2014]. No data exists that would allow us to accurately incorporate small-scale ( 100 km) heterogenity in model. If such heterogeneity is pervasive and at scales of ∼10 km or smaller, it would affect the style of melt transport [Katz and Weatherley, 2012], with fertile regions creating pathways for rapid melt transport through the melting region. It may be the case that this is captured by a high effective permeability, though this is probably not a testable hypothesis. Large-scale heterogeneity would leave the mantle homogenous over the scale at which we calculate E CO 2 , so the underlying melt transport model would be unaffected, though parameters would need to be adjusted according to the oceanic region. It seems likely that such variations would cancel in the integral for global CO 2 emissions rate.

Summary
The model presented above builds on the reduced model of Crowley et al. [2015] to calculate the efflux of a highly incompatible chemical component from a mid-ocean ridge, and how that efflux would vary with changes in sea level. It is based on a description of melt transport through a homogenous mantle and assumes perfect incompatibility of the component. This leads to a simple but physically consistent model of chemical transport through the melting region beneath a ridge. The model calculates total melt supply rate and global background emissions of CO 2 that are consistent with data and prior estimates [White et al., 1992, Bown and White, 1994, Coltice et al., 2004, Cartigny et al., 2008].
In the model, changing sea level affects the depth of first silicate melting; this alters the rate at which CO 2 enters the melting region, segregates from its mantle source, and (some time later) arrives at the ridge axis and is degassed into the ocean.
The emission rate of CO 2 is predicted to vary by up to 7 MtCO 2 /yr when the model is forced with reconstructed Pleistocene sea level variation. There is uncertainty in the predicted magnitude and timing (relative to sea level forcing) of this effect, as two parameters of the model -C solid CO 2 and K 0 -are weakly constrained by existing data. However, within reasonable ranges for model parameters, the amplitude of global CO 2 emission-rate variation will remain on the order of several MtCO 2 /yr. The total difference in the mass of CO 2 emitted from MORs during sea-level driven deviations from global baseline CO 2 is up to ∼80 Gt CO 2 for high permeability (see figure 10(d-i)) and much lower for lower permeability.
To our knowledge, this is the first study to attempt to quantify how CO 2 emissions are modulated by changing sea level. Our results indicate that the CO 2 emissions from mid-ocean ridges are temporally variable in response to sea-level change. These results align with the hypothesis of Huybers and Langmuir [2009]; however, whereas they assumed instantaneous melt transport, we show that the mid-ocean ridge CO 2 emissions response lags the sea-level forcing by approximately the travel time of CO 2 through the melting region. Further work is needed to assess the climatic impact of the variable CO 2 emissions predicted for mid-ocean ridges in this paper and for subarial volcanoes by Huybers and Langmuir [2009].

Value
Parameter description c 1200 J/kg K Specific heat C CO 2 0.375 kg/m 3 CO 2 mass per m 3 of mantle (derived from 125ppm by weight) g 10 m/s 2 Gravitational acceleration K 0 10 −13 -10 −11 m 2 Permeability at 1% porosity where symbols are as defined in the main text and we have added a pressure term for varying sea level. The mass change of a water column is ρ w S, provided the thermal expansion of the ocean is much less than S. Fortunately, oceanic thermal expansion is negligible, causing less than 1% of total glacial sea-level change [McKay et al., 2011]. Freshwater density is used for ρ w , as salt remains in the ocean when water is removed, and ocean water inputs are fresh. Coordinates are upwards-positive with origin at the MOR, therefore z is negative and increasing sea level S is positive. We model the solidus temperature following Katz [2008], Hewitt [2010], using a solidus that is linear in pressure and composition parameter. We consider pressure terms due to the mass of mantle above the solidus, and sea level deviations from reference conditions. This gives T Solidus = T S ref − γg (ρz − ρ w S) .
We ignore bulk ocean+atmospheric weight and topographic relief as these are accounted for in T S ref and isostatic compensation respectively. Setting the temperatures in equations (12), (11) equal and solving for depth we get Eqn. (5) in the main text.

A.2 Upper Boundary of the Melting Region
The upper boundary of the melting region, like the lower boundary, is a region where the solidus temperature matches the local mantle temperature. Local mantle temperature in the melting region can be modelled as the adiabatic temperature profile (Eqn. (11)) plus terms for a) thermal energy lost to latent heat during melting of (z − z m )ΠL/ρc, and b) conduction near the cold ceiling (i.e., the ocean floor) from the half-space cooling model. This approach assumes superposition of temperature fields without cross terms correcting for deviations from the assumptions of each model (e.g., half-space cooling assumes a constant background temperature with respect to depth). Matching this to solidus temperature from equation (12) gives where L is the latent heat capacity of the mantle. Equation (13) cannot be rearranged to give z in terms of x. However we can get x in terms of z, as shown in equation (14), an expression shown to accurately match melting region from numerical studies in Crowley et al. [2015] (justifying the previous superposition of temperature fields assumption). This gives us a reference to benchmark against: To get z in terms of x, we discard all depth-dependent temperature terms in equation (13) except half-space cooling. This approach is equivalent to plotting an isotherm in a half space cooling model. Thus Figure 11 compares equation (15) to the benchmark (eqn. (14)). For a boundary temperature of ∼1550 K -close to the solidus temperature of mantle -equation (15) is a poor approximation of the real melting region. However unphysically high boundary temperatures of 1640 K give good approximations of the melting region, as this compensates for the temperature effects not explicitly accounted for. Therefore, we use equation (15) with T l = 1643 K to approximate the upper boundary of the melting region.

A.3 Width of the Melting Region
Equation (14) shows a finite width of the melting region defined by dx l /dz = 0 (see figure 11 at z 55, x 130 km). The derivative of equation (14) is where: B = γρg − αg T /c − ΠL/ρc  [Crowley et al., 2015], compared to half-space cooling isotherms for 1565K (a realistic solidus temperature) and 1643K (a best fit isotherm). The 1643K isotherm is used as the upper boundary of the melting region for this paper. Calculated for mantle potential temperature of 1648K, surface temperature of 273K and a half-spreading rate of 3cm/yr. Plotted quantities all scale similarly with U 0 , thus isotherms have the same relative error across all spreading rates.
We have not found an analytical solution to equation (16) for dx l /dz = 0. However, computational investigation of equation (14) for varying U 0 shows that the depth of the widest point of the melting region is constant. Calculating x l at this fixed depth gives the maximum width of the melting region. Thus, where z w is the depth of the widest point of the melting region.