The Snowball Stratosphere
R. J. Graham Tiffany Shaw Dorian Abbot





















R.J. Graham1, Tiffany A. Shaw2, Dorian S. Abbot2
1Department of Physics,
  The University of Oxford, Oxford, UK
2Department of Geophysical Sciences, University of
  Chicago, Chicago, IL, USA


  * The simulated snowball stratosphere with or without ozone displays weaker circulation than the modern case.

  * The snowball stratosphere with ozone has weaker wave forcing and a stronger jet than the modern and has no sudden stratospheric warmings.

  * Changes to stratospheric circulation during the snowball do not impact pCO_2 estimated from proxy data unless pCO_2 was much greater than pO_2.




















According to the Snowball Earth hypothesis, Earth has experienced periods of low-latitude glaciation in its deep past. Prior studies have used general circulation models (GCMs) to examine the effects such an extreme climate state might have on the structure and dynamics of Earth's troposphere, but the behavior of the stratosphere has not been studied in detail. Understanding the snowball stratosphere is important for developing an accurate account of the Earth's radiative and chemical properties during these episodes. Here we conduct the first analysis of the stratospheric circulation of the Snowball Earth using ECHAM6 general circulation model simulations. In order to understand the factors contributing to the stratospheric circulation, we extend the Statistical Transformed Eulerian Mean framework. We find that the stratosphere during a snowball with prescribed modern ozone levels exhibits a weaker meridional overturning circulation, reduced wave activity, stronger zonal jets, and is extremely cold relative to modern conditions. Notably, the snowball stratosphere displays no sudden stratospheric warmings. Without ozone, the stratosphere displays slightly weaker circulation, a complete lack of polar vortex, and even colder temperatures. We also explicitly quantify for the first time the cross-tropopause mass exchange rate and stratospheric mixing efficiency during the snowball and show that our values do not change the constraints on CO_2 inferred from geochemical proxies during the Marinoan glaciation (∼635 Ma), unless the O_2 concentration during the snowball was orders of magnitude less than the CO_2 concentration.





























§ INTRODUCTION





The snowball Earth hypothesis proposes that during the Neoproterozoic glaciations (∼710 and ∼635 Ma) ice covered most or all of Earth's oceans and continents from the poles to very low latitudes, possibly to the equator <cit.>. During this time, silicate weathering would be greatly reduced, allowing CO_2 from volcanic outgassing to build up to the extremely high levels needed to melt the tropical ice and reverse the global glaciation <cit.>. Although there are alternative models that attempt to explain observations of low-latitude Neoproterozoic glaciations without invoking a full snowball <cit.>, this study focuses on the asymptotic case of a full, global glaciation.

How life survived the snowball events and how deglaciation occurred are important questions that haven't been fully resolved <cit.>. Atmospheric circulation during a snowball influences both of these issues through heat transport; the lofting of dust, sometimes referred to as cryoconite, which affects surface albedo and the formation of holes in the ice; and clouds, which can generate greenhouse warming. So far most investigations have focused on the snowball troposphere. For example, previous work has shown that 1.) the tropopause in a snowball is lower <cit.>; 2.) the Hadley circulation is more intense <cit.>, thermally indirect in the annual mean <cit.>, and strongly influenced by vertical diffusion of zonal momentum <cit.>; and 3.) baroclinic eddy dry static energy transport is proportionally more important than in its modern counterpart <cit.>. The only investigation of the snowball's stratosphere that we have encountered is that of <cit.>, which examined the radiative effect of ozone on the snowball atmosphere's temperatures and found that reduced ozone leads to a cooler surface. The stratospheric circulation during a snowball, however, has not been investigated in detail. 


The dynamics of the snowball stratosphere and stratosphere-troposphere exchange have important consequences for the interpretation of a widely cited piece of geochemical evidence for the snowball. The atmospheric Δ^17O (= δ^17O-0.52×δ^18O), a measure of the mass-independent fractionation of oxygen isotopes 17 and 18, reflects changes in the concentration of CO_2 via a stratospheric photochemical reaction network linking O_3, O_2, and CO_2 <cit.>. With other quantities held constant, greater CO_2 content leads to more negative Δ^17O. There are large negative anomalies in Δ^17O fractionation recorded in sulfates that were produced during and immediately after the second proposed Neoproterozoic glaciation (the Marinoan glaciation, 635 Ma), which have been interpreted as evidence for a very high CO_2 partial pressure from volcanism in the absence of silicate weathering <cit.>. The inferred quantity of CO_2 depends on the circulation within the stratosphere and the stratosphere-troposphere exchange because the signal observed in rocks derives from fractionation in the upper stratosphere that dilutes as it propagates into the troposphere. 

The modern stratospheric circulation is predominantly wave-driven and strongly impacts stratospheric temperatures and ozone abundance <cit.>. In particular, waves generated in the troposphere interact with the winter hemisphere's stratospheric polar vortex and converge in the stratosphere driving the Brewer-Dobson circulation (poleward motion in midlatitudes and descent over the pole) via a Coriolis torque. When the wave-driving is particularly strong it can reverse the sign of the stratospheric polar vortex from eastward to westward, an event called a sudden stratospheric warming <cit.>. Future projections suggest an acceleration of the stratospheric circulation under warming induced by greenhouse gas forcing <cit.> due to enhanced wave breaking in the lower stratosphere. How the different factors affecting the stratospheric circulation, e.g. wave driving, are impacted during cold conditions like the snowball has not been investigated.

Here, we quantify the stratospheric circulation during a snowball using a general circulation model (GCM). With modern ozone, the simulated snowball stratosphere displays a weaker circulation, weaker wave activity, stronger jets, an absence of sudden stratospheric warmings, and colder temperatures compared to conditions in our modern simulation. In the absence of ozone, the stratospheric temperatures are colder, circulation is weaker, wave-driving is weaker, and stratospheric zonal jets are non-existent. We also find that these changes to the circulation probably do not change geochemical proxy-based estimates of CO_2 during the snowball that were calculated with the implicit assumption that the snowball stratosphere behaved like the modern, unless the snowball's O_2 concentration was much lower than its CO_2 concentration. In the remainder of the paper, we describe the simulations and methods we used to diagnose and analyze the properties of the snowball stratosphere (Section <ref>), summarize our results (Section <ref>), and speculate about implications of this work while also describing future directions (Section <ref>).












§ METHODS




 §.§ Model Configuration

We ran the ECHAM6 <cit.> GCM with boundary conditions set to represent the pre-industrial modern Earth and several variants of an idealized hard snowball Earth. We ran each simulation for a five year spin-up period, checked for equilibrium to within 1 W m^-2 at the top of the atmosphere, and ran the simulations for an additional five years to generate the data used in this study. All simulations have a 360-day year, zero eccentricity, and obliquity equal to 23.4^∘, and they use a modern continental configuration with a 50m deep slab ocean. We do not include ocean energy transport because it does not affect the snowball circulation, and our simulated modern stratospheric circulation agrees well with ERA reanalysis without it. The slab ocean is covered in sea ice in the snowball simulations. We use a prescribed seasonally-varying modern ozone profile, except in the simulation that lacks ozone. We consider a simulation without ozone because O_2 levels are poorly constrained during the Neoproterozoic and may have been one or two orders of magnitude lower than today <cit.>. The snowball simulations have CO_2 concentrations of 280 ppmv and 10^4 ppmv in order to test the effect of CO_2 build up over time during the snowball due to weakened silicate weathering. We also cover the land surface with glaciers in the snowball simulations and use a constant albedo of 0.60 for both glaciers and sea ice. The glaciers do not change the topography of the land. We ran one of the snowball simulations with flattened topography to test the influence of topography on circulation and eddy activity. A summary of the different simulations can be found in Table <ref>.



 §.§ Stratospheric circulation


The stratospheric circulation is typically quantified using the transformed Eulerian mean (TEM) circulation, which accounts for the impact of wave driving on the circulation <cit.>. However, the TEM is only approximate, as it assumes small amplitude eddies. A more accurate representation of the stratospheric circulation is the dry isentropic circulation, which does not require any eddy amplitude approximation (see e.g. <cit.>). <cit.> introduced the Statistical Transformed Eulerian Mean (STEM) as a simple and accurate approximation of the dry isentropic circulation that allows one to separate the mean and eddy components. The main assumption for the STEM framework is that the eddy variance is finite and Gaussian. In the limit of small eddy variance STEM converges to the TEM formulation.

Like the TEM circulation, the STEM circulation can be decomposed into mean and eddy components:


    Ψ_STEM(θ,ϕ)   = Ψ_eddy(θ,ϕ)+Ψ_mean(θ,ϕ)
       =∫_-∞^θ∫_0^∞2π a cosϕ/gv'θ'(θ-θ)/√(2πθ'^2^3)exp[-(θ-θ)^2/2θ'^2]dpdθ
       +∫_-∞^θ∫_0^∞2π a cosϕ/gv/√(2πθ'^2)exp[-(θ-θ)^2/2θ'^2]dpdθ

where ϕ is latitude, v is the meridional wind velocity, θ is potential temperature, g is surface gravity, p is pressure, overbars represent zonal and monthly averages, and primes represent deviations from those averages. Ψ_mean approximates the Eulerian mean circulation, including the three-cell structure in isentropic coordinates (Hadley, Ferrel, and polar cells, see Figures <ref> and <ref>, although the polar cells are unresolved at that contour interval). Ψ_eddy approximates the circulation driven by eddy heat fluxes, which is analogous to a Stokes drift in isentropic coordinates.

As mentioned previously, diffusion has been found to dominate the tropospheric circulation in simulations of the snowball <cit.>. The mean circulation Ψ_mean depends on the mean meridional flow v. Here we extend the STEM formulation to specifically account for advection of mean meridional momentum, eddy momentum flux convergence, and unresolved forces (which includes diffusion). Following <cit.> we can use the zonal-momentum budget in the free atmosphere to decompose the mean meridional flow into components (v = v_adv + v_EMFC + v_F) related to mean advection (v_adv), eddy momentum flux convergence (EMFC, v_EMFC), and residual (v_F):


    ∂u/∂ t = 0    = fv+ζv-ω∂u/∂ p - 1/acos^2ϕ∂(u'v'cos^2ϕ)/∂ϕ-∂u'ω'/∂ p+F^u_D
    v_adv   = 1/f(-ζv+ω∂u/∂ p)
    v_EMFC   =1/f(1/acos^2ϕ∂(u'v'cos^2ϕ)/∂ϕ+∂u'ω'/∂ p)
    v_F   = -F_D^u/f

where f = 2Ωsinϕ, ζ is relative vorticity, ω is vertical pressure velocity, p is pressure, u is the zonal wind, and F_D^u is the inferred residual component of the zonal momentum budget, which we assume to be dominated by diffusion. With this decomposition, the STEM formulation for the mean circulation can be written as the sum of advection, EMFC, and residual contributions, Ψ_mean = Ψ_adv + Ψ_EMFC + Ψ_F, where:



    Ψ_adv   = ∫_-∞^θ∫_0^∞2π a cosϕ/gv_adv/√(2πθ'^2)exp[-(θ-θ)^2/2θ'^2]dpdθ
    Ψ_EMFC   = ∫_-∞^θ∫_0^∞2π a cosϕ/gv_EMFC/√(2πθ'^2)exp[-(θ-θ)^2/2θ'^2]dpdθ
    Ψ_F   = ∫_-∞^θ∫_0^∞2π a cosϕ/gv_F/√(2πθ'^2)exp[-(θ-θ)^2/2θ'^2]dpdθ



This decomposition therefore allows us to quantify the importance of diffusion for the stratospheric circulation. Note that because the meridional flow does not take into account the momentum balance at the surface these circulations do not close at the surface.

Finally, to better understand changes to the strengths of the “shallow branches” and “deep branch” of the eddy-driven circulation (e.g. <cit.>), following <cit.> we also decompose v'θ ' into contributions from planetary-scale (k=1-3) and synoptic-scale (k = 4-30) eddies and calculate the circulation for each. With these decompositions of the eddy-driven and mean components of Ψ_STEM, we can understand the drivers of the stratospheric circulation.









 §.§ The Eliassen-Palm flux


In order to quantify the importance of wave driving in the stratosphere we follow <cit.> and calculate the Eliassen-Palm (EP) flux convergence in the vicinity of the winter stratospheric jet:


    ∫_p_1^p_2∫_ϕ_1^ϕ_2   ∇·F
       =a^-1cosϕ∫_p_1^p_2 F_(ϕ) dp|_ϕ_2 - a^-1cosϕ∫_p_1^p_2 F_(ϕ) dp|_ϕ_1
        + ∫_ϕ_1^ϕ_2cosϕ F_(p) dϕ|_p_2 - ∫_ϕ_1^ϕ_2cosϕ F_(p) dϕ|_p_1

where F_(ϕ) = acosϕ[-u'v'+(∂u/∂ p)v'θ '/(∂θ/∂ p)], F_(p) = acosϕ[(f+ζ)v'θ '/(∂θ/∂ p)-u'ω'], ϕ_1 is 40.0^∘, ϕ_2 is 90.0^∘, p_1 is the average pressure at the tropopause as defined by the 380 K isentrope in the modern case and the 285 K isentrope in the snowball cases, and p_2 is 1.0 hPa. With the integrated convergence of the EP flux, we can compare the propagation of wave activity between the troposphere and stratosphere in our modern and snowball simulations. This provides information useful for explaining differences in the structure and behavior of the zonal-mean flow of the stratosphere in our simulations.



 §.§ Tropospheric Δ^17O

Δ^17O_O_2 is a measure of the degree of non-mass-dependent depletion of heavy oxygen (^17O and ^18O) in O_2. This depletion occurs because the photochemical formation of ozone from O_2 in the upper stratosphere tends to produce O_3 molecules with mass independent enrichment of ^17O and ^18O <cit.>. This ozone fractionation signal is transferred to CO_2 through isotopic exchange reactions that move the excited species O(^1D), which is released in UV photolysis of O_3 and therefore enriched in heavy oxygen, to CO_2 <cit.>. Transfer of heavy oxygen atoms from O_2 to O_3 to CO_2 leads to a depletion of heavy oxygen in O_2 (a negative Δ^17O_O_2), and a larger reservoir of CO_2 allows for the storage of more heavy oxygens and correspondingly more negative Δ^17O_O_2 <cit.>. 

The interpretation of anomalously negative Δ^17O_O_2 in ancient minerals as evidence for high CO_2 during and immediately after the Neoproterozoic snowball requires that O_2 bearing this signal propagated down from the stratosphere into the troposphere without entirely losing its mass-independent fractionation. In the rest of the study, we will use Δ^17O to refer to the fractionation of tropospheric O_2 specifically unless otherwise noted. In <cit.>, the heavy oxygen fractionation of tropospheric oxygen in steady state is given as:


    Δ^17O   = δ^17O-0.52×δ^18O
       = -Φ(ρ)γΘτ/1+ρ+γΘτ


where ρ is the ratio of pO_2/pCO_2 (pO_2 and pCO_2 are the concentrations of the respective gases), Φ(ρ) (≡0.519×(64+146×(ρ/1.23))/(1+ρ/1.23)-7.1738) is the difference of Δ^17O between steady state O_2 and CO_2 in the photochemical reaction, τ is the residence time of oxygen (≡pO_2/F_op, and F_op is the rate of oxygen production), γ is the cross-tropopause mass exchange rate (defined later in section <ref>), and Θ is the stratospheric mixing efficiency (defined later in section <ref>). 

Qualitatively, equation <ref> means that the degree to which the fractionation signal can propagate into the troposphere from the stratosphere is determined by 1.) the magnitude of fractionation in the steady-state reaction network in the upper stratosphere (a function of ρ); 2.) how well the depleted O_2 can mix from the upper stratosphere to the lower stratosphere (captured by Θ); 3.) the rate of cross-tropopause mass exchange (γ) to carry the fractionation signal from the lower stratosphere into the troposphere; and 4.) the residence time of oxygen (τ), since O_2's residence time is determined by its rate of photosynthetic production, and photosynthesis produces oxygen that lacks mass-independent fractionation, so higher rates of O_2 production (e.g. smaller τ) will tend to wash out the stratospheric fractionation signal in the troposphere. That is why ρ, γ, Θ, and τ all appear in the numerator. ρ appears in the denominator because very large ρ (e.g. pO_2 ≫ pCO_2) limits the fraction of heavy oxygens from O_2 that can be sequestered in CO_2 molecules. γ, Θ, and τ appear in the denominator because the magnitude of the tropospheric Δ^17O anomaly cannot exceed the value reached in the upper stratosphere, so there must be a limit past which more efficient transport of altered oxygen from the upper stratosphere to the troposphere (e.g. increases in γ or Θ) or reduced photosynthetic flux of unaltered oxygen (e.g. increases in τ) stop affecting the troposphere's value. We focus on estimating γ and Θ in this study.




  §.§.§ Stratosphere-troposphere exchange


The cross-tropopause mass exchange rate γ can be represented as the ratio between the flux of mass across the tropopause and the total mass of the atmosphere, giving a quantity with units of yr^-1. The γ used in the snowball calculations in <cit.> to infer high CO_2 during the snowball is derived from data in <cit.>, which describes the modern atmosphere. We directly calculate the upward mass flux across the tropopause in our simulations by taking the difference between the sum of all local maxima and the sum of all local minima of the STEM circulation along an isentrope that roughly coincides with the highest portion of the tropopause in potential temperature coordinates (380 K in the modern simulation, and 285 K in the snowballs):


    γ   = ∑Ψ_STEM(θ_tropo)_max-∑Ψ_STEM(θ_tropo)_min/M_atm
       =F_STE/M_atm


where F_STE is the mass flux across the tropopause, θ_tropo is potential temperature chosen to represent the tropopause, and M_atm is the mass of the atmosphere as determined by the pressure it exerts at Earth's surface. The globally averaged surface pressure in each of the simulations is 9.855× 10^4 Pa. Here, the surface pressure refers to the actual pressure on the ground, including the effects of topography, and not the mean sea level pressure. M_atm is calculated by dividing the surface pressure by the gravity g. We calculate the mass flux across the isentrope that coincides with the highest portion of the tropopause in our simulations (instead of directly across the tropopause) because those potential temperatures approximately coincide with the bottom of the atmospheric region in which isentropes fall entirely within the stratosphere, e.g. the “stratospheric overworld” <cit.>. There is a large amount of adiabatic mass exchange with the troposphere along the isentropes in the “lowermost stratosphere” that cross the tropopause. Since that portion of the stratosphere is more tightly coupled to the troposphere by mass exchange than it is to the “overworld” <cit.>, it is better considered a part of the troposphere for the purpose of calculating Δ^17O.



  §.§.§ Stratospheric mixing efficiency

The stratospheric mixing efficiency, denoted Θ in this study, is defined in <cit.> as the “fraction of the steady-state O_2 in the O_2–CO_2– O_3 reaction network within the total O_2 in the stratosphere.” In <cit.>, Θ is inferred from equation (<ref>) using modern values of Δ^17O (-0.34), pCO_2 (375 ppmv), pO_2 (209,500 ppmv), and τ (1244 years). To examine whether this parameter may have been different during a snowball, changing estimates of pCO_2 made using the Δ^17O values presented in <cit.>, we estimate Θ using data from our simulations. 

ECHAM6 does not contain interactive chemistry, and there is no standard formula to directly calculate Θ. To estimate Θ, we calculate the fraction of mass at a given altitude that has passed through the region of the stratosphere where rates of photochemistry would peak, but has not yet passed through the troposphere, and we take Θ to be directly proportional to this quantity at the tropopause. We call this quantity ξ and treat it as the concentration of a passive tracer, and we use “A” to denote the scaling factor that relates ξ(tropopause) to Θ, such that Θ = A×ξ(tropopause). Assuming γ×Θ from our modern simulation to be equal to the γ×Θ given in <cit.> allows us to calculate the proportionality constant A, since in that case A = γ_Cao×Θ_Cao/γ×ξ. We can then use this value of A and the ξ values calculated from our snowball simulations to estimate how Θ changes under snowball conditions.

We assume that the altitude where the prescribed O_3 peaks in concentration is the region of the stratosphere where isotopically altered O_2 is produced, since this process is ultimately due to UV photolysis of ozone (see discussion above in section <ref>). The potential temperature of this region is 1000 K in each of the simulations with ozone, so we use this value as the “source” altitude. We do not account for the fact that the polar night region of the stratosphere would not be photochemically active, since less than 10% of the stratosphere is in polar night at a given time, and this is a rough estimate of Θ. The same value is used to calculate Θ in the simulation without ozone for ease of comparison, despite the fact that a stratosphere lacking ozone would not produce this fractionation signal. We use the no ozone simulation as a limiting case to study the properties of stratospheres with reduced ozone abundances. We determine the “sink” altitude's potential temperature in each case by calculating the spatial and temporal average of the potential temperature at the tropopause. The rate of change of ξ in a region with a given potential temperature θ can be calculated by multiplying each mass flux into and out of that region by the ξ of the region the flux comes from, giving the fluxes of “photochemically altered” mass into and out of the region with that potential temperature from the regions directly above and below it. This allows us to write down a conservation equation at each potential temperature:


    dξ (θ)/dt   = F_above×ξ(θ + dθ)-F_above×ξ(θ)+ F_below×ξ(θ-dθ) - F_below×ξ(θ)


where F_above is the flux of mass between stratospheric regions with potential temperatures θ  and θ + dθ and F_below is the flux of mass between regions with potential temperatures θ and θ-dθ. Each F is calculated as the difference between the maximum and minimum of the STEM circulation at the boundary between two regions. Figure <ref> shows a schematic of a stratospheric layer with potential temperature θ exchanging mass with layers with potential temperatures θ+dθ and θ-dθ via F_above and F_below. Assuming steady-state, we can then derive an equation for ξ:


    ξ(θ)    = F_above×ξ(θ + dθ) + F_below×ξ(θ-dθ)/F_above+F_below


 We solve this system of equations (equation (<ref>) for each zonal- and time-averaged potential temperature) using matrix inversion, with the boundary conditions that ξ is equal to zero immediately below the tropopause (because ξ is a measure of the fraction of mass that has reached a given point directly from the upper stratosphere without passing through the troposphere) and ξ(1000 K) – where the prescribed ozone peaks – is equal to one. Multiplying the value of ξ at the tropopause by the scaling factor A gives Θ.



§ RESULTS



We begin by applying the STEM circulation, including its extension with advection, eddy momentum flux convergence, and diffusion, to reanalysis data and the modern simulation. Next, we quantify changes in the stratospheric circulation between the modern and snowball simulations. We then connect the circulation changes to changes in wave driving, the stratospheric polar vortex, sudden stratospheric warmings, and stratospheric temperatures. Finally, we compare the values of γ and Θ under modern and snowball conditions.



 §.§ Stratospheric circulation



Since this is the first time the STEM framework and its extension have been used to analyze the stratospheric circulation, we apply it to the ERA-Interim reanalysis dataset (see e.g. <cit.> for discussion of the ERA-Interim data; see Fig. <ref> for our decomposition of the data). The full STEM circulation, Ψ_STEM, successfully reproduces the dry isentropic circulation in the stratosphere and all the features of the TEM stratospheric circulation (see Figure 5 in <cit.>). Further, through the decomposition into different dynamical components (see equations <ref> to <ref>), we see that the stratospheric circulation is largely driven by eddies (see Fig. <ref>, Ψ_eddy & Ψ_EMFC), as expected <cit.>. Unresolved processes, possibly gravity wave drag or analysis increments associated with the data assimilation scheme, are responsible for a non-negligible portion of the stratospheric circulation between 0^∘ and 30^∘ (see Fig. <ref>, Ψ_F). We see similar results for the modern simulation (Fig. <ref>). The simulation displays a slightly stronger atmospheric circulation than the ERA-Interim reanalysis data because it lacks ocean heat transport (see e.g. <cit.>).

Consistent with previous work by <cit.>, the winter-hemisphere Hadley cell is much stronger in the snowball simulations than in the modern case (Tables <ref>, <ref>; Figs. <ref> & <ref>). In contrast, the stratospheric circulation is notably weaker in both summer and winter (Tables <ref>, <ref>, <ref>; Figs. <ref> & <ref>). Increasing CO_2 from 280 ppmv to 10,000 ppmv in the snowball simulation strengthens the stratospheric flow (Fig. <ref>). Removing ozone in the snowball simulation increases the depth of the stratospheric circulation (Fig. <ref>) and raises the height of the tropopause, calculated according to the standard WMO definition. Flattening the topography further weakens the snowball stratosphere's circulation and shifts Ψ_EMFC to lower altitudes compared to the mountainous simulations (Fig. <ref>).

Largely consistent with one of the main results of <cit.>, Ψ_F is a major driver of Ψ_mean in the troposphere of our snowball simulations, accounting for a substantial portion in the winter cell, and nearly the entire mean circulation in the summer cell (see Tables <ref> & <ref>; the summer cell is too weak to be visible in the snowball circulation figures). The situation is similar in the winter stratosphere, with Ψ_F accounting for between 22% and 80% of the strength of Ψ_STEM at its extremum (Table <ref>). Ψ_eddy is also a major driver of the winter snowball stratosphere's circulation, especially poleward of 30^∘ (see Fig. <ref>), but specifically at the point where Ψ_STEM reaches its  maximum in the  stratospheric overworld, Ψ_F plays a large role. In the summer hemisphere snowball stratosphere, Ψ_STEM is dominated by Ψ_eddy, but Ψ_F is the main component of Ψ_mean (the summer circulation is too weak to be seen in the figures of snowball stratospheric circulation at the contour interval we used; see Table <ref>). Ψ_F seems to be responsible for some of the highest-altitude transport that occurs in the snowball stratosphere (Fig. <ref>).


Most of the reduction in the strength of the snowball stratosphere's circulation can be attributed to changes in the wave-driven components. Planetary-scale eddy-driven circulation is strongly reduced in the 280 ppmv CO_2 snowball with modern ozone and topography compared to the modern simulation (see Fig. <ref>). This is likely at least partly due to a reduction in planetary-scale Rossby waves because of the lack of land-sea contrast on a fully ice-covered snowball. The flattened snowball simulation displays even weaker planetary-scale wave-driven circulation because of a lack of Rossby waves that are generated as wind blows over topography (Fig. <ref>). It is also clear from Fig. <ref> that planetary-scale eddies contribute to the “deep branch” of the stratospheric circulation and synoptic-scale eddies to the “shallow branch” in both modern and snowball simulations. The synoptic-scale component of Ψ_eddy also weakens in the snowball simulations, particularly in the summer hemisphere, although not by as much as the planetary component. The “shallow” synoptic-scale component of the flattened snowball's circulation is more vigorous than that of the snowball with modern topography (Fig. <ref>).








 §.§ Wave activity



When considering the polar winter stratosphere as a box, the reduction in upward wave propagation in the snowball relative to modern leads to decreased EP flux convergence in the box because the flux out of the subtropical boundary is mostly unchanged in the simulations with ozone (Fig. <ref>). The simulation without ozone displays a similar convergence to the other snowball simulations despite considerably weaker influx across the bottom boundary because its subtropical boundary outflux is also strongly reduced. Since the stratospheric circulation is driven by the EP flux convergence, a reduction of vertical EP flux with small changes in meridional EP flux weakens the circulation as shown in Fig. <ref>.













 §.§ Stratospheric jets


The stratospheric polar vortex occurs due to the potential temperature gradient generated by solar heating of stratospheric ozone in the tropics and strong radiative cooling by CO_2 in the polar night region, and the strength of the vortex is modulated by the convergence of EP flux <cit.>. Because EP flux convergence acts as a torque on the zonal-mean flow and decelerates it, the weakened convergence over the polar cap in the snowballs with modern ozone is consistent with a strengthened polar vortex compared to the modern simulation (right column of Fig. <ref>). In contrast, the snowball simulation without ozone does not display a vortex, since it lacks the temperature gradient generated by ozone heating.

The differences in jet speeds among the three snowball simulations with modern ozone are also consistent with the modulation of vortex strength by EP flux convergence. The small, 1 kg m s^-4 increase in EP flux convergence between the 280 ppmv and 10,000 ppmv CO_2 snowball simulations (compare the blue and orange circles in  Fig. <ref>) mirrors a small, ∼10 m s^-1 reduction in jet speed (Fig. <ref>). The larger 2.5 kg m s^-4  decrease in EP flux convergence between the 280 ppmv CO_2 snowball simulation and the no mountain snowball simulation (compare the blue and green circles in Fig. <ref>) accompanies a 30 m s^-1 increase in jet speed (Fig. <ref>). 

EP flux convergence over the polar cap is also responsible for creating sudden stratospheric warming (SSW) events where the zonal-mean zonal wind reverses sign through deceleration of the polar vortex <cit.>. In <cit.>, an event is defined to occur when the zonal-mean zonal winds at a latitude of 60 and a pressure of 10 hPa become easterly during the period of November to March. At least twenty days must pass after an event before another can be identified, and winds must return to westerly for at least ten days prior to April 30th to avoid classification as a “final warming.” <cit.> find that the modern Earth experiences approximately 0.62 SSWs per year, while our modern simulation produces about one SSW per year, suggesting that the simulation has a bias toward weak zonal winds. The frequency of SSWs goes from about one per year for the modern simulation to none at all for the snowball simulations, which is striking, given the bias in the modern simulation. However, with only 5 years of simulation data, we cannot rule out the possibility that SSWs might occur at low frequencies in the snowball simulations. This is similar to the situation in the modern Earth's Southern Hemisphere stratosphere, where weak wave driving leads to a polar vortex that is stronger than in the Northern Hemisphere, with very infrequent SSWs <cit.>.


 §.§ Stratospheric Temperature


The weakening of the stratospheric circulation, vertical EP flux, and  EP flux convergence as well as the strengthening of the stratospheric polar vortex all imply colder stratospheric temperatures. In particular, in the absence of strong adiabatic warming induced by a vigorous stratospheric circulation, we would expect the temperatures to be closer to radiative equilibrium. This is indeed the case in the snowball simulations with modern ozone (left column of Fig. <ref>). At any given time more than 65% of the Earth's surface in the snowball simulations with modern ozone is beneath stratospheric air that has a temperature less than 195K, the threshold for type 1 polar stratospheric cloud formation <cit.>. In contrast, only around 30% of the modern simulation's surface lies beneath stratospheric air colder than 195K, and none falls below 188 K, the threshold for type 2 polar stratospheric cloud formation <cit.>, except occasionally in the Antarctic during polar winter. In the snowball simulation without ozone, the entire stratosphere falls below 188 K because there is no solar heating of ozone to warm it.



 §.§ Cross-tropopause mass exchange rate, stratospheric mixing efficiency, and Δ^17O


The cross-tropopause mass exchange rate γ of 0.12 yr^-1 calculated from our modern simulation (Table <ref>) is close to the γ of 0.13 yr^-1 calculated from data in <cit.> and used in <cit.>. In turn, the modern Θ of 0.018 is similar to the value of 0.017 used by Cao and Bao (by construction, see section <ref>)). Each snowball simulation has a γ that is smaller than that of the modern simulation, while each snowball simulation has a Θ value that is larger than that of the modern simulation (Table <ref>). 

<cit.> estimated the CO_2 concentration during and after the snowball from Neoproterozoic Δ^17O by assuming modern γ and Θ. Since Δ^17O is a function of the product of γ and Θ (see equation (<ref>)) the decreases in γ and increases in Θ between the modern simulation and snowball simulations end up partially canceling. The largest deviation from modern in γ×Θ is displayed by the snowball simulation without ozone, with a value 27% larger than modern. To see how the relationship between pCO_2/pO_2 and Δ^17O depends on γ and Θ calculated from our simulations, we vary γ×Θ across the values in Table <ref> while holding the other variables in equation (<ref>) constant at the values used in <cit.> and listed in Table 1 of <cit.> (see our Fig. <ref>). 

The variations in γ×Θ can, depending on the tropospheric Δ^17O value to be explained, produce a maximum increase in estimated CO_2 of up to 91× or a maximum decrease by a factor of 3.6, measured relative to estimates produced with the modern γ×Θ. However, the region in Fig. <ref> that displays the largest variation in pCO_2/pO_2 as a function of γ×Θ is also a region where pCO_2 is orders of magnitude greater than pO_2, which is arguably not likely to have been the case during the Neoproterozoic. <cit.> all estimate that Neoproterozoic oxygen levels may have been equal to or greater than 10% of modern pO_2, implying an atmospheric concentration ≥20,000 ppmv, and <cit.> use boron and calcium isotopes to estimate a maximum late Neoproterozoic CO_2 concentration of 90,000 ppmv. In the region of the parameter space where pCO_2 ≤ 10×pO_2, the values of pCO_2/pO_2 remain close together over the range of γ×Θ produced by our simulations, which helps validate the use of modern γ and Θ to estimate snowball CO_2 in <cit.>. 





§ CONCLUSION & DISCUSSION




 §.§ Conclusions

In this study we analyzed GCM simulations of the stratosphere of the snowball Earth. We focused on changes in the stratospheric circulation, EP flux convergence, polar vortex, and temperatures relative to modern conditions. The stratospheric circulation was examined using the Statistical Transformed Eulerian Mean framework of <cit.>. We extended the STEM framework to include components accounting for advection, eddy momentum flux convergence, and diffusion, since <cit.> demonstrated the importance of diffusion for the snowball Hadley cell.




We found that the stratospheric circulation during a snowball would have been weakened compared to the modern. The reduced circulation strength in the winter stratosphere of the Northern Hemisphere was linked to a weakening of the torque induced by eddy convergence in the stratosphere. The decreased torque was likely caused primarily by a reduction in Rossby waves due to the lack of land-sea contrast in the snowball simulations. In the snowball simulations, wave activity at both planetary and synoptic scales was generally weakened relative to modern. The weakening of eddy convergence led to a much stronger polar vortex in snowball simulations with modern ozone, while the lack of differential stratospheric solar heating in the case without ozone eliminated the vortex entirely. Reduced eddy convergence also eliminated sudden stratospheric warming events in all snowball simulations. In the simulations with ozone, the weak eddy heat flux led to extremely low stratospheric temperatures, particularly over the poles where the powerful vortex isolated the stratospheric air. Since the simulation without ozone had much less UV absorption, the temperatures in the stratosphere were even lower.
 
The anomalous Δ^17O signal discovered in Neoproterozoic sulfates and interpreted as evidence for elevated CO_2 during a snowball in <cit.> based on the model presented in <cit.> is dependent on the product of the stratosphere-troposphere exchange (γ) and the stratospheric mixing efficiency (Θ). In this study, we quantified γ and Θ under snowball conditions from GCM output for the first time. We showed that reductions in γ and increases in Θ under snowball conditions partially cancel one another in our simulations. This makes intuitive sense, as a reduction in the flux of material from the troposphere into the stratosphere (a smaller γ) should make it easier to homogenize the isotopic composition of the stratosphere (an increase in Θ) and vice versa. This is most clearly illustrated by considering the limit where γ is zero–the composition of the stratosphere would homogenize without a flux of material from the troposphere, which is equivalent to having a Θ of one (perfect stratospheric mixing efficiency). The use of modern values for those variables in <cit.> should not introduce large errors in estimates of snowball CO_2 from Δ^17O, unless there was orders of magnitude more CO_2 than O_2 in the snowball atmosphere. There are larger uncertainties in the estimate of pCO_2 from Δ^17O due to the lack of information about what fraction of oxygen in the mineral samples came from the atmosphere (see <cit.>). 


 §.§ Discussion

Because our models did not have interactive chemistry, the implications of our findings for Δ^17O are tentative. The photochemistry that generates the fractionation signal is determined by the temperatures, bulk velocities, and spatial distributions of participating O_3, O_2, and CO_2 molecules. Given the extreme cold, powerful jet (in the modern ozone cases) or weakened jet (in the no ozone case, which can be considered an end-member of reduced ozone scenarios), and weakened Brewer Dobson circulation, a detailed quantification of the effects of changes to these properties on estimates of the magnitude of Δ^17O is warranted. This will better constrain how CO_2 concentration influenced oxygen fractionation under the exotic conditions of the snowball. Chemical transport models would be useful in exploring the evolution of the coupled O_3-O_2-CO_2 photochemical network during the Neoproterozoic glaciations.

The extremely cold temperatures in the stratospheres of our snowball simulations might allow for the formation of “polar stratospheric clouds” (PSCs) far from the poles (see e.g. <cit.> for discussion of PSCs). If present, optically thick PSCs would provide a strong greenhouse effect <cit.>. If conditions allowed for the formation and persistence of optically thick PSCs, especially over the tropics, then their greenhouse effect might have contributed to deglaciation, similar to the role of tropospheric clouds explored in <cit.>.



The changes in the simulated stratosphere of the snowball Earth are analogous but opposite to changes to the modern stratosphere that are predicted to occur in response to warming due to anthropogenic climate change. For example, the snowball involves weaker Brewer Dobson circulation, weaker eddy fluxes, zero SSWs, a stronger stratospheric vortex (under modern ozone), and colder temperatures compared to modern. In contrast, during warmer climates, the circulation is predicted to accelerate due to increased eddy fluxes <cit.>. In a warmer climate there are competing effects of a colder stratosphere and more wave driving <cit.>; however, in a model with a strengthened Madden Julian Oscillation, the stratosphere is predicted to exhibit a weaker polar vortex, more sudden stratospheric warming events, and warmer temperatures <cit.>. 


Given the bi-directional coupling between stratospheric and tropospheric dynamics that has been observed and modeled under modern conditions, an altered snowball stratosphere may impact deglaciation scenarios and geochemical proxy interpretation. Interesting directions for future investigation include the modeling of stratospheric cloud formation and the inclusion of dynamical ozone chemistry in 3D simulations of the snowball atmosphere. It would also be useful to examine the stratospheric properties of partially glaciated climates like the Jormungand <cit.>, the Slushball <cit.>, or a stable “Waterbelt” <cit.>, which might have significantly different dynamics and trace chemical abundances, with the potential to influence interpretation of the Δ^17O CO_2 proxy.  








 This work was completed using resources provided by the University of Chicago Research Computing Center. We thank two anonymous reviewers and Dann Mitchell for their insightful comments. We acknowledge support from the University of Chicago College Research Fellows Program, NASA grant number NNX16AR85G, which is part of the “Habitable Worlds” program, the Packard Foundation, the Clarendon Scholarship, and Modalert. Data to reproduce figures can be found at DOI:10.6082/uchicago.1920.


agufull08

































































































