Assessing Volcanic Controls on Miocene Climate Change

The Miocene period saw substantially warmer Earth surface temperatures than today, particularly during a period of global warming called the Mid Miocene Climatic Optimum (MMCO; ∼17–15 Ma). However, the long‐term drivers of Miocene climate remain poorly understood. By using a new continuous climate‐biogeochemical model (SCION), we can investigate the interaction between volcanism, climate and biogeochemical cycles through the Miocene. We identify high tectonic CO2 degassing rates and further emissions associated with the emplacement of the Columbia River Basalt Group as the primary driver of the background warmth and the MMCO respectively. We also find that enhanced weathering of the basaltic terrane and input of explosive volcanic ash to the oceans are not sufficient to drive the immediate cooling following the MMCO and suggest that another mechanism, perhaps the change in ocean chemistry due to massive evaporite deposition, was responsible.

growth of Antarctic ice sheets (Levy et al., 2016) and pCO 2 decrease (Rae et al., 2021).High rates of organic carbon (C org ) burial during the MMCO are more likely a response to the warming and sea level rise rather than a factor in the subsequent cooling episode (Sosdian et al., 2020).Another explanation is the deposition of calcium sulphate evaporites, a process which increases ocean alkalinity and CO 2 uptake (Shields & Mills, 2021).This is supported by evidence for evaporite deposition across the Mediterranean during the Mid-Miocene (De Leeuw et al., 2010).One unexplored possibility is that volcanic eruptions themselves drove the cooling at the end of the MMCO.It is known basaltic terranes, such as the CRBG, are highly erodible and reactive, thus may increase the global rate of silicate weathering and CO 2 draw-down (Dessert et al., 2001(Dessert et al., , 2003)).These terranes may also supply large amounts of nutrients to the oceans, where in the right conditions they stimulate productivity and sequester carbon via the biological pump (Dessert et al., 2003;Li et al., 2016).Furthermore, explosive volcanic eruptions, and deposition of volcanic ash in the oceans, are known to lead to enhanced drawdown and burial of marine C org (Longman et al., 2019;Longman, Gernon, et al., 2021) , potentially providing a separate sink for CO 2 .
In this study, we use the new SCION (Spatial Continuous Integration) Earth Evolution model (Mills et al., 2021) to investigate changes to the global carbon and nutrient cycles throughout the Miocene, with a particular focus on the MMCO and MMCT.SCION uses a 3-dimensional interpolation scheme to link physical climate model outputs (from the Fast Ocean Atmosphere Model -FOAM; Jacob, 1997) to long-term biogeochemical cycles.The model allows the investigation of changes to surface processes like silicate weathering in 2D over geological time under a realistic steady-state climate dynamically linked to the model carbon cycle (Mills et al., 2021).SCI-ON is a 'forwards' model which makes predictions for climate and global elemental cycles using only tectonic and evolutionary boundary conditions, and relies on a large number of General Circulation Model (GCM) runs performed at discrete points in geological time and for a range of CO 2 levels (Donnadieu et al., 2006;Goddéris et al., 2014).Advances in modeling have been matched by recent advances in compilations of palaeo-proxy data (Rae et al., 2021;Scotese et al., 2021;Westerhold et al., 2020), and so comparisons between the two approaches should be possible at levels of accuracy previously unobtainable, helping us to assess the model predictions.

Model Setup and Scenarios
We use the SCION model version 1.1 for this work, which contains an important update to the first model version (1.0: Mills et al., 2021).In the original model, the Cenozoic climate was represented by GCM runs at 52, 30 and 0 Ma.In this version, we add GCM runs for 15 Ma, which were part of the set used by Goddéris et al. (2014) but had not been converted for the SCION framework due to the whole-Phanerozoic focus of the first SCION paper and some minor issues in scaling the digital elevation model to the same resolution as the climate outputs, which have now been resolved.The model is run forwards in time for the whole Phanerozoic, as in Mills et al. (2021), but we only show the Oligocene-to-present estimates here.We concentrate on the model predictions for atmospheric CO 2 , global average surface temperature and whole-ocean δ 13 C of DIC, and use 1000-member ensembles to test the variation in key parameters as in Mills et al. (2021).An extended description of the model may be found in Text S1 of Supporting Information S1.The full model flux data are included in Figure S1 of Supporting Information S1 and the model code can be downloaded at https://github.com/bjwmills/SCION.

Model Scenarios
Figure 1 shows the evolving continental configuration, including the location of the CRBG, alongside the modified chemical inputs assumed throughout the model runs.The baseline model assumes a decreasing rate of CO 2 input throughout the last 30 million years, based on the declining material subduction rate (Domeier & Torsvik, 2019) and the reduction in the extent of continental arcs and rifts (Brune et al., 2017;Mills et al., 2017).The uncertainty range is described in Mills et al. (2021) and represents the boundaries of the estimates based on arc length alone versus those based on material subduction flux alone.
To estimate CRBG degassing flux (orange in Figure 1d) we use the timeframe of 16.7 to 15.9 Ma, and derive a maximum CO 2 flux estimate in a manner similar to Armstrong McKay et al. (2014).This method compares the size of the CRBG and that of the Siberian traps, probably the most well-studied LIP in terms of its 'cryptic degassing' budget -carbon degassed from the surrounding sediments rather than the basalt itself, which forms the major part of the CO 2 budget (perhaps 90%; Armstrong McKay et al., 2014;Svensen et al., 2009).A conservative estimate places the Siberian traps volume at roughly 1 × 10 6 km 3 (Ivanov et al., 2013), five times the volume of the CRBG (2 × 10 5 km 3 ; Reidel et al., 2013), meaning the contacted surface area of surrounding sediment was around three times that of the CRBG (assuming spherical shape).Thus, an assumed cryptic degassing flux of 7.2 × 10 12 mol C yr −1 from the Siberian traps (Dal Corso et al., 2020;Svensen et al., 2009) would imply a flux of 2.38 × 10 12 mol C yr −1 from the CRBG, for a total maximum possible CO 2 degassing of 2.64 × 10 12 mol yr −1 , assuming the CRBG was emplaced into relatively carbon-rich surroundings (Armstrong McKay et al., 2014).
We include the emplacement of the CRBG as a highly weatherable terrane covering the entire gridbox around latitude 47°N, 60°W (Figure 1b).The gridbox area is 280,000 km 2 , similar to the estimated initial area of the CRBG (Ernst et al., 2021;Mills et al., 2014;Reidel et al., 2013).The CRBG is included in the land map for 15 Ma but not in the present-day simulation, assuming that it's enhanced weathering effect has waned sufficiently (Vance et al., 2009).The weathering effects of the LIP are assumed to begin at 16.7 Ma and are linearly interpolated out of the land map between 15 and 0 Ma.We calculate the weathering effects in a similar manner to previous research (Goddéris et al., 2017;Lefebvre et al., 2013).Total present-day mafic rock area is around 6.8 million square kilometers (Mkm 2 ), contributing ∼35% of the silicate weathering flux (Dessert et al., 2003).The total weatherable land area which may contain silicates (i.e., ignoring carbonates and evaporites) is about 120 Mkm 2 (Hartmann et al., 2014).So, if mafic lithologies represent around ∼5% of weatherable silicates, we determine them to be ∼7 times more weatherable than the homogeneous silicates which comprise the remainder of the model land surface.Thus, we multiply the silicate weathering flux for the CRBG gridcell by a factor 7. Sensitivity tests are detailed in Table S1 of Supporting Information S1.
To consider nutrient supply from the major silicic volcanic episodes in central Europe during the Mid-Miocene (Lukács et al., 2018), we use an ash depletion model adapted from Lee et al. (2018).This approach uses the GeoROC database to estimate the initial composition of erupted volcanic material, before comparing it to measured compositions of altered tephra to estimate the amount of phosphorous loss (Longman, Mills, et al., 2021).This is then used to estimate productivity increases associated with the input of P by including it as an additional P flux to the model ocean.There are two clear pulses of volcanic input at 14.36 and 14.88 Ma, associated with the eruption of the Harsany and Demjen tephras (Lukács et al., 2018).To estimate the volume of these eruptions and the amount of P supplied, we use Monte Carlo simulations considering errors on a range of variables.For the two relevant tephras, we estimate their volume as a proportion of the total 4,000 km 3 erupted across the Mid Miocene, using their relative extent (compared to the total extent of all tephras) as observed in outcrops (Lukács et al., 2018).For Monte Carlo simulations, we estimate error on these values to be 20%.We then convert these values to P supply using ash density, P content and P depletion factors (and associated errors) from Longman, Mills, et al. (2021).Finally we assume 90% of the erupted tephra was deposited in the ocean, based on paleogeographic reconstructions (Lukács et al., 2018).These formations were deposited in what was the Pannonian Sea, a back-arc basin which persisted in central Europe until the Pliocene Epoch (Balázs et al., 2016).All variables and their associated errors used in the Monte Carlo modeling may be found in Table S2 of Supporting Information S1.For simplicity we model these eruptions as a single Gaussian input curve which begins at the timing of the Harsany tephra and peaks at the timing of the Demjen tephra, as the model is not designed to represent timescales of less than 100 kyrs.

Results and Discussion
The SCION model baseline run (Figures 2a-2c; no LIP and no tephra) simulates the long-term temperature decline of the Miocene well, but overestimates CO 2 before ∼10 Ma (Figure 2b).This is likely due to the low CO 2 climate sensitivity of the FOAM climate model (Donnadieu et al., 2006;Goddéris et al., 2014), meaning relatively high CO 2 levels are required to raise the surface temperature to a point of carbon cycle balance (Mills et al., 2021).Climate and carbon cycle balance at 15 Ma has previously been simulated in the GEOCLIM model using the same climate model runs that we use here (Lefebvre et al., 2013).The previous investigation used a present-day CO 2 degassing rate rather than an elevated rate of degassing, resulting in lower steady state CO 2 concentrations than the current model -∼500 ppm.Here, the increase over preindustrial levels was due partly to the positioning of weatherable basaltic terranes in arid areas, not considered in the present work.Returning to the SCION predictions, a clear mismatch between the baseline output and proxy data is visible during the MMCO (∼17-14 Ma) and the following MMCT (13.9 -12 Ma) where the model does not capture temperature, CO 2 and the carbon isotope composition change associated with climatic variability (Figure 2).
The addition of degassing from the emplacement of the CRBG improves this mismatch considerably across the MMCO, especially regarding global temperature change (Figures 2d-2f).Previous modeling used a smaller CO 2 input, even when a cryptic element was considered, insufficient to force the MMCO in its entirety (Armstrong McKay et al., 2014).The difference here is that our calculations are based on new chronological constraints indicating the majority of CRBG emplacement in as little as 750 kyr (Kasbohm & Schoene, 2018), and our method of comparison to the Siberian traps considers the difference in contacted sediment surface area directly.The injection of CO 2 associated with this event is sufficient to drive the MMCO inception, and constraining the majority of CO 2 degassing to a short period after 16.7 Ma places the CRBG emplacement at the same time as the onset of the MMCO (Kasbohm & Schoene, 2018).Such a finding supports previous assertions of a causative link between Mid Miocene volcanism and climate change (Hodell & Woodruff, 1994;Kasbohm & Schoene, 2018;Kürschner et al., 2008;Steinthorsdottir et al., 2021).
The model estimate of global temperature increase at the MMCO onset is ∼1°C, lower than some proxy reconstructions (e.g., 2.3°C rise; Scotese et al., 2021).However, this is well within proxy reconstruction variability and is also somewhat affected by a low climate sensitivity in FOAM compared to more complex modern climate models (e.g., Zhu et al., 2019).It is possible that other greenhouse gases were also released as part of the CRBG emplacement, with evidence for hydrothermal alteration of organic-rich sediments and associated CH 4 release (Bindeman et al., 2020), potentially enhancing the greenhouse effect.Comparison of the model output with reconstructions of CO 2 further supports an assertion of CRBG degassing as a primary driver of Mid Miocene climatic change (Figure 2), with simulated CO 2 of ∼1,200 ppm at the MMCO peak in line with previous reconstructions (Rae et al., 2021), but above a number of previous estimates (e.g., Steinthorsdottir et al., 2020Steinthorsdottir et al., , 2021)).
The spatial pattern of modeled surface air temperatures at 15 Ma (Figure 3a) shows general agreement with a data set of MMCO proxy reconstructions of temperature (Goldner et al., 2014).The largest offsets between model and proxy data are in polar regions, where the model simulates anomalously low temperatures (Figure 3a).This is a common issue when modeling warm periods of the Earth history such as the Eocene and the Pliocene (Haywood  , 2020;Huber & Caballero, 2011), and it remains a challenge to successfully model the Miocene polar regions with fully coupled ocean-atmosphere models (Burls et al., 2021).
We assume a value of zero for the δ 13 C composition of the CRBG-degassed CO 2 .While mantle-derived carbon would be expected to have a negative fractionation compared to the PDB standard, most of the degassed carbon is assumed to have come from sediments with varied fractionations.The SCION model has an average sediment δ 13 C value of zero, used here in the absence of more detailed constraints.It is possible that Miocene volcanic carbon contained more positive δ 13 C than usual, although there is little experimental evidence of such δ 13 C values in mantle carbon (Deines, 2002).Under the current scenario, the model produces a small positive δ 13 C excursion during the CRBG degassing due to elevated temperature, weathering, and productivity.This suggests the hypothesis of Sosdian et al. (2020) may be correct; high atmospheric CO 2 across the MMCO led to a positive shift in the fractionation factor associated with photosynthesis, offsetting any negative mantle δ 13 C emissions (Armstrong McKay et al., 2014;Sosdian et al., 2020).This additional factor is not considered in the model, but the error window on δ 13 C shown encompasses likely changes due to photosynthetic fractionation (Mills et al., 2021).
By also modeling enhanced silicate weathering from the CRBG, we consider the impact of LIP emplacement holistically across the Mid Miocene (Figures 2g-2i).Basaltic terranes such as the CRBG are highly erodible and the products which enter the oceans contain abundant silicate cations and are nutrient-rich (Dessert et al., 2001(Dessert et al., , 2003;;Gaillardet et al., 1999;Hartmann et al., 2009;Li et al., 2016).As such, it has been suggested that as well as increasing silicate weathering rates, these nutrients may enhance productivity, and raise the levels of marine organic carbon (C org ) burial, a process which removes CO 2 from the ocean-atmosphere system and leads to cooling on long timescales (Berner et al., 1983;Li et al., 2016;Li & Elderfield, 2013).
The inclusion of CRBG weathering in the SCION model has little impact on the model results (Figures 2g-2i).This is because although the terrain is assumed to contribute seven times more silicate weathering (and P delivery) than the rest of the model homogeneous lithology, it still only represents a small fraction of the land surface and is emplaced in an area of relatively low temperature in the high latitudes.The position of the CRBG in the model is in an area where temperature is well-predicted (Figure 3a), and runoff rates are reasonably high (Figure 3b).Accordingly, the weathering rate of the CRBG is estimated to be more than two times higher than any local gridcell.Despite this clear contribution to high-latitude weathering, its overall weathering contribution is similar to a single equatorial gridcell.
It has been proposed that MMCO length was extended by a delayed silicate weathering feedback (Babila & Foster, 2021), but our model results suggest the warm interval length is consistent with the single injection of mantle CO 2 , with no clear requirement for a delay in the onset of CRBG weathering to reproduce proxy data (Figure 2).Additionally, the scale of sequestration associated with silicate weathering of the CRBG does not appear to be sufficient to drive the transition to cooler climate after ∼13.9 Ma, the end of the MMCO (Steinthorsdottir et al., 2020), with proxy evidence suggesting a much more rapid shift to cooler temperatures than the SCION reconstruction (Figure 2g).
We now consider other potential controls on Miocene cooling.In central Europe, extremely large tephra deposits (∼4,000 km 3 tephra per eruption) are preserved (Lukács et al., 2018).Since these primarily erupted into the Paratethys Ocean, and as the alteration of tephra is known to release nutrients (Browning et al., 2015;Frogner et al., 2001;Jones & Gislason, 2008;Longman et al., 2019), we consider these as a potential source of bioavailable nutrients to the ocean system.Dates for these eruptions indicate two pulses, the first at 18.2-16.8Ma, and the second between 14.9 and 14.4 Ma (Lukács et al., 2018), with evidence of ashfall from southern Germany (Arp et al., 2021) to western Russia (Danišík et al., 2021) indicative of the scale of the events.Using a separate modeling approach (see Methods), we are able to estimate the scale of P supply from these eruptions (Figure 1).The input of additional P to the ocean plays very little role in global-scale climate change in our model (Figures 2j-2l).This contrasts with a separate study, which linked the Ordovician glaciations to tephra-derived P supply (Longman, Mills, et al., 2021), indicating the volcanic episodes observed in central Europe across the mid-Miocene were not large enough to play a major role in global climate cycles -the P supply we calculate here is more than an order of magnitude smaller than the Ordovician events.
Another mechanism which may have enhanced the cooling feedback following the CRBG emplacement is the relatively high sea levels during the MMCO, a result of Antarctic ice shelf loss (Gasson et al., 2016;Lear et al., 2010;Steinthorsdottir et al., 2020).Higher sea levels would have flooded large areas of shallow continents, and resulted in larger areas of continental shelves (Sosdian et al., 2020), locations known to be centers of C org burial (Bianchi et al., 2018;Burdige, 2007).It is possible high CO 2 atmospheres as a result of volcanism during the MMCO provided a source of carbon to the upper ocean, with nutrients additionally supplied by weathering of the CRBG (Sosdian et al., 2020).Easily available nutrients, high atmospheric CO 2 , and high temperatures combined with the large extent of continental shelves, likely resulted in the deposition of large areas of C org -rich sediments in the circum-Pacific area (e.g., the Monterey formation; Flower & Kennett, 1993;Pisciotto & Garrison, 1981).There is evidence to suggest the high C org contents are not related to exceptionally high productivity, rather near-perfect conditions for preservation (Föllmi et al., 2005;Laurent et al., 2015).Whatever the exact mechanism, the timing of this event (∼16-14 Ma; Sosdian et al., 2020) indicates it cannot have played a role in the MMCT, but was rather a feature of the warm temperatures in the aftermath of the CRBG emplacement.Our modeling supports this assertion, with global C org burial rates increasing by 10% during the MMCO.
Another proposed cooling mechanism at the end of the MMCO is the closure of the Panama gateway (Montes et al., 2015), which led to the formation of the Atlantic Meridional Overturning Circulation (AMOC), which in the modern ocean acts as the single largest carbon sink in the Northern Hemisphere (Gruber et al., 2002).Research has suggested this occurred between 15 and 12 Ma, a period which could have implications for the end of the MMCO and the MMCT (Bacon et al., 2015;Montes et al., 2015).However, most evidence suggests ongoing throughflow between the Pacific and the Atlantic until the Pliocene (Karas et al., 2017;Lessios, 2015;O'Dea et al., 2016;Panitz et al., 2018).
The timing of the end of the MMCO also correlates well with evidence for intensification of physical weathering in the Himalaya Mountains linked to the onset of the winter monsoon (Ali et al., 2021;Bretschneider et al., 2021;Lee et al., 2020).This supports studies which proposed that enhanced weathering of the Himalayas led to Miocene cooling but were limited by temporal uncertainties (Allen & Armstrong, 2012;Tada et al., 2016).Such a process would have enhanced CO 2 drawdown, and may have led to cooling (Allen & Armstrong, 2012), but it is unlikely this process was sizable enough to have controlled the end MMCO cooling (Clift & Jonell, 2021).Earlier work using these FOAM simulations has shown that the South Asian Monsoon amplified silicate weathering rates at the present-day, relative to 15 Ma (Lefebvre et al., 2013), and this feature persists in the SCION model.More detailed testing of monsoonal effects on climate on the Myr scale would require a larger set of continental configuration and climate model runs.
Another potential cause for the rapid cooling observed at the MMCT is widespread evaporite deposition (Shields & Mills, 2021).Deposition of thick calcium sulphate evaporites depletes ocean calcium, impacting calcium carbonate precipitation and increasing marine alkalinity.This alkalinity increase leads to drawdown of CO 2 into the ocean on a million-year scale, before the silicate weathering feedback balances the alkalinity and calcium budget and restores equilibrium (Shields & Mills, 2021).It has been demonstrated that around 10 19 mol of gypsum deposition would roughly halve oceanic Ca concentration, potentially driving ∼1.5 degrees of cooling over 1 Myr (Shields & Mills, 2021).The mid-Miocene did indeed see widespread evaporite deposition, culminating in the Badenian Salinity Crisis at around 13.8 Ma and the deposition of evaporite facies across much of eastern Europe (De Leeuw et al., 2010;Peryt, 2006).The total weight of sulphate evaporites deposited is uncertain, but the estimated budget for the whole Miocene (including the later Messinian Salinity Crisis) is around  1.5 × 10 19 mol (Hay et al., 2006), so is compatible with a large Langhian-Serravallian deposition contributing to the MMCT.The impact of evaporite disposition on climate cannot be explicitly modeled in SCION because it employs steady-state carbonate chemistry, but it is feasible that the disconnection between model and proxy data 15-13 Ma is linked to this process (Figures 2j and 2k).

Conclusions
Using a new coupled biogeochemical-climate model, we investigate the driving forces behind the climate of the Miocene.We demonstrate the emplacement of the Columbia River Basalt Group, and associated CO 2 degassing led to the warming of the Mid Miocene Climatic Optimum (MMCO; 17-15 Ma).The modeled response of the Earth system to the injection of magmatic CO 2 is very similar in magnitude and extent of the warming and CO 2 increase observed in proxy records from this time.The cooling following the MMCO cannot be easily explained by either weathering of the basalt group, or a greater phosphorus supply associated with large-scale volcanism.The abundant organic carbon deposition observed in Pacific facies across the MMCO (the Monterey event) is also unlikely to have driven to the cooling, rather it appears this was a feature of high MMCO temperatures leading to greater productivity.We suggest the most likely cause of the rapid cooling period is the emplacement of giant evaporites in central-eastern Europe, which depleted the ocean calcium reservoir and raised ocean alkalinity.

Figure 1 .
Figure 1.Summary of topographical change for the three GCM model runs used in this study: 30 Ma (a), 15 Ma (b) and 0 Ma (c).Green and brown grid cells indicate topography below and above 1 km respectively.Panel B shows the location of the Columbia River Basalt Group (CRBG).D displays the long-term CO 2 degassing rate (orange line) and the two proposed forcings tested by this work, the input of excess CO 2 from the emplacement of the CRBG, and phosphorous input to the oceans from the deposition of tephra in central Europe.

Figure 2 .
Figure 2. Model outputs of global average surface temperature (GAST), atmospheric CO 2 and carbonate δ 13 C, comparing various SCION runs through the Miocene.(a-c) display the baseline SCION run, with no additional forcings.(d-f) show the baseline with additional LIP CO 2 degassing.(g-i) show baseline with additional LIP CO 2 degassing and basaltic terrane weathering.(j-l) show baseline with additional LIP CO 2 degassing, basaltic terrane weathering and phosphorous supply from tephra input.The surface temperature approximation is calculated from the oxygen isotope data set of Westerhold et al. (2020), under the temperature conversion of Hansen et al. (2013).CO 2 proxy data is from boron isotope measurements (indicative of pH) from Rae et al. (2021).Carbon isotope data is from Westerhold et al. (2020).

Figure 3 .
Figure 3. Gridded map outputs for the complete SCION model scenario shown in Figures 2j-2l at 15 Ma.Displayed here are surface air temperature (a), continental runoff (b) and silicate weathering rates (c).Panel A includes a comparison of model surface air temperature at 15 Ma with proxy reconstructions of mean annual temperature (filled circles; Goldner et al., 2014), with map colors the same as the individual proxy reconstructions.All panels show the location of the Columbia River basalt using a dashed line.