Abstract
Extratropical cyclones are major contributors to consequential weather in the mid-latitudes and tend to develop in regions of enhanced cyclogenesis and progress along climatological storm tracks. Numerous studies have noted the influence that terrestrial snow cover exerts on atmospheric baroclinicity which is critical to the formation and trajectories of such cyclones. Fewer studies have examined the explicit role which continental snow cover extent has in determining cyclones’ intensities, trajectories, and precipitation characteristics. While several examinations of climate model projections have generally shown a poleward shift in storm tracks by the late 21stcentury, none have determined the degree to which the coincident poleward shift in snow extent is responsible. A method of imposing 10th, 50th, and 90th percentile values of snow retreat between the late 20th and 21st centuries as projected by 14 Coupled Model Intercomparison Project Phase Five (CMIP5) models is used to alter 20 historical cold season cyclones which tracked over or adjacent to the North American Great Plains. Simulations by the Advanced Research version of the Weather Research and Forecast Model (WRF-ARW) are initialized at 0 to 4 days prior to cyclogenesis. Cyclone trajectories and their central sea level pressure did not change substantially, but followed consistent spatial trends. Near-surface wind speed generally increased, as did precipitation with preferred phase change from solid to liquid state. Cyclone-associated precipitation often shifted poleward as snow was removed. Variable responses were dependent on the month in which cyclones occurred, with stronger responses in the mid-winter than the shoulder months.
Introduction
The influence of snow cover
Northern hemisphere snow cover is, at its seasonal maximum, the largest component of the terrestrial cryosphere and exerts considerable influence on the mid-latitude atmospheric circulation through a diverse set of mechanisms which have a general cooling effect (e.g. Leathers et al. 1995; Vavrus 2007; Dutla 2011). Near the surface, the presence of snow cover typically lowers air temperature due to the snow’s high albedo (Baker et al. 1992) and its properties as an effective sink of sensible and latent heat (Grundstein and Leathers 1999), which contribute to an increase in static stability (Bengtsson 1980) and a reduction of moisture flux into the atmosphere (Ellis and Leathers 1999). This inhibition of upward moisture flux may be responsible for the negative correlation between snow cover and precipitation observed by Namias (1985) and modelled by Walland and Simmonds (1996) and Elguindi et al. (2005).
Studies have also shown that continental snow cover extent (SCE) is sometimes responsible for modulating upper-level circulation (e.g. Namias 1962; Walland and Simmonds 1996; Cohen and Entekhabi 1999; Gutzler and Preseton, 1997; Notaro and Zarrin, 2011) and that accurately initializing snow cover can improve subseasonal forecast skill considerably (Jeong et al. 2013; Thomas et al. 2016). It is because of this apparent relationship between snow cover and atmospheric circulation that determination of the regional dependence and the temporal scales at which snow cover drives responses in the atmosphere is of fundamental importance to both short- and long-term forecasting.
Observations and hypotheses about the influence of established SCE on the characteristics of ensuing synoptic weather systems may have begun with Lamb (1955). However, one of the first analyses of this relationship was provided by Namias (1962) who hypothesized that the abnormally extensive North American (NA) snow cover of the winter of 1960 had contributed to the more frequent and intense cyclone development observed along the Atlantic coast by enhancing baroclinicity between the continent and the much warmer ocean. Dickson and Namias (1976) subsequently showed that periods of great continental warmth or cold in the American Southeast had a direct influence on the strength of the baroclinic zone near the coast and would affect the average frequency and positions of extratropical cyclones, drawing them further south when the region was colder. Likewise, Heim and Dewey (1984) showed that extensive NA snow cover contributed to a greater frequency of cyclones in the southern Great Plains and Southeast and a reduction in the amount of cyclones tracking further north. From 1979-2010 in NA, a greater frequency of cold season mid-latitude cyclones was observed in a region 50-350 km south of the southern snow extent boundary (snow line) by Rydzik and Desai (2014) who noted a similar distribution of low-level baroclinicity relative to the snow line.
Modeling studies have indicated a similar relationship between snow extent and extratropical cyclone statistics. Ross and Walsh (1986) studied the influence of the snow line on 100 observed North American cyclone cases which progressed approximately parallel to the baroclinic zone within 500-600 km of the snow line. By measuring forecast error from a barotropic model, they were able to determine that the baroclinicity associated with the snow boundary was an important factor in cyclone steering and intensity. Walland and Simmonds (1997) performed global climate model (GCM) experiments with forced anomalously high and low extents of realistic snow cover distributions, ultimately finding a reduction in NA cyclone frequency when snow cover was more extensive, with cyclones frequently occurred further south, similar to the observations of Heim and Dewey (1984). Elguindi et al. (2005) used a 25-km-resolution nested domain over a portion of the Great Plains in the Penn State-National Center for Atmospheric Research (NCAR) Mesoscale Model (MM5) and simulated eight well-developed cyclone cases with snow cover added throughout the domain, initializing 48 hours prior to each cyclone’s arrival to the inner domain. All perturbed cyclone case simulations underwent an increase in central pressure and decrease in total precipitation with slight shifts in the cyclone trajectory which were highly variable and inconsistent. However, this study only used a limited number of cases and only perturbed simulations by adding snow to the entirety of the inner domain rather than altering the position of SCE.
In North America, the net effects of snow cover are nowhere more pronounced than the Great Plains region, which has the highest local maximum of snow albedo (Robinson and Kukla 1985; Jin et al. 2002) and where the strongest correlation between snow cover and negative temperature anomalies has been observed in NA (Heim and Dewey 1984; Robinson and Hughes 1991). The Great Plains region represents one of the largest disparities between local maximum snow albedo and background land surface albedo on the continent, indicating the greatest albedo gradient across a snow line (Figure 1). The land surface is characterized by high inter- and intra-annual snow cover variability (Robinson 1996) and low surface roughness. Winter cyclones track over the Great Plains with high frequency due in part to areas of enhanced cyclogenesis in the lee of the Rocky Mountains (Reitan 1974; Zishka and Smith 1980). The two most prolific cyclogenetic zones over the NA landmass account for the two types of cyclone tracks studied here: the Alberta Clipper track, which typically begins in Alberta, Canada and proceeds to the southeast (Thomas and Martin 2007), and the Colorado Low track, which starts near southeast Colorado and often proceeds northeast toward the Great Lakes region (Zishka and Smith 1980). Because of their spatial extent and great frequency in the region, extratropical cyclones contribute substantially to the hydrology of the Great Plains, accounting for greater than 80% of the total winter (December through February) precipitation throughout much of the region (Hawcroft et al. 2012).
The scope of this study
Typically, simulations of projected future climate states are implemented with global climate models (GCMs) which are limited by expansive resolutions over a global domain and, as a result, do not allow the models to adequately resolve the fine details necessary to accurately reproduce phenomena like precipitation and the diurnal cycle. Harding et al. (2013) demonstrated that dynamically downscaling Coupled Model Intercomparison Project Phase Five (CMIP5) simulations to 30 km resolution in the NCAR WRF model improved simulation of precipitation, especially extreme precipitation events, in the Central U.S. Many modeling studies have applied global and regional climate models to study the projected behavior of extratropical cyclones in the late 21st century (e.g. Maloney et al. 2014 and Catto et al. 2019), but few if any have examined the contribution made solely by the projected changes in SCE. While many observational and modelling studies have analyzed the effects of extensive distributions of snow cover on cyclone behavior (e.g. Namias 1962; Heim and Dewey 1984; Elguindi et al. 2005), few have explicitly studied the effects of reductions in snow cover besides Walland and Simmonds (1984), who did not experiment with projected SCEs. A few studies have suggested the importance of the snow extent boundary to cyclone behavior (e.g. Ross and Walsh 1986; Rydzik and Desai 2014), but there haven’t been modelling studies that experiment with shifts explicitly applied to these boundaries. While the pronounced effects of snow cover in the Great Plains has long been well understood and while regional modelling with snow forcing has been applied to the area (e.g. Elguindi et al. 2005), regional climate studies in the Great Plains focusing on projected snow extent retreat have not been performed. Finally, the greatest deficiency in all such regional, case-oriented modelling studies performed to date is the dependence on a small number of simulations. Examining several simulations across greater numbers of cases not only provides the statistical robustness of a large dataset but also the chance to examine the seasonality of the snow-cyclone relationship.
Snow retreat is of particular relevance given the likely changes in projected snow cover under anthropogenic climate change (Brown and Mote 2008; Peacock 2012; Notaro et al 2014; Krasting et al 2013). All studies point to reductions in North American persistent snow cover extent duration. However, there are also areas of increased snow cover and varying sensitivities depending on both temperature and precipitation trends. Simulations of climatological snow cover redistribution consistent with GCMs and studies of its impact on subsequent extratropical cyclones has not been done.
The purpose of this study was to determine whether changes in underlying snow cover on the Great Plains result in consistent, discernable influence on cyclone steering, intensity, and precipitation by conducting a broad survey of numerous cyclone simulations with snow cover perturbed up to 96 hours prior to cyclogenesis. Snow cover was perturbed with varying degrees of areal extent reductions and at multiple initialization times in order to determine if there is any spatial or temporal relationship between snow cover perturbation and changes in cyclone intensity or track. The analysis attempted to broadly define what direct effect, if any, North American snow cover reductions due to future climate change will have on extratropical cyclone events. This particular study did not intend to examine individual cases and outliers to explain dynamical relationships between the surface boundary conditions and the air aloft. An in-depth investigation of two simulations from this study is presented in Breeden et al. (submitted).
We hypothesize that, because cyclones preferentially track along the margin of snow extent (Ross and Walsh 1986; Rydzik and Desai 2014), cyclone trajectories in simulations with poleward-shifted snow lines will deviate poleward in kind. Because as much as 30% of the moisture in extratropical cyclones is obtained by surface evaporation (Trenberth 1998) and local precipitation recycling is significant in the Great Planes (Bagley et al. 2012), it is also expected that the removal of snow from the domain will result in appreciable increases in cyclone precipitation.
Methods
Experimental design and data
In order to test the effect of snow line position on extratropical cyclones, 20 cold season NA cyclones (Fig. 2) between 1986-2005 were simulated using the Advanced Research core of the NCAR Weather Research and Forecasting model (WRF-ARW) version 4.0.3 (Skamarock et al., 2019) with perturbed SCE. Four cyclone cases were subjectively selected from each of the months from November through March based on manual observational evaluation of all mid-latitude cyclones identified by low-pressure centers through this period in daily surface and upper-level weather charts. The criteria of selected cases required storm trajectories over or adjacent to the Great Plains study area which resemble either the Alberta Clipper track or that of the Colorado Low with lifetimes of at least 2 days, based on presence of well-defined central minimum pressure. Cases were chosen until a sufficient variety of differences in the lifetime minimum sea-level pressure (SLP) and magnitude of upper level forcings in the form of 500 hPa height curvature and vorticity advection by the thermal wind were found. Cases were simulated with observed initial conditions and validated against observations using the 32-km spatial resolution North American Regional Reanalysis (NARR; Mesinger et al. 2006) to ensure that WRF could accurately simulate each case.
Alterations to the SCE of each case were made by applying average poleward snow line retreat (PSLR) from the 20-year periods of 1986-2005 (historical) to 2080-2099 (projected) for each of the five months examined in this study. Projected PSLR was determined by examination of the grid cell snow mass change in 14 models of the 5thphase of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012) wherein daily snow mass data were available and experiments were conducted with two Representative Concentration Pathway forcings: RCP4.5 and RCP8.5 (van Vuuren et al. 2011)(Table 1). Grid cells were identified as snow-covered if their simulated snow mass was at least 5 kg m-2, which corresponds to typically 5 cm of snow depth (assuming a 10:1 snow to water ratio), sufficient to cover the surface. We did test other thresholds and did not find a strong sensitivity to this choice in the projected snow cover maps. The southernmost such grid cells were considered to comprise the snow line if the 5 degree span to the north of a cell had an average snow mass exceeding that threshold. This search radius was employed in order to exclude outlying isolated southern patches of snow. To limit artifacts that arise from small-scale variability in snow cover, a 600 km moving window average was then applied to all derived southern extent of snow cover, hereafter referred to as the “snow line”. For each month, the 20-year average snow line of the historical and projected periods was calculated, and the amount of projected PSLR was determined from west to east in 30 km-wide bins across North America. Different iterations, realizations, and physics options belonging to experiments of the same model were combined in a “one model, one vote” scheme. With PSLR calculated for both RCP forcings for each of the 14 models, each month contained 28 PSLR values from which the 10th, 50th, and 90th percentiles were determined (Table 2).
The modeling effort involved simulating each of the 20 mid-latitude cyclone cases with five degrees of snow line perturbation, each at five different initialization times, from zero to four days prior to cyclogenesis, yielding a total of 500 distinct simulations. One hundred simulations were generated without changes made to snow cover (control). The remaining 400 runs imposed projected snow line changes of varying magnitude (10th, 50th, and 90th percentiles; P10 ,P50 , P90 ) or complete snow removal across the domain in order to determine the degree to which the position of the snow line influences storms as opposed to that attributable solely to snow removal. Snow lines for perturbed simulations were determined by applying values of PSLR to corresponding 30 km bins of the snow lines, as determined based on the method above, for each case and removing all snow south of the new snow line except at altitudes greater than 2,000 m, where snowpack may persist even in warmer climates, based on conclusions by Rhoades et al. (2018). It should be noted that the removal of all snow south of the assigned snow line creates a discontinuous step function in snow depth, a hard margin which is not necessarily characteristic of real snow extent boundaries.
WRF model configuration
WRF-ARW simulations were executed in a domain comprising the continental United States (CONUS), central and southern Canada, northern Mexico, and much of the surrounding oceans. The WRF-ARW has previously been shown to be reliable in simulating seasonal temperature and precipitation dynamics over the United States (Wang and Katamarthi 2014), with biases in line with other mesoscale numerical weather models (Mearns et al 2012). We ran WRF-ARW with 30 km horizontal resolution to best capture synoptic scale transport, a 150 km buffer zone on each side, and 45 vertical levels (Fig. 3). Initial and lateral boundary conditions were derived from 3-hour NARR data provided in grib format by NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, at https://www.esrl.noaa.gov/psd/. Version 4.0 of WRF offers a “CONUS” suite of physics options which was used in this experiment, and appears to accurately reproduce large-scale circulations (Hu et al. 2018). The NOAH Land Surface Model (Noah LSM; Mitchell et al. 2001) was altered to reduce surface snow accumulation to zero during simulation in order to avoid snow deposition prior to the arrival of the cyclone of interest into the area without removing precipitation in the atmosphere. The Noah LSM uses a single layer snow model and calculates snow albedo according to the method developed by Livneh et al. (2010), which calculates the albedo of the snow-covered portion of a grid cell as
\begin{equation} \alpha_{\text{snow}}=\ \alpha_{\max}A^{t^{B}}\nonumber \\ \end{equation}
where αmax is the maximum albedo for fresh snow in the given grid cell (established by data from Robinson and Kukla, 1985), t is the age of the snow in days, and A andB are coefficients which are, respectively, 0.94 and 0.58 (0.82 and 0.46) during periods of accumulation (ablation). CoefficientsA and B were set to accumulation phase for simulations in every month except for March, when the snow was considered to be ablating. Sensible and latent heat fluxes are calculated from surface and snow using an energy balance approach based on snow pack temperature and moisture input.
Analytical methods
A number of tracking methods have been proposed for cyclones, as reviewed in Rydzik and Desai (2014). Here, cyclones are tracked by defining the center as the local SLP minimum and following it as the cyclone proceeds. Because each cyclone case was known in advance from subjective selection, identifying the genesis of each cyclone involved searching a known area at a known time for SLP minima. Recording changes in storm trajectory between two simulations of the same cyclone case is done by calculating the mean trajectory deviation (MTD), which is the sum of the absolute north-south deviation distance between the two storm centers (perturbed-control) at each corresponding time step divided by the number of time steps. Because each model time step is 3 hours, MTD is expressed in km (3h)-1.
Examination of precipitation amount and type involved isolating storm associated precipitation using the method introduced by Hawcroft et al. (2012). For each time step, it is assumed that precipitation attributable to any cold season cyclone simulation occurs within a 12° radial cap of the storm center. Analyzing the precipitation quantity of a cyclone’s lifetime required determining precipitation amounts and types from within the radial cap at each time step and ignoring those values outside of it.
To study broad changes in wind speed, we determined the integrated kinetic energy (IKE) of each simulated cyclone using a variant of the method first proposed by Powell and Reinhold (2007). IKE is determined by integration of the KE in the volume (V ) of the bottom model layer based on wind speeds (U ) and assuming a constant air density (ρ ) of 1 kg m-3,
\begin{equation} IKE=\ \int_{V}{\frac{1}{2}\rho U^{2}\text{dV}}\nonumber \\ \end{equation}
As with storm associated precipitation, IKE is only calculated within a 12° radius of the storms’ pressure minima. ΔIKE represents the normalized ratio of control to corresponding perturbed simulations.
Results
Snow cover trends
Before snow extent changes could be applied to model initialization data for perturbation experiments, it was necessary to conduct a survey of the 14 selected CMIP5 models to determine mean PSLR from the 1986-2005 period to 2080-2099 in the span east of the Rocky Mountains and west the Atlantic coast of North America (105 West to 55 West). The mean PSLR of both RCP experiments for each of the models is shown for each of the cold season months in Figure 4. The results show snow retreat differences among models is large although some trends are clear. All models for both experiments in all months show a projected poleward shift in snow cover extent, with a minimum average retreat in January of 51 km and a maximum in November of 1,025 km. The models show that the shoulder months of November, December and March experience greater PSLR than those in the middle of winter.
Generally, simulations of the RCP8.5 experiment yielded greater PSLR than RCP4.5, although all months have multiple exceptions. According to the Student’s t-test, RCP8.5 PSLR is greater than RCP4.5 within the 99.9% confidence interval except in March, when the confidence level shrinks to 99%. February has the lowest standard deviation of PSLR across models at 175 km, which is comparable to January and March with 185 km and 183 km, respectively. The early winter months have the higher standard deviations at 219 km and 210 km for December and November, respectively.
Cyclone trajectory
All control cases were selected to those where the control run well depicted observed cyclone trajectory and net precipitation. The 400 perturbation cases were then compared to these 100 control runs. Cyclone shifts in response to imposed snow cover extent (SCE) shifts, expressed as mean trajectory deviation (MTD), were quite small, often less than the domain grid spacing of 30 km (55% of the time), and only infrequently did they exceed two entire grid spaces (12%), indicating that cyclones in perturbed simulations followed their control counterparts faithfully with only minor exceptions. The different in cyclone track in the perturbed snow cover cases relative to the control was modest, usually smaller than the mean difference of the control case to observed cyclone tracks in reanalysis and not related to cyclone typo (e.g., Alberta clipper and Oklahoma panhandle cyclone). Although the simulated responses of these cyclones to short-term reductions in snow extent may be regarded as minute, they are not always devoid of significance.
Plotted together according to total area of snow removed (Fig. 5a), the 400 MTDs of each perturbed simulation cyclone present a mild but significant positive linear relationship (R2 = 0.232, p < 0.01). The strength of the relationship increases when limited to simulations initialized at least two days out (R2 = 0.292, p < 0.01), which implies an adjustment timescale for cyclone dynamics to respond to SCE changes. Perturbed simulations initialized at the time of cyclogenesis or one day prior have diminished MTDs when compared to the magnitude of the responses for simulations initialized two days out and greater where the signal stabilizes, with gradual increases in the mean at three and four days prior (Fig. 5b). Figure 5c also reveals this relationship by calendar month, thereby summarizing the MTD response according to initialization time as well as perturbation degree. There was some seasonality to the results with weaker responses in the late autumn to early winter (Nov-Dec, average MTD (μ)=29.6 km 3 hr-1), the strongest responses in mid-winter (Jan-Feb, μ=36.1), and moderate responses in late winter to spring (Mar, μ=32.7), implying albedo was not the dominant mechanism driving changes.
Despite the prevalence of cyclone trajectory deviation among perturbed simulations, except for a few outliers, poleward deflection across cases in response to a retreating snow line was not substantial (Fig. 5d). The mean values and quartiles for each month never exceeded a single grid space, even when only initialization times greater than 2 days out are considered. Mean poleward deviation was nearly as likely to be negative as positive for most cases. Another way to examine cyclone steering is by calculating the tendency of trajectories to deviate toward the perturbed snow line, which may lie to the south of the cyclone trajectory. This method, however, also falls short of producing a robust signal. Like the poleward shift, the snow line oriented shift only indicated a positive signal a little over half the time and rarely with substantial quantities.
Storm center SLP
Across all simulations, 70% of perturbed simulation cyclones decreased average lifetime central low SLP compared to the corresponding control simulation, however slightly, and every perturbed simulation cyclone experienced a significant difference in central SLP compared to control at some point in their lifetime. Most central SLP differences present in perturbed simulations, like those in the analysis of the MTDs, are small. The magnitude of mean cyclone lifetime central SLP change never exceeded 2.2 hPa, though maximum differences could exceed 10 hPa. Figure 6a summarizes lifetime mean central SLP changes averaged for each month according to perturbation degree and initialization time. There is a robust dependence on perturbation degree for the latter three months of the cold season, with exceptional responses in the no snow simulations. November and December cyclones have a much weaker, if at all discernable, response to the degree of PSLR. For all perturbed simulations, November and December cyclones only undergo a mean lifetime deepening 53% of the time, while mean lifetime deepening occurs 81% of the time for the latter three months. In these latter months, responses become more pronounced when simulations are initialized 2-4 days prior to cyclogenesis, although responses within this period are similar. Cyclones in transit over the region where snow had been removed deepened, on average, 2.5 times as much as others and nearly 4 times as much as those which remained over snow (p < 0.01 ).
The maximum instance of deepening for cyclones in perturbed simulations decreased almost 6 hPa (Fig. 6b), although the deepening was more strongly responsive to initialization time and also tended to stabilize at greater than two days out. The dependence on month is not as great for percentile-based snow removals, although it is stark in the no snow simulations with the same latter months have a much greater response. Maximum deepenings are more robustly correlated with MTD (Fig. 6c,R2 = 0.395, p < 0.01). The relationship is notably less robust when examining mean lifetime pressure change (R2 = 0.209, p< 0.01), although still statistically significant.
Kinetic energy
Across all 400 perturbed simulations, 72% of perturbed simulations experienced a positive mean intensification over their lifetime compared to the control runs. Figure 7a shows that, like the other previously examined variables, changes in integrated kinetic energy (IKE) relative to control caused by perturbation of the snow fields in the vicinity of cyclones is highly subject to both initialization time and perturbation degree. One notable difference regarding IKE is that there is an exceptional tendency for it to abate at the higher perturbation degrees, particularly in the shoulder months. At the 90thpercentile of snow cover reduction, almost all November storms experience a mean reduction in IKE, and November storms undergo an average 1% decrease in IKE when initialized 4 days prior to cyclogenesis. The fact that those very same cyclones experience a nearly equivalent increase in IKE when initialized three days out indicates a nonlinear relationship between IKE and initialization time. For no snow simulations, short spin-up times generally reduced IKE, although simulations with one day of spin-up or greater increased IKE substantially. The December cyclones appear to be the only exception to these observations.
The seasonality of changes to IKE are made more apparent in Figure 7b, which shows that, on average, December simulations experienced a small decrease in intensity, although some outliers decreased intensity by over 3%. Simulations in every other month intensified on average, though the signal was considerably weaker in November than in the mid-winter and spring months. Some outliers in February and March intensified by over 9% when all snow was removed, and those months still had dramatic intensifications in the 50th and 90th percentile experiments.
Precipitation
Among perturbed simulations, 86% of cases experiences an increase in domain-integrated precipitation (Fig. 8a). Of the variables examined in this study, precipitation had the strongest response to removed snow cover and the greatest sensitivity to initialization time. Precipitation in perturbed simulations had weak responses when no spin-up time was allowed except in no snow simulations. Once again, December cases had the weakest responses to the snow cover perturbations with the lowest mean change in domain-integrated precipitation (Fig. 8b). While January cases had the highest mean response in precipitation, the November cases had the highest total increases in precipitation within individual cases and the greatest spread among cases. In many perturbed simulations, the phase of the precipitation changed from snow to rain, in southern latitudes and often near the original snow line. While grid cells with such phase changes never exceed 2% of the cells in the study domain, the overall increase in precipitation across the domain contributed to a substantial generation of new rain in perturbed simulations.
Changes in the volume of precipitation were very regionally dependent. In response to a poleward retreating snow line, cyclone-associated precipitation increased substantially across regions where snow was removed and across northern latitude regions downstream of the Great Plains, while southern regions experienced decreases in total precipitation (Fig. 9). The locations and amounts of enhanced precipitation appear to have been largely dependent on whether snow had been removed in that area, but new precipitation was often generated over snow near the perturbed snow line. Removing all snow from the domain resulted in significant quantities of new precipitation, particularly in the latitudes north of the U.S.-Canada border with the province of Quebec receiving an average of 1 mm of extra precipitation per grid cell and the Southeast United States (Virginia, North Carolina, South Carolina, Georgia, Florida, Alabama) experiencing an average decrease of 0.05 mm per grid cell.
Discussion
The retreat of southern snow extent calculated by comparing averages of historical (1986-2005) and late twenty-first century (2080-2099) snow lines is substantial. Surprisingly, applying it to historical cyclone cases for simulations with spin-up times of 4 days or less fails to result in striking changes to cyclone trajectory or central minimum SLP, notwithstanding the conclusions of other studies which would suggest a more direct and influential relationship. The changes made to underlying snow cover did, however, produce noteworthy responses to cyclones’ total kinetic energy and the storm-associated precipitation within a broad radius of the storm center.
Storm-associated precipitation had the most robust positive relationship to snow removal with the highest percentage of perturbed simulations yielding greater amounts of either solid or liquid precipitation. Even simulations with decreased domain-wide precipitation had relatively little reduction relative to increase. This outcome agrees with a large number of previous works which find an increase in precipitation amount and intensity in the Northern Hemisphere by the late 21st century (e.g. Catto et al. 2019). This study, however, does not find this result due solely to the Clausius-Clapeyron relationship whereby a warming climate drives increases in airborne water vapor but primarily due to the removal of snow from the surface and its effect on atmospheric thermodynamics. This finding supports observations made by previous authors (e.g. Namias 1962 and Elguindi et al. 2005) that snow cover suppresses precipitation from overhead extratropical cyclones. This may be due to the lack of moisture flux (Trenberth 1998 and Ellis and Leathers 1999), the increase of static stability (Bengtsson 1980), or more likely both. We can therefore reasonably assume that the increases in precipitation shown here represent only a portion of the increased precipitation for which climate change will be responsible and that the poleward migration shown is likely to be more intense.
The cyclone integrated kinetic energy (IKE), a measure of 10 m wind speed associated with the storm, also had noteworthy responses to snow removal. A large majority of cyclones in perturbed simulations intensified and there exists a positive relationship between IKE and snow removal area, suggesting that it may be related to surface energy budget. These results contradict those of GCM studies such as Ulbrich et al. (2009) and Seiler and Zwiers (2015) which find reductions in extratropical cyclone wind speed by the late 21stcentury. These results may differ due to changes in upper level baroclinicity caused by climate change or even due to the differences in our determination of intensity.
In contrast to our original expectations, trajectory deviations were minimal. MTD measures the amount of deviation from control in perturbed simulation cyclone trajectories and averages over the cyclones’ lifetimes. Because the majority of cyclones in perturbed simulations did not deviate from control for most of their courses, most MTDs are measured as less than the length of the domain grid spacing of 30 km and only 49 of the 400 perturbed cyclone simulations deviated by an average of more than two grid spaces. Directional MTDs considering deflection toward the North Pole or the perturbed snow line were inconsistent and notably minimal with few outliers. The fact that both of these metrics often yielded such minor results indicates that many perturbed simulation cyclones deviated primarily stochastically from the control trajectory.
The study by Elguindi et al. (2005) wherein snow was added to a Great Plains nested domain two days prior to cyclone arrival generated similar trajectory outcomes with deviations in perturbed cases only rarely exceeding 100 km. The trajectory deviations in these tests, like our own, varied substantially and defied any discernable trend. It is reasonable to infer that differences in trajectories between control and perturbed cyclones in both studies are likely chaotic reactions to considerable energy disturbances caused by step changes to surface conditions over extensive areas, rather than functional responses to the specific positioning of snow cover. This conclusion is unexpected, given the significant cyclone responses to snow anomalies found by multiple observational studies (e.g. Dickson and Namias 1976; Heim and Dewey 1984; Rydzik and Desai 2014) as well as modelling done by Ross and Walsh (1986) and Walland and Simmonds (1997). The most obvious distinction here is temporal scale, hinting that a similar study conducted at the seasonal timescale may reveal a more robust relationship.
Like MTDs, changes to cyclones’ central low SLP due to a retreating snow line were minimal. This, however, differed from the results of Elguindi et al. (2005) who found an average positive difference of 4 hPa in response to expanded snow cover, a threshold which only seven simulations in this whole study exceeded, all of which as part of the no snow sensitivity experiment. Perhaps this can be attributed to the fact that they added snow as opposed to removing it or to the physics of the MM5 model compared to WRF-ARW. Even with the disparity in the magnitude of pressure changes, their discovered trend of central pressure increasing when snow is added is complemented by the findings of this study where snow removal generally contributed to a decrease in central pressure. The enhanced frequency of central low SLP decreasing while in transit over regions where snow had been removed corroborates the conclusion of Elguindi et al. (2005) that snow cover prevents the deepening of mid-latitude depressions by reducing warm sector temperature and moisture gradients, weakening surface convergence and fronts. The relationship shown between MTD and pressure changes indicates that, while the two may not be directly linked, they do respond similarly to perturbed simulations.
Seasonality
We find consistent trends across all examined variables affirming that cyclone responses to poleward-shifted snow lines depend upon when in the cold season the cyclones occur. Generally, responses of virtually every investigated variable are greater in the mid-season months of January and February and weaker in the shoulder months, although there are often greater responses in March than in November or December. This is counter-intuitive from an inspection of PSLR as seen in Figure 4. If anything, there appears to be an inverse relationship between amount of mean PSLR and response of cyclones to the correspondingly-shifted snow lines. However, it has been shown that the surface temperature effect of snow cover is strongest in late winter (Walsh et al. 1982) and March snow cover has reduced efficacy due to its properties during ablation (Livneh et al. 2010).
December consistently has the most abnormal responses, even contradicting consistent trends in other months; for example, December is the only month with a mean reduction in IKE but is also the month with the weakest solar radiation, implying a weaker albedo gradient effect. Still, it is not entirely understood why December in particular has these properties, although the apparent problem in attempting to make determinations about the seasonality of the data is that each month represents only four separate cyclone cases which are then averaged together.
Conclusions
Twenty cold season extratropical cyclones over or near the North American Great Plains were generated in a series of simulations in order to gauge the dependence of their trajectories, intensities, and associated precipitation on underlying snow cover. When a realistic retreat of snow cover consistent with climate warming scenarios was applied to these cases, a majority of cyclones experienced an average decrease in pressure and increases in precipitation, but only limited changes in trajectory and modest increases in kinetic energy. These results contradict expectations gained from observational studies such as that of Namias (1962) and Rydzik and Desai (2014) but reflects the results of modelling done by Elguindi et al. (2005), reflecting a continued disagreement among models and observations.
It is yet unknown why the cyclone trajectories did not adhere more closely to shifted snow lines, as the findings of other studies would have suggested. Weaker responses to the removal of snow cover at the time of cyclogenesis suggest that the presence or absence of the snow margin has a minor, though not entirely imperceptible, immediate effect. There is little to imply that the effect on trajectory deviation, pressure change, or precipitation plateaus for simulations initialized at four or more days out and so the question of the full extent of the snow margin’s influence cannot be answered until longer case study simulations are executed.
Lingering questions remain on mechanisms of snow cover on sea-level pressure, differences among cases in surface energy-balance and radiative properties and its influence on cyclone dynamics, and upper-level dynamics. Some of these, especially upper-level dynamics, are studied in individual cases in detail in a companion paper by Breeden et al. (submitted). The simulation model outputs provide a rich data set for future evaluation and a provided at the archive below for public access.