Hydrological Cycle Changes Explain Weak Snowball Earth Storm Track Despite Increased Surface Baroclinicity

Simulations show that storm tracks were weaker during past cold, icy climates relative to the modern climate despite increased surface baroclinicity. Previous work explained the weak North Atlantic storm track during the Last Glacial Maximum using dry zonally asymmetric mechanisms associated with orographic forcing. Here we show that zonally symmetric mechanisms associated with the hydrological cycle explain the weak Snowball Earth storm track. The weak storm track is consistent with the decreased meridional gradient of evaporation and atmospheric shortwave absorption and can be predicted following global mean cooling and the Clausius‐Clapeyron relation. The weak storm track is also consistent with decreased latent heat release aloft in the tropics, which decreases upper tropospheric baroclinicity and mean available potential energy. Overall, both hydrological cycle mechanisms are reflected in the significant correlation between storm track intensity and the meridional surface moist static energy gradient across a range of simulated climates between modern and Snowball Earth.


Introduction
Storm tracks are regions where midlatitude cyclones occur most frequently and several modern theories connect their intensity to (near-) surface baroclinicity (meridional temperature gradient, Chang et al., 2002;Held, 2018;Shaw et al., 2016). The surface baroclinicity-intensity connection is supported by the seasonality of storm track intensity in modern reanalysis data (see Figure S5 in O'Gorman, 2010) and idealized simulations of annual mean storm track intensity across a range of climates without ice (see Figure 4 in O'Gorman & Schneider, 2008) (see Figure 1 in Caballero & Langen, 2015). It is also assumed in energy balance models (EBMs) based on surface temperature (Mbengue & Schneider, 2018;North, 1975) and surface moist static energy (Frierson et al., 2007). However, it has been known for some time that the surface baroclinicity-intensity connection fails when applied to simulations of past cold, icy climates such as the Last Glacial Maximum (LGM; Hall et al., 1996;Li & Battisti, 2008) and Snowball Earth (Pierrehumbert, 2005). The Snowball Earth hypothesis proposes that during Neoproterozoic glaciations, ice covered most of the planet, possibly to the equator (Hoffman et al., 2017).
Several mechanisms have been proposed to explain why the Northern Hemisphere (NH) wintertime Atlantic storm track during the LGM is weaker than modern despite increased surface baroclinicity. All of the mechanisms are based on dry zonally asymmetric dynamics; for example, they appeal to changes in stationary wave amplitude (Hall et al., 1996), reflection (Lofverstrom et al., 2016), and structure (Riviere et al., 2018) due to the orographic forcing of the continental ice sheets or changes in upstream seeding (Donohoe & Battisti, 2009). Zonal asymmetries are weaker in Snowball Earth simulations (Graham et al., 2019); however, the weakening of the NH wintertime zonal mean Snowball Earth storm track is too large (∼4 PW, see Figure 7 in Pierrehumbert, 2005) to be explained by zonally asymmetric LGM mechanisms. Thus, alternative mechanisms are needed.
The hydrological cycle offers alternative zonally symmetric mechanisms to explain the weak Snowball Earth storm track despite increased surface baroclinicity. The cold, icy Snowball Earth conditions decrease latent heat release aloft, decrease evaporation, and increase surface albedo. The decrease of latent heat release aloft affects atmospheric thermal structure (static stability and meridional temperature gradient aloft) and mean available potential energy (MAPE), which impacts storm track intensity as measured by eddy kinetic energy (EKE). O'Gorman and  and Schneider et al. (2010) showed that storm track intensity changes from warm and wet to cold and dry climates in idealized simulations were controlled primarily through the impact of latent heat release on static stability. Thus, the weak Snowball Earth storm track may be connected to static stability changes rather than surface baroclinicity changes. However, Yuval and Kaspi (2020) showed that static stability changes can oppose surface baroclinicity resulting in a small net storm track intensity change in dry idealized simulations, highlighting the difficulty of predicting storm track changes based on static stability alone.
The transition from ocean and land to ice everywhere decreases evaporation (surface latent heat flux) and increases surface albedo, which affects the moist static energy (MSE) budget of the atmosphere. Shaw et al. (2018) developed an MSE framework for zonal mean storm track intensity and showed that the increased meridional latent heat flux gradient in response to global mean surface warming strengthens the NH wintertime storm track. Tan and Shaw (2020) used the locking technique to confirm that global mean warming dominates the meridional latent heat flux gradient response rather than the deviation from global mean warming (i.e., the surface baroclinicity response). Thus, the weak Snowball Earth storm track may be connected to the decreased meridional latent heat flux gradient in response to global mean surface cooling.
Here we quantify whether zonally symmetric mechanisms associated with the hydrological cycle can explain why the Snowball Earth storm track is weaker than modern despite increased surface baroclinicity. We simulate a hard Snowball Earth (ice everywhere) and quantify the impact of the hydrological cycle via changes in surface boundary conditions using the MSE framework and atmospheric thermal structure using MAPE. We predict a weaker Snowball Earth storm track assuming that hydrological cycle changes follow global mean cooling and the Clausius-Clapeyron relation. Finally, we quantify the connection between storm track intensity and other variables such as meridional surface MSE gradient and upper tropospheric baroclinicity.

MSE Framework
The MSE framework for storm track intensity can be derived from the zonal mean atmospheric MSE budget with the global mean removed: Kang et al., 2008), where F TE = ⟨[v ′ m ′ ]⟩ and F SC = ⟨[vm]⟩ are the MSE flux by transient eddies and stationary circulation (mean meridional circulation plus stationary eddies), respectively, ∇ ⋅ F NE = NE is in flux form where NE is the net energy input with its global mean removed, ⟨⋅⟩ denotes a mass-weighted vertical integration, [⋅] denotes a zonal average, and⋅ denotes a monthly average, that is, an average over a particular month and not the climatological monthly average. Similar results are obtained using a 10-day high-pass filter (Shaw et al., 2018). The global mean is removed to emphasize meridional gradients (see Donohoe & Battisti, 2012;Donohoe et al., 2020;Kang et al., 2008;Shaw et al., 2018).
An equation for storm track intensity I (unit PW) is derived by multiplying (1) by 2 a 2 cos where a is the radius of the Earth and is latitude and integrating between the pole and the storm track position s (where ∇ ⋅ F TE = 0):

Geophysical Research Letters
10.1029/2020GL089866 (3) (Shaw et al., 2018). Consequently, a storm track intensity change ( I) is connected to changes in net energy input ( I NE ) or MSE flux by the stationary circulation ( I SC ): In order to quantify and predict the impact of hydrological cycle changes on storm track intensity, we decompose net energy input as follows: where THF is the turbulent heat flux (sensible SH plus latent LH heat flux), SWABS is the shortwave absorption (TOA minus surface shortwave radiation), GHE is the greenhouse effect (surface minus TOA longwave radiation), h/ t is atmospheric storage where h = c p T + L v q is specific enthalpy and all terms are written in flux form, for example, ∇ ⋅ F THF = THF and THF has its global mean removed. The storm track interacts with the terms on the right hand side of (5); thus, it is difficult to infer causality when the MSE framework is used diagnostically. However, Shaw et al. (2018) and Barpanda and Shaw (2020) showed that external parameters in the MSE framework (e.g., insolation and mixed layer depth) can form the basis of predictions, scaling estimates, and thereby causal understanding.
We can predict the impact of hydrological cycle changes on storm track intensity between the Snowball Earth and modern climates using the MSE framework and the Clausius-Clapeyron relation. The fractional humidity change between different climates following the Clausius-Clapeyron relation is where q * s is the saturation specific humidity, e s is the saturation water vapor pressure, and the subscripts M and S refer to the modern and Snowball Earth climates, respectively (Boos, 2012;Held & Soden, 2006;Schneider et al., 2010). Using representative global mean surface temperatures for the Snowball Earth (T S ≈ 245 K) and modern (T M ≈ 285 K) climates, the humidity decrease implied from (7) is approximately 95% ( ≈ −0.95).
Here we predict the impact of the hydrological cycle via surface boundary condition (latent heat flux and surface albedo) changes on Snowball Earth storm track intensity. While we expect that hydrological cycle changes will also impact the greenhouse effect (via a colder global mean temperature), it is not clear how the meridional structure of temperature will change because it is affected by the storm track itself. Thus, we cannot easily predict I GHE . Along similar lines, it is not easy to predict I SC and I h/ t . Hence, we assume that I GHE , I SC , and I h/ t are small and focus on predicting I THF and I SWABS . We use the simulations discussed below to justify our assumptions and test our predictions by diagnosing all the terms in the MSE framework.
Assuming that the change in turbulent heat flux between the Snowball Earth and modern climates is dominated by latent heat flux and the latent heat flux follows changes in saturation specific humidity, then where the p subscript refers to prediction and H is relative humidity (Schneider et al., 2010). The predicted intensity change is where ∇ ⋅ F LH,p = LH p and LH p has its global mean removed. Since < 0, the meridional gradient of LH p is positive (Figures S1a and S1b) implying I THF, p < 0, that is, a weakening of the storm track.

10.1029/2020GL089866
Shortwave absorption depends on downward TOA ( S ↓ T ) and surface ( S ↓ S ) shortwave radiation and the surface ( s ) and planetary ( p ) albedos, that is, Assuming that Snowball Earth shortwave absorption depends on surface ice albedo, shortwave optical depth [ S = − ln(S ↓ S,S ∕S ↓ T,S )], and unchanged TOA shortwave radiation, that is, The predicted intensity change is where ∇ ⋅ F SWABS,p = SWABS p and SWABS p has its global mean removed. If we assume that the Snowball Earth shortwave optical depth decreases following the Clausius-Clapeyron relation, that is, S ≈ (1 + ) M = 0.05 M , and the surface ice albedo is s,S = 0.6 in AGCM and s,S = 0.8 in AQUA, then the meridional gradient of SWABS p is positive (Figures S1c and S1d) implying I SWABS, p < 0, that is, a weakening of the storm track. We note that SWABS p is dominated by the change in shortwave optical depth. For both predictions (10) and (12), the storm track position is assumed to be fixed to its modern value following previous work (Shaw et al., 2018).

MAPE Framework
Storm track intensity defined using EKE is linearly related to MAPE in idealized simulations and reanalysis data (Gertler & O'Gorman, 2019;O'Gorman, 2010;O'Gorman & Schneider, 2008). In the simulations discussed in section 2.3, EKE is linearly related to MAPE, which we define following equation (3) in O'Gorman and Schneider (2008): where = R∕c p , R is dry gas constant, c p is specific heat at constant pressure, g is gravitational acceleration, p 0 = 1, 000 hPa is a reference pressure, {⋅} is an average over the NH extratropics (20 • to 90 • N), p t is the WMO tropopause pressure, and is potential temperature. We integrate from 850 hPa to p t because it yields good agreement with the parcel-moving algorithm of Stansifer et al. (2017, see Figure S2).
Following O'Gorman and Schneider (2008), we approximate the variance as where y (⋅) ≡ (⋅)/ /a is the meridional gradient in spherical coordinates and L Z is the meridional width of the baroclinic zone. Accordingly, MAPE is approximated as and changes can be decomposed into baroclinicity ( Baro.) and stability ( Stab.) contributions, that is, Since MAPE depends on the time and zonal mean temperature throughout the atmosphere, it cannot be predicted a priori for Snowball Earth.

Simulations
We focus on the NH wintertime (December, January, and February) storm track and simulate hard Snowball Earth and modern climates using two configurations of the ECHAM6 general circulation model (Stevens et al., 2013). The first configuration is a slab ocean atmospheric general circulation model, hereafter referred to as AGCM. We simulate the modern climate with modern topography, obliquity, greenhouse gases (CO 2 = 280 ppmv), a 50-m mixed layer depth, and zero ocean energy transport. We simulate a hard

10.1029/2020GL089866
Snowball Earth by imposing ice everywhere (ocean covered with sea ice and land covered with glaciers) in the modern simulation with an ice albedo of 0.6. The AGCM simulations are identical to those in Graham et al. (2019). Climate change in the AGCM is the difference of the Snowball Earth and modern simulations.
The second configuration is a slab ocean aquaplanet general circulation model with modern obliquity, greenhouse gases, zero ocean energy transport, and thermodynamic sea ice (a motionless single slab, Giorgetta et al., 2012), hereafter referred to as AQUA. The modern AQUA simulation has a 50-m mixed layer depth (Barpanda & Shaw, 2020;Donohoe et al., 2014). Here we vary the mixed layer depth as a simple way of simulating a range of climates between modern and Snowball Earth. In particular, we vary the mixed layer depth from 50 to 5 m for a total of 14 simulations (we use AQUA data from both wintertime hemispheres for a total of 28). When the mixed layer depth is ≲20 m, we obtain a hard Snowball Earth with sea ice at the equator ( Figure S3) because the surface heat capacity becomes small enough to trigger a runaway ice-albedo feedback. The ice-albedo is not prescribed in AQUA and can be greater than 0.6 because of the presence of snow on ice. Climate change in AQUA is the difference of Snowball Earth (15-m mixed layer depth) and modern (50-m mixed layer depth) simulations. The difference between the other mixed layer depth simulations and the modern simulation captures a range of climates between Snowball Earth and modern.
Since it is very common in the literature to perform slab ocean aquaplanet simulations without thermodynamic sea ice (temperature can be below freezing without ice forming, Bordoni & Schneider, 2008;Kang et al., 2008;Lee et al., 2008;O'Gorman & Schneider, 2008), we perform a 15-m mixed layer depth simulation without sea ice and with a surface (ocean) albedo of 0.06 ( Figure S4a) to quantify the importance of icy boundary conditions.

Simulated Snowball Earth Climate
Relative to the modern climate, the Snowball Earth simulations have more ice (surface albedo ≥ 0.6, Figure  1a), larger NH surface baroclinicity (equator minus pole surface temperature difference increases by 27 K in AGCM and 29 K in AQUA, Figure 1b The importance of icy boundary conditions for the Snowball Earth climate is quantified by comparing the 15-m mixed layer depth AQUA simulations with and without sea ice. When sea ice is disabled, surface baroclinicity and storm track intensity are both weaker than the modern simulation ( Figure S4). This demonstrates that a weak storm track-increased surface baroclinicity climate change can only be achieved with ice.

Impact of Hydrological Cycle via Changes in Surface Boundary Conditions
According to the MSE framework, the weaker AGCM Snowball Earth storm track ( I, Figure 2a Figure 2a) also strengthens the storm track slightly consistent with stationary eddy storm track compensation (see Barpanda & Shaw, 2017;Manabe & Terpstra, 1974, Figures 1d and 1e). The contribution from atmospheric storage is small (− I h/ t , Figure 2a). Similar results are seen throughout the NH ( Figure S5).
The weaker AQUA Snowball Earth storm track is also consistent with decreased turbulent heat flux ( I THF , Figure 2b) and shortwave absorption ( I SWABS , Figure 2b). Latent heat flux ( I LHF , Figure 2b) also dominates over sensible heat flux (difference between I THF and I LHF , Figure 2b). The larger shortwave absorption contribution in AQUA as compared to AGCM (compare I SWABS , Figures 2a and 2b) is discussed below. The greenhouse effect contribution in AQUA is opposite of that in AGCM (compare − I GHE , Figures 2a and 2b) due to land and clouds ( Figure S6). Finally, the stationary circulation contribution in AQUA is  Figures 2a and 2b). The overprediction occurs because the latent heat flux decrease is less than 95% in the Southern Hemisphere ( > −0.95) and the change in sensible heat flux is non-negligible (Figures S1a and S1b). Assuming the shortwave optical depth decreases following the Clausius-Clapeyron relation [ S ≈ 0.05 M , see (12)] and there is a surface ice albedo everywhere overpredicts the simulated weakening in AGCM (compare I SWABS to I SWABS, p , Figure  2a) but accurately predicts it in AQUA (compare I SWABS to I SWABS, p , Figure 2b). The AGCM overprediction is once again associated with the response in the Southern Hemisphere (Figures S1c and S1d) where the decrease of summertime shortwave optical depth is less than 95% of its modern value (i.e., S > 0.05 M ) and there are important meridional changes in planetary albedo.
Across the range of climates between modern and Snowball Earth in AQUA, the storm track intensity change follows the turbulent heat flux contribution ( I THF , Figure 2c). The shortwave absorption contribution exhibits a nonlinear relationship with intensity changes ( I SWABS , Figure 2c) and the greenhouse effect contribution exhibits the opposite sign (− I GHE , Figure 2c). Thus, turbulent heat fluxes are not only important for the weakening of the storm track between the Snowball Earth and modern climates but are also important across a range of climates in between. Overall, the cold, icy Snowball Earth conditions decrease the meridional gradient of latent heat flux and shortwave absorption, which leads to a weaker meridional net energy input gradient and thereby a weaker storm track. The weakening can be predicted following global mean cooling and the Clausius-Clapeyron relation and assuming a surface ice albedo everywhere because the greenhouse effect, stationary circulation, and atmospheric storage contributions are small. The weaker storm track is also consistent with weaker meridional surface (Figure 1c) and lower tropospheric (Figures 2d and 2e) MSE gradients, which also follow the Clausius-Clapeyron relation as discussed below.

Impact of Hydrological Cycle via Changes in Thermal Structure
Next, we examine how the cold, dry Snowball Earth conditions impact storm track intensity, as measured by EKE, via changes in atmospheric thermal structure that affect MAPE. The weaker storm track in the AGCM Snowball Earth simulation relative to modern is consistent with weaker MAPE ( MAPE, square, Figure 3a), which is dominated by the change in baroclinicity ( Baro., square, Figure 3a) rather than stability ( Stab., square, Figure 3a). When examining the range of climates between modern and Snowball Earth in AQUA, the storm track intensity change relative to modern is significantly correlated with MAPE ( R = 0.98, MAPE, Figure 3a) and is mostly dominated by changes in baroclinicity ( Baro., Figure 3a) rather than stability ( Strat., Figure 3a).
Since surface baroclinicity does not account for the weaker Snowball Earth storm track relative to the modern climate (Figure 1b), the MAPE results suggest that upper tropospheric baroclinicity is important. Indeed, the temperature difference between the Snowball Earth and modern simulations involves weaker upper tropospheric baroclinicity connected to greater tropical cooling aloft (Figures 3b and 3c).

Connecting Storm Track Intensity to Other Variables
The surface baroclinicity-intensity connection is appealing because it connects storm track intensity (a turbulent quantity) to surface baroclinicity (a mean quantity), which can be inferred from paleoproxy data. However, storm track intensity is not significantly correlated with surface ( R = 0.01, Figure 4a) or near-surface ( R = 0.20, Figure 4b) baroclinicity across the AQUA and AGCM simulations.
The hydrological cycle changes suggest that surface latent heat flux and latent heat release aloft (upper tropospheric baroclincity) may be more important than surface temperature. Recall that the weaker Snowball Earth storm track is consistent with a weaker meridional gradient of latent heat flux and shortwave absorption, which follow global mean cooling and the Clausius-Clapeyron relation [see (7) -(12)]. The meridional surface MSE gradient in Snowball Earth also follows the Clausius-Clapeyron relation ( Figure S7a). Consistently, storm track intensity is significantly correlated with equator minus pole surface MSE ( R = 0.94, Figure 4c) and equator minus pole surface MSE estimated from the Clausius-Clapeyron relation ( R = 0.94, Figure S7b) across the AQUA and AGCM simulations. Furthermore, equator minus pole surface MSE is significantly correlated with equator minus pole surface latent heat flux ( R = 0.96).
Storm track intensity is also significantly correlated with upper tropospheric baroclinicity across the AQUA and AGCM simulations ( R = 0.98, Figure 4d). Upper tropospheric baroclinicity is significantly correlated with equator minus pole surface MSE ( R = 0.98) consistent with the impact of surface moisture on latent heat release aloft in the tropics via moist adiabatic adjustment ( Figure S7c). Upper tropospheric and surface baroclinicity are not significantly correlated ( R = 0.10). Note that in midlatitudes, the storm track itself shapes the temperature structure aloft in addition to convection.
Overall, the meridional surface MSE gradient encapsulates both mechanisms associated with hydrological cycle changes between the Snowball Earth and modern climates. The significant correlation between storm track intensity metrics (transient eddy MSE flux or EKE) and the meridional surface MSE gradient shows that surface moisture changes are important for determining Snowball Earth storm track intensity relative to the modern climate and can dominate over changes in surface baroclinicity.

Conclusions and Discussion
Snowball Earth is an exotic climate with a radically different hydrological cycle than the modern Earth. Here we examined whether zonally symmetric mechanisms associated with the hydrological cycle can account for the weaker Snowball Earth storm track despite increased surface baroclinicity. The results show that the weak Snowball Earth storm track is consistent with the decreased meridional gradient of turbulent heat flux (evaporation) and shortwave absorption (shortwave optical depth), which decreases the meridional gradient of net energy input. The weaker storm track can be predicted assuming a surface ice albedo everywhere and the hydrological cycle changes follow global mean cooling and the Clausius-Clapeyron relation. The weaker storm track is also consistent with decreased latent heat release aloft in the tropics, which decreases MAPE and upper tropospheric baroclinicity via convective adjustment. Overall, both mechanisms (decreased meridional gradient of latent heat flux and shortwave absorption and decreased latent heat release aloft) are reflected in the significant correlation between storm track intensity and the meridional surface MSE gradient across a range of climates between Snowball Earth and modern. They are also reflected in the significant correlation between storm track intensity defined using transient eddy MSE flux Geophysical Research Letters 10.1029 or EKE ( R = 0.99). Ultimately, the results show that the decreased meridional surface moisture gradient between the Snowball Earth and modern climates is important for storm track intensity and dominates over increased surface baroclinicity.
Our results are consistent with Lapeyre and Held (2004) who showed that lower-layer MSE (rather than lower-layer temperature) is most appropriate for diffusive models of energy transport. Since MSE transport reflects both thermodynamic (moisture) and dynamic (EKE) components, future work must focus on understanding the importance of the moisture versus EKE decrease for the decreased Snowball Earth MSE transport, including the connection to diffusivity changes. The results imply that EBMs, which connect storm track intensity to surface baroclinicity, must also include the dependence of storm track intensity on the meridional surface moisture gradient. Hence, EBMs based on the atmospheric MSE budget, that is, an equation for surface MSE (Frierson et al., 2007;Hwang & Frierson, 2010), should be used instead of EBMs based on the TOA energy budget, that is, an equation for surface temperature (Mbengue & Schneider, 2018;North, 1975). Finally, our results are also consistent with previous work that highlighted the importance of upper tropospheric baroclinicity changes for the storm track intensity response to increased CO 2 (Harvey et al., 2015;O'Gorman, 2010).
Here we showed that zonally symmetric changes in the hydrological cycle can explain the weak Snowball Earth storm track despite increased surface baroclinicity. Previous work explained the weak North Atlantic storm track during the LGM using dry zonally asymmetric mechanisms connected to orographic forcing. Assessing whether hydrological cycle mechanisms can also account for the weak LGM North Atlantic storm track requires extending the MSE framework for storm track intensity to the North Atlantic basin as well as quantifying the importance of orography versus latent heat flux and surface albedo changes (cf. Roberts & Valdes, 2017), which is a work in progress. If the Snowball Earth results hold for other paleoclimates, then estimates of paleosurface MSE gradients (combining paleoproxy temperatures and the Clausius-Clapeyron relation) could potentially be used to estimate paleostorm track intensity.

Data Availability Statement
Data and code necessary to reproduce the figures in this paper are available via The University of Chicago's institutional repository Knowledge@UChicago (https://doi.org/10.6082/uchicago.2201). T. A. S. acknowledges support from NSF (AGS-1742944) and R. J. G. acknowledges scholarship funding from the Clarendon Fund and Jesus College, Oxford. The authors thank three anonymous reviewers for their helpful comments. T. A. S. thanks Paul O'Gorman for providing code to calculate MAPE using the parcel-moving algorithm and Osamu Miyawaki for providing code to calculate the moist adiabat. The simulations in this paper were completed with resources provided by the University of Chicago Research Computing Center.