Extratropical Cyclone Response to Projected Reductions in

: Extratropical cyclones develop in regions of enhanced baroclinicity and progress along


40
Northern Hemisphere snow cover is, at its seasonal maximum, the largest compo-41 nent of the terrestrial cryosphere and exerts considerable influence on the mid-latitude 42 atmospheric circulation through a diverse set of mechanisms [1,2]. Snow cover depresses 43 near-surface air temperature due to increased albedo at day and greater radiational cool- 44 ing at night [3]. Its properties lead to an effective sink of sensible and latent heat [4], 45 contributing to an increase in static stability [5,6] and a reduction of moisture flux into the 46 atmosphere [7]. This inhibition of upward moisture flux may be responsible for the nega- 47 tive correlation between snow cover and precipitation observations [8] and models [9,10]. 48 Studies have also shown that the total extent of continental snow cover is sometimes 49 responsible for modulating upper-level circulation [9][10][11][12][13][14][15] and that accurately initializing 50 snow cover can improve subseasonal forecast skill considerably [16,17]. It is because of 51 this apparent relationship between snow cover and atmospheric circulation that determi-52 nation of the regional dependence and the temporal scales at which snow cover drives 53 responses in the atmosphere is of fundamental importance to both short-and long-term 54 forecasting. 55 Observations and hypotheses about the influence of established snow cover extent 56 on the characteristics of ensuing synoptic weather systems may have begun with Lamb 57 [18]. However, one of the first analyses of this relationship was provided by Namias [11], 58 who hypothesized that the abnormally extensive North American snow cover of the win-59 ter of 1960 had contributed to the more frequent and intense cyclone development ob-60 served along the Atlantic coast by enhancing baroclinicity between the continent and the 61 much warmer ocean. Dickson and Namias [19] subsequently showed that periods of great 62 continental warmth or cold in the American Southeast had a direct influence on the 63 strength of the baroclinic zone near the coast and would affect the average frequency and 64 positions of extratropical cyclones, drawing them further south when the region was 65 colder. Likewise, Heim and Dewey [20] showed that extensive North American snow 66 cover contributed to a greater frequency of cyclones in the southern Great Plains and 67 Southeast and a reduction in the frequency of cyclones tracking further north due to dis-68 placement of storm track. From 1979-2010 in North America, cold season mid-latitude cy-69 clones were more frequently observed in a region 50-350 km south of the southern snow 70 extent boundary (snow line) [21]. That study also noted a similar distribution of low-level 71 baroclinicity around the snow line. 72 Modeling studies have indicated a similar relationship between snow extent and ex-73 tratropical cyclone statistics. Ross and Walsh [22] studied the influence of the snow line 74 on 100 observed North American cyclone cases which progressed approximately parallel 75 to the baroclinic zone within 500-600 km of the snow line. By measuring forecast error 76 from a barotropic model, they determined that baroclinicity associated with the snow 77 boundary was an important factor in cyclone steering and intensity. Wallace and Sim-78 monds [9] performed global climate model (GCM) experiments with forced anomalously 79 high and low extents of realistic snow cover distributions, ultimately finding a reduction 80 in North American cyclone frequency when snow cover was more extensive with cyclones 81 frequently occurring further south, similar to the observations of Heim and Dewey [20], 82 owing to changes in baroclinicity induced by meridional temperature gradients. Alexan-83 der et al. [14] conducted GCM studies demonstrating that reduced snow cover led to 84 greater absorbed solar radiation and increased latent, sensible, and longwave fluxes from 85 the surface, leading to continental scale warming, especially in fall and spring. 86 The most extensive study to date is the one of Elguindi et al. [10], who used a 25-km-87 resolution nested domain over a portion of the Great Plains in the Penn State-National 88 Center for Atmospheric Research (NCAR) Mesoscale Model (MM5) and simulated eight 89 well-developed cyclone cases with snow cover added throughout the domain, initializing 90 48 hours prior to each cyclone's arrival to the inner domain. All perturbed cyclone case 91 simulations underwent an increase in central pressure and decrease in total precipitation 92 with slight shifts in the cyclone trajectory which were highly variable and inconsistent but 93 related to reduction in surface temperature and moisture gradients at the surface and 94 across fronts. However, this study was limited to one control and one perturbed run per 95 case. Further, the perturbed cases are an extreme scenario of adding snow to the entirety 96 of the inner domain rather than altering the position of continental snow cover extent 97 which we hypothesize here to have a larger impact on storm dynamics through changes 98 low-level baroclinicity in the storm track region. Here, we seek to address these 99 deficiencies in prior studies with a focus on realistic snow cover retreat experiments across 100 multiple perturbations and a larger number of cases.

101
In North America, the net effects of snow cover are nowhere more pronounced than 102 the Great Plains region, which has the highest local maximum of snow albedo [23,24] and 103 where the strongest correlation between snow cover and negative temperature anomalies 104 has been observed [20,25]. In the Great Plains region exists one of the largest disparities 105 between local maximum snow albedo and background land surface albedo on the conti-106 nent, suggesting the greatest albedo gradient across a snow line (Figure 1). The land sur-107 face is characterized by high inter-and intra-annual snow cover variability [26] and low 108 surface roughness. Winter cyclones track over the Great Plains with high frequency due 109 in part to areas of enhanced cyclogenesis in the lee of the Rocky Mountains [27][28][29]. The 110 two most prolific cyclogenetic zones over the North American landmass account for the 111 two types of cyclone tracks studied here: the Alberta Clipper track, which typically begins 112 in Alberta, Canada and proceeds to the southeast [30], and the Colorado Low track, which 113 starts near southeast Colorado and often proceeds northeastward towards the Great Lakes 114 region [28]. Because of their spatial extent and frequency in the region, extratropical cy-115 clones contribute substantially to the hydrology of the Great Plains, accounting for greater 116 than 80% of the total winter (December through February) precipitation throughout much 117 of the region [31].  Robinson and 120 Kukla, 1985) and background surface albedo calculated by the WRF Preprocessing System with in-121 put from the WRF vegetation parameter lookup table. In principle, these values represent the max-122 imum albedo gradient across a hypothetical local snow line. The large region of enhanced albedo 123 difference in the center of the continent represents the Great Plains study area. While large snow 124 albedo gradients are also found in the desert Southwest and Mexico, it is rare for snow to persist in 125 most locations there.

126
Typically, simulations of projected future climate states are implemented with global 127 climate models, which are limited by expansive resolutions over a global domain and 128 large timesteps which do not allow the models to simulate phenomena such as cyclone-129 associated precipitation and its diurnal cycle as accurately as regional weather models. 130 Harding et al. [32] demonstrated that dynamically downscaling Coupled Model Intercom-131 parison Project Phase Five (CMIP5) simulations to 30 km resolution in the NCAR WRF 132 model improved the simulation of precipitation, especially extreme precipitation events, 133 in the Central U.S. Many modeling studies have applied global and regional climate 134 models to study the projected behavior of extratropical cyclones in the late 21 st century 135 [33,34] but few if any have examined the contribution made solely by the projected 136 changes in snow cover extent. While the pronounced effects of snow cover in the Great 137 Plains have long been well understood and while regional modelling with snow forcing 138 has been applied to the area [10], regional climate studies in the Great Plains focusing on 139 projected snow extent retreat have not been performed.

140
Snow retreat is of relevance given the ongoing and projected changes in snow cover 141 under anthropogenic climate change including over North America [35,36,37,38,39]. All 142 studies point to reductions in North American persistent snow cover extent duration. 143 However, there are also areas of increased snow cover and varying sensitivities depend-144 ing on both temperature and precipitation trends. Here, we seek to determine whether 145 changes in underlying snow cover across the Great Plains result in a consistent, discerna-146 ble influence on cyclone steering, intensity, and precipitation by conducting a broad sur-147 vey of cyclone simulations with snow cover perturbed 4 days prior to cyclogenesis. Unlike 148 prior studies, here snow cover was perturbed with multiple degrees of areal extent reduc-149 tions to determine if there is any spatial or temporal relationship between snow cover 150 perturbation and changes in cyclone intensity or track. The analysis attempted to broadly 151 define what direct effect, if any, North American snow cover reductions due to future 152 climate change will have on extratropical cyclone events. A greater in-depth investigation 153 of mechanisms within two simulations from this study is presented in Breeden et al. [6]. 154 We hypothesize that, because cyclones preferentially track along the margin of snow 155 extent [21,22], cyclone trajectories in simulations with poleward-shifted snow lines will 156 deviate poleward in kind. Because a significant proportion of moisture in extratropical 157 cyclones can be obtained by surface evaporation, even in winter [40] and local precipita-158 tion recycling has been shown to be significant in the Great Plains [41], it is also expected 159 that the removal of snow from the domain will alter hydrological cycling [42], leading to 160 increases in cyclone-associated precipitation.

186
These 20 cases were then simulated with observed initial conditions and validated 187 against sea-level pressure observations using the 32-km spatial resolution North Ameri-188 can Regional Reanalysis (NARR) [44] to ensure that WRF could accurately simulate tra-189 jectory and intensity of each control case within a reasonable tolerance (trajectory of low 190 pressure center within ~200 km through interior of Great Plains model domain). Of the 191 original 20 selected cases, this qualitative validation criteria led us to select 15 of those 192 cases (Table 1) (Table 204 2). Grid cells were identified as snow-covered if their simulated snow mass was at least 5 205 kg m -2 , which corresponds to typically 5 cm of snow depth (assuming a 10:1 snow to water 206 ratio), sufficient to cover the surface. We did test other thresholds and did not find a strong 207 sensitivity to this choice in the projected snow cover maps. The southernmost such grid 208 cells were considered to comprise the snow line if the 5-degree span to the north of a cell 209 had an average snow mass exceeding that threshold. This threshold was employed in or-210 der to exclude outlying southern patches of snow. To limit artifacts that arise from small-211 scale variability in snow cover, a 600 km moving window average was then applied to all 212 derived southern extent of snow cover, hereafter referred to as the "snow line". Using these criteria, then for each month, the 20-year average snow line of the histor-216 ical and projected periods was calculated, and the amount of projected snow line retreat 217 was calculated from west to east in discreet bins across North America. As cases take place 218 on a 30 km horizontal resolution model grid, each 30 km zonal coordinate was assigned 219 values of poleward snow retreat from the averaged, derived values for each month and 220 perturbation condition. Since each GCM experiment has differing numbers of realizations, 221 initializations, and physics options, we combined these in a "one model, one vote" scheme 222 for calculating snow retreat. With snow line retreat calculated for both RCP forcings for 223 each of the 14 models, each month then contained 28 poleward snow line retreat values 224 from which the 10 th , 50 th , and 90 th percentiles were extracted (  Table 3. Models and RCP experiments which were used to determine 10 th , 50 th , and 90 th percentile 227 values for late twentieth to late twenty-first century snow line retreat in eastern North America. The modeling effort involved simulating each of the 15 validated mid-latitude cy-230 clone cases with a control case (observed snow cover as initialized in North American 231 Regional Reanalysis boundary condition) and three snow line perturbation boundary con-232 ditions (10 th , 50 th , and 90 th percentile, P10, P50, P90), at four days prior to cyclogenesis, where 233 liquid snow cover on ground was reduced to zero. We also conducted simulations with 234 initialization times from 0-3 days prior to cyclogenesis to evaluate consistency and inter-235 nal variability, but here we focus on the results of four days prior to allow for atmospheric 236 adjustment while still allowing for examples where the model realistically captures cy-237 clone trajectory in the control cases (Supplemental Fig. S4). An additional set of simula-238 tions with all snow removed was also performed as a test response case and archived in 239 the model output repository, but not directly analyzed here (except in the tests of initiali-240 zation times).

241
Control simulations were run for all 15 cases. The remaining runs imposed the three 242 projected snow line changes to determine the degree to which the position of the snow 243 line alone influences cyclone behavior. Snow lines for perturbed simulations were deter-244 mined by applying values of snow line retreat to corresponding 30 km bins of the snow 245 lines, as determined by the method above, for each case and removing all snow south of 246 the new snow line except at altitudes greater than 2,000 m, where snowpack may persist 247 even in warmer climates, as concluded by Rhoades et al. [47]. It should be noted that the 248 removal of all snow south of the assigned snow line creates a discontinuous step change 249 in snow depth, a hard margin which is not necessarily characteristic of real snow extent 250 boundaries. . We ran WRF-ARW 257 with 30 km horizontal resolution, a 150 km buffer zone on each side, and 45 vertical levels 258 (Fig. 3). Initial and lateral boundary conditions were derived from 3-hour NARR data pro-259 vided by NOAA at https://www.esrl.noaa.gov/psd/. Version 4.0 of WRF offers a "CO-260 NUS" suite of physics options which was used in this experiment and appears to accu-261 rately reproduce large-scale circulations [50]. The Noah Land Surface Model (Noah LSM) 262 [51] was altered to reduce surface snow accumulation to zero during simulations in order 263 to avoid new snow accumulation prior to the arrival of the cyclone of interest into the area 264 without removing precipitation in the atmosphere. The Noah LSM uses a single layer 265 snow model and calculates snow albedo according to the method developed by Livneh et 266 al. [52], which calculates the albedo of the snow-covered portion of a grid cell as  [52]. 274 We did not incorporate ensemble simulations as might be needed in global general 275 circulation model experiments (e.g., [53]), as natural variability in a regional climate 276 model is partly constrained by nudging at lateral boundaries. We did not apply any spec-277 tral nudging to the large-scale fields aloft beyond the nudging along the lateral boundary 278 conditions. This approach allows the model to fully generate atmospheric feedback re-279 sponses to modified snow experiments without being nudged back to reanalysis 280 state. Our analysis of cases with earlier initialization times show consistency in key statis-281 tics (e.g. cyclone trajectory deviation and maximum pressure change) after 2 days of ini-282 tialization (Supplemental Figure S4), so we focus our results here on 4 days of initializa-283 tion. We used a single physics parameterization (WRF CONUS convection-permitting 284 physics suite, https://www2.mmm.ucar.edu/wrf/users/physics/ncar_convec-285 tion_suite.php) that reasonably replicated the observed storm trajectory in each of the 286 control case. These are the same physics parameters as the National Weather Service 287 (NWS) uses to create their United States national forecasts. 288

289
Several tracking methods have been proposed for cyclones, as reviewed in Rydzik 290 and Desai [21]. Here, cyclones are tracked by defining the center as the local SLP minimum 291 and following it as the cyclone proceeds. Recording changes in storm trajectory between 292 two simulations of the same cyclone case is done by calculating the mean meridional tra-293 jectory deviation (MTD), which is the sum of the absolute north-south deviation distance 294 between the two storm centers (perturbed-control) at each corresponding time step di-295 vided by the number of time steps. Because each model time step is 3 hours, MTD is ex-296 pressed in km 3h -1 .
As with storm-associated precipitation, integrated KE is only calculated within a 12° ra-311 dius of the storms' pressure minima. ΔIKE then represents the normalized ratio of control 312 to corresponding perturbed simulations.    Generally, simulations of the RCP8.5 experiment yielded significantly greater snow 334 line retreat than RCP4.5 (p<0.01). February has the lowest standard deviation of snow line 335 retreat across models over both RCP experiments at 175 km, which is comparable to Jan-336 uary and March with 185 km and 183 km, respectively. The early winter months have the 337 higher across-model standard deviations at 219 km and 210 km for December and No-338 vember, respectively, implying less consistent agreement among models for snow line re-339 treat in autumn over mid-winter.

Cyclone trajectory 341
As noted, in 15 of the 20 cases, the control run faithfully reproduced the observed 342 cyclone trajectory, generally true for the more well-defined cyclones. The perturbation 343 cases were then compared to these control runs (Figures 4-6). Cyclone track shifts in re-344 sponse to imposed snow cover extent shifts, expressed as mean trajectory deviation 345 (MTD), were quite small, often less than the domain grid spacing of 30 km (71% of the 346 time), and only infrequently did they exceed two entire grid spaces, indicating that cy-347 clones in the perturbed simulations followed their control counterparts faithfully with 348 only minor exceptions (Figures 7). MTD reflects an average deviation, but as the figures 349 show, the actual track at any one time step may deviate much more. The difference in 350 cyclone track was also not related to cyclone type (e.g., Alberta clipper and Colorado low 351 cyclone). Although the simulated responses of these cyclones to reductions in snow extent 352 may be regarded as small, they are not always devoid of significance. The largest trajec-353 tory deviations occurred in the late season clipper system (6 March 1987) case with 50 th 354 percentile snow removal (Figure 6b). Deepest pressure changes occurred also in a late 355 season Rocky Mountain low (22 March 1994) but with 90 th percentile snow removal (Fig-356  ure 6f). This case had the greatest absolute pressure deepening among all other cases for 357 all levels of snow removal (Figure 7).   Table S1). Deviations from the snow line were nearly equally dis-380 tributed poleward and equatorward.

Cyclone characteristics 386
Across all simulations, 82% of perturbed simulation cyclones decreased in average 387 central SLP compared to the corresponding control simulation, however slightly, and 388 every perturbed simulation cyclone experienced a significant difference in central SLP 389 compared to the control at some point in their lifetime (Supplemental Table S2). Most 390 central SLP differences present in perturbed simulations, like those in the analysis of the 391 MTDs, are small. The magnitude of the mean change in cyclone lifetime central SLP aver-392 aged -0.21 hPa and average simulation maximum pressure change was -0.48 hPa with 393 only two cases where pressure deepened by at least 2 hPa above the control (Figure 4b). 394 There is a significant linear relationship of decrease in minimum SLP on snow removal 395 extent (r 2 =0.17, p<0.01). No significant difference in slope is found by season, but the larg-396 est pressure changes are found in late winter. Cyclones in transit over the region where 397 snow had been removed deepened, on average, 2.5 times as much as others and nearly 4 398 times as much as those which remained over snow (p<0.01). Maximum pressure decrease 399 are also correlated with MTD (r 2 =0.395, p<0.01). 400 Most perturbed simulations experienced a positive intensification of integrated KE 401 over their lifetime compared to control runs. One exception is for integrated KE to de-402 crease at the higher snow removal amounts, particularly in early season. At the 90 th per-403 centile of snow cover reduction, almost all early season storms experience a mean reduc-404 tion in integrated KE of 1-3%. Simulations in every other season and snow removal sce-405 nario intensified by a similar amount.

406
Among perturbed simulations, 86% of cases experience an increase in cyclone-asso-407 ciated precipitation (Figure 9a). Mid-season cases had the weakest responses to the snow 408 cover perturbations with the lowest mean change in domain-integrated precipitation. 409 However, there is no significant relationship between snow removal extent and precipi-410 tation change. In many perturbed simulations, the phase of the precipitation changed 411 from snow to rain in up to 2% of grid cells, in southern latitudes and often near the original 412 snow line, linked to increase in air temperature, but again no clear relationship with snow 413 removal existed among the experiments (Figure 9b). More of the overall increase in pre-414 cipitation across the domain fell as rain than snow. 415 Figure 9. Relationship of snow removal to a) total storm-associated precipitation and b) fraction of 418 frozen precipitation, colored by season. Linear relationships were not significant. 419 Changes in the volume of precipitation were regionally dependent. In response to a 420 poleward retreating snow line, cyclone-associated precipitation increased substantially 421 across regions where snow was removed and across northern latitude regions down-422 stream of the Great Plains. Meanwhile, southern regions experienced decreases in total 423 precipitation, and intensity of overall precipitation changes increasing with greater snow 424 removal, even if the average change was unaffected (Figure 10). The locations and 425 amounts of enhanced precipitation appear to have been largely dependent on whether 426 snow had been removed in that area, but additional precipitation was also generated over 427 snow near the perturbed snow line and not as commonly over areas continuing to remain 428 snow covered. The "all snow removed" cases were used to quantify the overall effect of 429 snow removal on precipitation. In those cases, the average increase in precipitation per 430 grid cell experiencing an increase was 1 mm, while the average decrease in grid cells ex-431 periencing reduction was 0.05 mm.

443
The changes made to underlying snow cover did produce responses to cyclones' total 444 kinetic energy and the storm-associated precipitation within a broad radius of the storm 445 center. These effects are further explored in Breeden et al. [6], where two cases studied 446 here are diagnosed in further detail. A potential vorticity budget analysis reveals that 447 snow removal led to lower heights in the center of cyclones, which strengthened the cir-448 culation and potentially enhanced moisture advection, particularly from the Gulf of Mex-449 ico. However, the snow removal also led to an opposing effect in responses of stability 450 and relative vorticity, which muted the overall response of the mid-latitude cyclone 451 deepening and precipitation enhancement. These competing physical mechanisms at least 452 partly explain the modest responses of the examined extratropical cyclones in the current 453 paper to imposed snow line retreat.

454
Storm-associated precipitation had the most robust response to snow removal with 455 the highest percentage of perturbed simulations yielding greater amounts of either solid 456 or liquid precipitation. This outcome agrees with a large number of previous works which 457 find an increase in precipitation amount and intensity in the Northern Hemisphere by the 458 late 21 st century (e.g. [34]). This result is not a reflection of the Clausius-Clapeyron rela-459 tionship whereby a warming climate drives increases in evaporation and atmospheric wa-460 ter vapor, but rather it is due to the removal of snow from the surface and its effect on 461 lower atmosphere temperature. This finding supports observations made by previous au-462 thors [10,11] that snow cover has a negative relationship with precipitation from overhead 463 extratropical cyclones and agrees with climate model simulations noting increases in pre-464 cipitation extremes from extratropical cyclones as a result of thermodynamic responses to 465 increased diabatic heating [54,56,57]. The presence of snow cover has been linked to a 466 reduction of moisture flux [7,40] and an increase of static stability [5,6], thus its removal 467 promotes evaporation and instability, an effect consistent with our results. This effect is 468 also consistent with findings of the impact of the Great Lakes on surface turbulent flux 469 enhancement and SLP decreases in winter [58]. We can therefore presume that the in-470 creases in precipitation shown here represent only a portion of the increased precipitation 471 for which climate change will be responsible and that the poleward migration shown is 472 likely to be more intense.

473
The cyclone-integrated kinetic energy, a measure of boundary-layer wind speed as-474 sociated with the storm, also had noteworthy responses to snow removal. While the rela-475 tive magnitude is small, this represents a large total net increase of energy over the storm 476 lifetime. A large majority of cyclones in the perturbed simulations intensified, suggesting 477 that it may be related to changes in the surface energy budget that influence boundary 478 layer wind profiles. These results are a consequence of the snow removal, but the effect 479 may be mitigated in future climates by reduced baroclinicity arising from polar amplifi-480 cation [59,60], leading to overall reductions in extratropical cyclone wind speed by the late 481 21 st century. We hypothesize that the kinetic energy increases in the absence of cyclone 482 deepening in terms of central pressure is occurring because of increased anticyclonic de-483 velopment at the peripheries, although this hypothesis requires further exploration.

484
Our modeled trajectory deviations were typically smaller than those found in obser-485 vational studies. MTD measures the amount of deviation from control in perturbed sim-486 ulation cyclone trajectories and averages over the cyclones' lifetimes. Because the majority 487 of cyclones in perturbed simulations did not deviate from control for most of their courses, 488 most MTDs are measured as less than the length of the domain grid spacing of 30 km. 489 Directional MTDs considering deflection toward the North Pole or the perturbed snow 490 line were inconsistent with few outliers.

491
The study by Elguindi et al. [10] wherein snow was added to a Great Plains nested 492 domain two days prior to cyclone arrival generated similar trajectory outcomes with de-493 viations in perturbed cases only rarely exceeding 100 km. The trajectory deviations in 494 these tests, like our own, varied substantially. It is reasonable to infer those differences in 495 trajectories between control and perturbed cyclones included both a forced response from 496 the change in boundary condition, along with chaotic reactions to considerable energy 497 disturbances caused by changes to surface conditions over extensive areas in advance of 498 the cyclone, obscuring the functional responses to the specific positioning of snow cover. 499 However, our approach to average across cases and look at sensitivity in multiple cases 500 helps distinguish the forced response from natural variability while also maintaining a 501 snow line where we might expect enhanced low-level baroclinicity. The finding of limited 502 trajectory changes stands in contrast to significant cyclone responses to snow anomalies 503 found by multiple observational studies (e.g. [19,20,21]  Like MTDs, changes to cyclones' central low SLP due to a retreating snow line were 507 minimal. This, however, differed from the results of Elguindi et al. [10] who found an av-508 erage positive difference of 4 hPa in response to expanded snow cover, a value which only 509 two simulations in this whole study exceeded. Perhaps this can be attributed to the fact 510 that they added snow as opposed to removing it or to the physics of the MM5 model 511 compared to WRF-ARW. Additional tests of physics parameterizations including of land 512 surface model parameters, boundary-layer scheme, and radiation scheme, may help re-513 solve differences and provide additional insight into sensitivity of results to boundary 514 layer thermodynamics. Even with the disparity in the magnitude of pressure changes, 515 their discovered trend of central pressure increasing when snow is added is comple-516 mented by the findings of this study where snow removal generally contributed to a de-517 crease in central pressure. The deeper central low SLP while in transit over regions where 518 snow had been removed corroborates the conclusion of Elguindi et al. (2005), who noted 519 that snow cover prevents the deepening of mid-latitude depressions by reducing warm 520 sector temperature and moisture gradients and weakening surface convergence and 521 fronts. This finding is also consistent with climate model projections of increasing cyclone 522 intensity and precipitation extremes in future warmer climates, especially over land 523 [61,62].

524
Seasonal differences showed a stronger response in mid and later season, although 525 we admit there were a limited number of cases (5) per season studied (early, mid, late) to 526 make strong claims. Generally, responses of most investigated variables were greatest in 527 the mid-season, weaker in late season (with the exception of mean MTD), and weakest in 528 early season. This is counter-intuitive from an inspection of snow line retreat and incon-529 sistent with the argument that albedo feedback is the dominant driver. If anything, there 530 appears to be an inverse relationship between amount of mean retreat and response of 531 cyclones to the correspondingly-shifted snow lines. However, it has been shown that the 532 surface temperature effect of snow cover is strongest in late winter, likely a result of 533 stronger solar radiation enhancing albedo effects [22]. These responses suggest a greater 534 focus on seasonality may help refine our understanding of snow extent on mid-latitude 535 cyclones.

537
Fifteen cold season extratropical cyclones over or near the North American Great 538 Plains were examined in a series of simulations with varying percentile of snow retreat 539 consistent with 21 st century climate model projections and differing model initialization 540 lead times in order to gauge the dependence of their trajectories, intensities, and associ-541 ated precipitation on underlying snow cover and experimental design. When a realistic 542 retreat of snow cover consistent with climate warming scenarios [63] was applied to these 543 cases, with more than two days lead time, a majority of cyclones experienced a small de-544 crease in SLP (82%), consistent increases in precipitation (86%), modest increases in kinetic 545 energy, and limited changes in overall trajectory (mean of < 30 km). These results are 546 muted compared to expectations gained from observational studies such as that of Na-547 mias [11] and Rydzik and Desai [21], but reflect the results of modelling done by Elguindi 548 et al.
[10], manifesting a continued disagreement among models and observations. Com-549 pared to earlier studies, our study confirms that the model results are robust to assump-550 tions about varying levels of snow cover retreat, initialization time, and type of mid-lati-551 tude cyclone.`

552
It is yet unknown why the cyclone trajectories did not adhere more closely to shifted 553 snow lines, as the findings of observational studies suggest. Weaker responses to the re-554 moval of snow cover at the time of cyclogenesis suggest that the presence or absence of 555 the snow margin has a minor or counteracting effects or modeled land-atmosphere feed-556 backs are weak or constrained by other processes. Trajectory deviation, pressure change, 557 and precipitation intensify with greater lead time, stabilizing in most statistics after 2 days 558 lead time for model initialization and boundary forcing change, similar to Elgundi et al 559 [10]. Model fidelity for simulating of snow albedo has also been called into question [64], 560 and new approaches have been developed [65]. The full extent of the snow margin's in-561 fluence cannot be answered until longer case study simulations are executed.

562
Lingering questions remain on mechanisms of snow cover SLP, differences among 563 cases in surface energy-balance and radiative properties and its influence on cyclone dy-564 namics, and upper-level dynamics. Some of these, especially upper-level dynamics, are 565 studied in individual cases in detail in a companion paper [6]. Given prior reported sig-566 nificant snow-atmosphere feedback hotspots in central North America, continued analy-567 sis of its impact to storm tracks is warranted, including over a larger number of cases 568 [66,67]. The simulation model outputs here provide a rich data set for future evaluation 569 and are provided at the archive below for public access [68].

570
Supplementary Materials: The following supporting information can be downloaded at: 571 www.mdpi.com/xxx/s1, Figure S1: Observed snow extents, 10 th percentile; Figure S2: Observed 572 snow extents, 50 th percentile; Figure S3: Observed snow extents, 90 th percentile; Figure S4: Impact of 573 model lead time; Table S1: Average MTD; Table S2: Pressure decrease.    Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual au-729 thor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to 730 people or property resulting from any ideas, methods, instructions or products referred to in the content.