Deglaciation-enhanced mantle CO 2 fluxes at Yellowstone imply positive climate feedback 1 2

12 The generation of mantle melts in response to decompression by glacial unloading has been 13 linked to enhanced volcanic activity and volatile release in Iceland1 and in global eruptive 14 records2,3. However, it is unclear whether this process is also important in magmatically-active 15 systems that do not show evidence of enhanced eruption rates. For example, the deglaciation of 16 the Yellowstone ice cap did not observably enhance volcanism4, yet Yellowstone may still have 17 released large volumes of CO2 to the surface due to the crystallization of melts at depth. Here we 18 develop models to simulate mantle melt production and volatile release associated with the 19 deglaciation of Yellowstone and Iceland. In agreement with previous work1, we find mantle melt 20 production in Iceland is enhanced 33-fold during deglaciation, generating an additional 3728 km3 21 of melt and releasing an additional 31–51 Gt of CO2. Beneath Yellowstone, we find mantle melt 22 production is comparably enhanced 19-fold during deglaciation, generating an additional 815 23 km3 of melt, though thicker lithosphere may prevent the transport of this melt to the surface. 24 These melts segregate an additional 135–230 Gt of CO2 from the mantle, representing a ~23– 25 39% increase of the global volcanic CO2 flux (if degassed during deglaciation). Our results 26 suggest deglaciation-enhanced mantle melting is important in continental settings with partially 27 molten mantle (potentially Greenland and West Antarctica) and may result in positive feedbacks 28 between deglaciation and climate warming. 29 30 Main 31 As an ice mass retreats and unloads the Earth’s surface, the underlying mantle rebounds 32 and undergoes a reduction in pressure. If the mantle is above the solidus, this decompression 33 generates additional melting relative to any background rate. Enhanced mantle melting can result 34 in increased volcanic activity1,5, which in turn may incite the release of aerosols into the 35 atmosphere, the acceleration of glacier flow by geothermal heating, and outburst flooding from 36 glacial lakes. The rapid flow of the Northeast Greenland Ice Stream has been attributed to 37 elevated geothermal heat fluxes (GHF) due to volcanism6 or the passage of the Iceland plume, 38 perhaps influencing the mass of the Greenland Ice Sheet over glacial-interglacial cycles7. 39 Beneath the West Antarctica Ice Sheet (WAIS), ice flow could be enhanced by elevated GHF 40 from subglacial volcanism8,9 or a mantle plume10. Understanding whether deglaciation enhances 41 continental and/or hotspot magmatism has implications for the retreat of the Greenland and West 42 Antarctic Ice Sheets. 43 Increased mantle melting also enhances the extraction of CO2 from the mantle. If released 44 to the surface, the additional magmatic CO2 can impact the Earth’s climate. During the last 45 deglaciation, subaerial volcanoes are thought to have erupted up to 1000–5000 Gt of additional 46 CO2 (refs. 2,3). Changes in sea-level associated with glacial-interglacial cycles may also enhance 47 CO2 emissions from mid-ocean ridge volcanoes11. However, little work has focused on the 48 enhancement of diffuse subaerial CO2 emissions from hydrothermal systems and dormant 49 volcanoes, despite their large present-day CO2 flux of 170 Mt/yr, representing roughly half of the 50 modern global volcanic CO2 flux12. 51 The link between deglaciation and enhanced mantle melting is most strongly established 52 in Iceland1,5,13, where increases in eruptive volumes coincide with the most rapid stage of the 53 Late Weischelian deglaciation of the Iceland ice sheet from 11–10 ka (BP). While shallower 54 crustal processes may also modulate the magmatic response to deglaciation, the importance of 55 enhanced mantle melting is evidenced by the magnitude of deglacial eruptive rates and the 56 coeval depletion of incompatible trace elements, first modelled by Jull and McKenzie1 (hereafter 57 JM96). 58 By comparison, deglaciation-enhanced melting in continental mantle has not been 59 quantified (with the exception of global ice mass loss scalings2), and observations of enhanced 60 volcanism during deglaciation in intraplate settings are primarily attributed to the triggering of 61 crustal magma chambers14,15. For example, Yellowstone is magmatically active and has 62 experienced rapid deglaciation. During the Pinedale (22–13 ka) and Bull Lake (140–150 ka) 63 glaciations, ice caps covered the Yellowstone caldera and beyond, extending 100 km in radius4. 64 While the Pinedale deglaciation occurred during a period of volcanic quiescence, the Bull Lake 65 deglaciation occurred during the most recent eruptive episode in Yellowstone, the Central 66 Plateau Member rhyolites (170–70 ka). Geological evidence suggests many of these eruptions 67 are syn-glacial16,17. The Central Plateau Member rhyolites were erupted from a large upper 68 crustal sill, maintained by an extensive deeper magmatic system potentially fed by a mantle 69 plume18. During the deglaciation interval there is no evidence that eruptive rates were 70 heightened, nor that the magmatic system was otherwise altered, relative to background 71 rates/trends. However, Yellowstone’s present-day magmatic CO2 flux (~5% of the modern 72 global flux19) is released not by eruptions, but by the crystallization of magmas at depth19. Thus, 73 it remains unclear whether mantle melting rates and associated volatile fluxes are significantly 74 enhanced under thicker continental lithosphere, particularly as glacially induced pressures are 75 attenuated with depth1, and by extension whether the singularly strong response of Iceland is 76 related to the unique juxtaposition of the Icelandic mantle plume and the Mid-Atlantic ridge. 77 In this study, we model deglaciation-enhanced mantle melting in both Iceland and 78 Yellowstone, to gain insight into local eruption rates and the potential for enhanced CO2 fluxes 79 from each system. We use the mantle convection code ASPECT20,21 to simulate changes in 80 pressure and melt production due to glacial unloading for Iceland and Yellowstone (see 81 Methods). The 2-D models are first run to steady-state to resemble present-day “background” 82 behavior (Figure 1) and are then loaded/unloaded using the reconstructed ice load for each 83 system. The models are unloaded by decreasing the ice sheet radius at a constant rate over a 84 prescribed deglaciation interval (1000 years for Iceland, 2000 years for Yellowstone), simulating 85 the retreat of the ice margin. The mantle melt production rate is the rate of melt fraction change 86 integrated spatially. We also calculate trace element concentrations and estimate the flux of CO2 87 segregated from the mantle by melts and the flux of CO2 exsolved to the surface. Finally, we 88 estimate the heat released by the emplacement of additional melts. 89

As an ice mass retreats and unloads the Earth's surface, the underlying mantle rebounds 32 and undergoes a reduction in pressure. If the mantle is above the solidus, this decompression 33 generates additional melting relative to any background rate. Enhanced mantle melting can result in increased volcanic activity 1,5 , which in turn may incite the release of aerosols into the 35 atmosphere, the acceleration of glacier flow by geothermal heating, and outburst flooding from 36 glacial lakes. The rapid flow of the Northeast Greenland Ice Stream has been attributed to 37 elevated geothermal heat fluxes (GHF) due to volcanism 6 or the passage of the Iceland plume, 38 perhaps influencing the mass of the Greenland Ice Sheet over glacial-interglacial cycles 7 . 39 Beneath the West Antarctica Ice Sheet (WAIS), ice flow could be enhanced by elevated GHF 40 from subglacial volcanism 8,9 or a mantle plume 10 . Understanding whether deglaciation enhances 41 continental and/or hotspot magmatism has implications for the retreat of the Greenland and West 42 Antarctic Ice Sheets. 43 Increased mantle melting also enhances the extraction of CO2 from the mantle. If released 44 to the surface, the additional magmatic CO2 can impact the Earth's climate. During the last 45 deglaciation, subaerial volcanoes are thought to have erupted up to 1000-5000 Gt of additional 46 CO2 (refs. 2,3 ). Changes in sea-level associated with glacial-interglacial cycles may also enhance 47 from each system. We use the mantle convection code ASPECT 20,21 to simulate changes in 80 pressure and melt production due to glacial unloading for Iceland and Yellowstone (see 81 Methods). The 2-D models are first run to steady-state to resemble present-day "background" 82 behavior ( Figure 1) and are then loaded/unloaded using the reconstructed ice load for each 83 system. The models are unloaded by decreasing the ice sheet radius at a constant rate over a 84 prescribed deglaciation interval (1000 years for Iceland, 2000 years for Yellowstone), simulating 85 the retreat of the ice margin. The mantle melt production rate is the rate of melt fraction change 86 integrated spatially. We also calculate trace element concentrations and estimate the flux of CO2 87  We estimate the concentration of CO2 in the melt and the flux of CO2 released to the 129 surface. We calculate the partitioning of CO2 into the melt using the retained melt fraction 130 formulation 22 , which can reproduce the magnitude of the observed 5,23 depletion in trace element 131 concentrations due to deglaciation (see Methods and Figure S9a). Our background CO2 fluxes 132 (orange lines in Figure 3d) are within the range inferred from helium fluxes 24 . During the 133 deglaciation, we calculate that for a mantle CO2 content of 300-500 ppm (see Methods), an 134 additional 31-51 Gt of CO2 is released over 1 kyr (dash-dotted black line, Figure 3d), 135 corresponding to a 13-fold increase over the background flux. This additional CO2 is likely not 136 released instantaneously, but is slowed by processes such as melt migration 25 . This value is of 137 the same order of magnitude as prior estimates 25 , which found an extra ~165 Gt CO2 was 138 released over the 11 kyrs following deglaciation for a mantle CO2 content of 285 ppm. 139 Finally, we examine the conditions under which the heat released by the emplacement of 140 the additional melts may reach the surface. The emplacement of our steady-state melt production 141 rate at a depth of 10 km releases 8.7 GW of heat (comparable to the 8 GW estimated in a similar 142 calculation 26 ). This flux may be transferred conductively to the surface over long time scales, and 143 is consistent with borehole measurements from outside the rift zone 26 . During the deglaciation, 144 we estimate the emplacement of the additional melts releases 281 GW at depth, for a total of 145 9×10 21 J over the entire interval. For comparison, the energy required to melt a 100,000 km 3 146 Icelandic ice sheet near its melting point is 30×10 21 J. 147

Deglaciation melting in Yellowstone 156
We next estimate how deglaciation affects mantle melt production rates associated with 157 the Yellowstone plume. Prior to unloading, the background mantle flow field represents a 158 combination of shearing from the westward motion of the North American plate and uplift from 159 the plume (Figure 4a, red arrows). Melts are produced over the depth interval from 90 to 70 km, 160 over a 300-km wide region (orange colors in Figure 4b). The background mantle melt production 161 rate of 0.022 km 3 /yr represents the rate of emplacement of basalts, assuming efficient melt 162 extraction. 163 During the deglaciation, we find that the enhancement of melting beneath Yellowstone is We also test the response of a transient upper mantle thermal anomaly without a plume 174 tail ( Figure S5) and higher melt production rates ( Figure S6). In the case lacking a plume tail, 175 unloading of the transient upper mantle thermal anomaly yields melt production rates that are 176 almost as high (93%) as the case with a plume tail ( Figure S7). In cases with higher melt production rates, greater volumes of additional melt are generated during deglaciation (see 178 Methods). 179 implying deglaciation signatures may not be resolvable ( Figure S9b). We infer that processes 212 governing melt migration through the lithosphere and crust mitigate volcanic activity despite 213 enhanced melting beneath Yellowstone. Understanding the transfer of the mantle melts to the surface is further complicated by the influence of unloading on the shallower magmatic system. 215 Various studies have examined how magma chambers can be triggered by deglaciation 14,15 . 216 Mantle melts may be pumped upwards as the continental lithosphere flexes during 217 deglaciation 31 . We suspect that relative to Iceland, the thickness of the lithosphere beneath 218 Yellowstone and the complexity of its magmatic system make it more difficult to efficiently 219 transport mantle melts to the surface.

Methods 297
We examine deglaciation-enhanced mantle melting beneath Iceland and Yellowstone using 298 the mantle convection code ASPECT 20,21 . The models are sufficiently idealized to facilitate 299 comparison between both settings, yet capture key geodynamic differences and match various 300 observations. We estimate CO2 and heat fluxes to understand the surface impact. 301 The mantle is assumed to behave as a Newtonian visco-elasto-plastic material with a 302 temperature-dependent viscosity. Viscosities are calculated for dry dislocation creep 47  boundary becomes a free surface that deforms in response to applied pressures. After the glacial 320 load is applied, the model is again allowed to stabilize to rule out the influence of the glaciation. 321 The Iceland ice sheet is simulated as a parabola 180 km in radius and 2 km high (as in JM96). 322 However, we assume the load retreats vertically from the margins, while JM96 kept the load 323 radius constant and horizontally thinned the ice sheet thickness. We compare the horizontally 324 thinned load from JM96 (constant radius, decreasing thickness), the vertically retreating load 325 (decreasing radius, constant maximum thickness) shown in Figure 2, and a horizontally and 326 vertically retreating smaller load (following refs. 35,51 ). In vertically retreating simulations, the 327 melt production rate increases through time as the zone of maximum decompression migrates towards the ridge axis where the load is centered ( Figure S4). For the Yellowstone ice cap, we 329 use a radius of 100 km and a height of 1.25 km, yielding a volume of 20,000 km 3 (ref. 4

). 330
Unloading the ice cap horizontally instead of vertically does not influence melt production rates 331 ( Figure S8a). The dimensions of the Yellowstone ice cap correspond to the most recent and well-332 constrained Pinedale deglaciation (15-14 ka), we assume the more relevant penultimate Bull 333 Lake glaciation (~150 ka) retreated similarly. Lengthening the duration of the deglaciation 334 reduces the melt production rate; however, the total volume of melt produced over the entire 335 deglaciation is unchanged and depends solely on the volume of ice lost ( Figure S8b). 336 The rate of melt fraction change depends on the material derivative of the pressure field, 337 which includes both instantaneous (elastic) changes in pressure and isostatic rebound. We also 338 include the dependence of the melt fraction rate on the temperature field due to the effects of 339 latent heat (as in ref. 52 ). We use a dry peridotite solidus 53 , implying our models underestimate 340 melt volumes under hydrated mantle conditions. Melt fractions (Figure 1c   We then calculate the total melt production rate ( Figure S2). Our 2-D model and melt 614 production rate should be equivalent to the values they obtain at the ridge axis, at the x-intercept 615 of their Figure 10a (red stars, Figure S2). JM96 do not specify which solidus they use, but they 616 do limit the region of melting using a 45° triangle truncated at depths 20-112 km. We obtain similar results with the dry solidus of Katz et al. 53 , limited over the same region (compare red 618 lines and stars in Figure S2). 619 We examine step-by-step the different assumptions made in the JM96 model and our 620 primary model, to explain why our estimate is more productive. The main difference arises from 621 the inclusion of enhanced melting from the wings of the melting region in our model (compare 622 red and dashed blue lines in Figure S2 The primary model run presented in the main text of this study (blue line in Figure S5) 676 has a spreading rate of 10 mm/yr and a mantle potential temperature of 1300 °C, excluding the 677 excess plume temperature of 175 °C. In the absence of a plume, these parameters yield a steady-678 state crustal thickness of 7 km. Faster spreading rates of 20 mm/yr (purple line) and warmer 679 Yellowstone 688

S3.1 Effect of mantle temperature 689
The primary model for Yellowstone presented in the main text has a background mantle 690 potential temperature 1320°C, plus an excess plume temperature of 80°C (i.e. 1400°C mantle 691 potential temperature at the plume center). These parameters yield a background mantle melt 692 production rate of 0.022 km 3 /yr (black lines in Figures S6a and 3c), within the range of estimated 693 crustal emplacement rates of basaltic mantle melts 19,29 . We explore the effect of increasing the 694 controversial 66 . We perform a run in which we remove the plume tail, to understand its effect on 710 our model. To do so, we artificially lower the temperature of the mantle at greater depths (>150 711 km), such that only an upper mantle thermal anomaly remains. Removing the plume tail and 712 reducing the influx of material lessens upwelling from the lower mantle. We find the rates of 713 pressure change in the melting region during unloading are similar, implying the 714 presence/absence of the plume tail itself does not affect our results (they are instead controlled 715 primarily by the viscosity of the upper mantle in the melting region, and overlying lithosphere). 716 As the thermal anomaly at the base of the lithosphere was originally set by the plume, this test is 717 not equivalent to explicitly modeling another mechanism (e.g., edge-driven subduction from the 718 slab). 719

S3.3 Effect of loading function 728
Given the thickness of the lithosphere and the great depth of melting beneath 729 Yellowstone, the pressure changes due to unloading are distributed throughout the melting zone. 730 As a result, the style of the ice cap retreat is unimportant ( Figure S8a). Ice caps which are 731 thinned horizontally yield nearly identical rates of melt production as ice caps which retreat 732 vertically, upon radial integration. If the ice cap retreats more slowly, the melt production rate 733 decreases proportionately but the total volume of extra melt produced is unchanged. Our