Persistent cold air outbreaks over North America in a warming climate

This study examines future changes of cold air outbreaks (CAOs) using a multi-model ensemble of global climate simulations from the Coupled Model Intercomparison Project Phase 5 and high resolution regional climate simulations. Overall, climate models agree on a dip in CAO duration across North America, but the percentage change is consistently smaller from western Canada to the upper mid-western US with historically more frequent CAO. By decomposing the changes of the probability density function of daily surface temperature into changes due to mean warming and changes in standard deviation (std) and skewness/higher order moments, the contributions of each factor to CAO changes are quantified. Results show that CAO changes can be explained largely by the mean warming, but the decrease in temperature std contributes to about 20% reduction of CAO from Alaska to northeastern US and eastern Canada possibly due to the Arctic amplification and weakening of storm track. A thermodynamical modulation of the skewness called the ‘0 °C mode’ effect is found to operate prominently along the 0 °C isotherm hemispherically and reduce CAO in western and northeastern US with winter snow cover by up to 10%. This effect also produces a manifold increase in CAO events over the Arctic sea ice. An increased frequency in atmospheric blocking also contributes to increases in CAO duration over Alaska and the Arctic region. Regional simulations revealed more contributions of existing snowpack to CAO in the near future over the Rocky Mountain, southwestern US, and Great Lakes areas through surface albedo effects. Overall, the multi-model projections emphasize that cold extremes do not completely disappear in a warming climate. Concomitant with the relatively smaller reduction in CAO events in northwestern US, the top five most extreme CAO events may still occur, and wind chill will continue to have societal impacts in that region.


Introduction
Extensive studies including those summarized in the IPCC Fourth Assessment Report (AR4) [1] have indicated a high probability for longer duration and more frequent heat wave as a result of increasing greenhouse gases projected for the future [2][3][4][5][6][7]. At the same time, cold extremes may not vanish, and they may cause more deaths than extreme heat. Keatinge [8] showed that in Britain, an average temperature increase of 2°C may lead to 2000 more heat related mortality, but reduce cold related mortality by about 20 000. Besides mortality, the extreme cold air outbreak (CAO) in 1983 and 1985 destroyed a majority of citrus trees in Florida that resulted in about $3 billion in economic losses in each event [9,10].
CAO has been attributed to atmospheric dynamics and large scale circulation anomalies associated with the North Atlantic Oscillation (NAO) and the Pacific-North American (PNA) teleconnection pattern during El Niño-Southern Oscillation (ENSO) events. In the negative phase of NAO, weaker zonal circulation could increase the chance of cold outbreaks over almost the entire northern hemisphere, including the US, Europe and Asia [11,12]. On the other hand, the onset of CAO over the eastern US is found to be accompanied by the development of the positive phase of the PNA [10], while the warm phase of ENSO may suppress extreme cold events throughout a major part of northern North America [13].
Recently, it has been hypothesized that global warming could perturb the polar vortex and enhance both warm and cold extremes [14,15]. The most recent severe CAO resulting from the polar vortex occurred in January 2014 [16] that affected large regions of the eastern United States. Using simulations from seven General Circulation Models (GCMs), Vavrus et al [17] found that although the frequency of CAOs generally decreases in the Northern Hemisphere, certain regions could experience more CAOs in the future. However, a more recent study [18] found that Arctic amplification could reduce winter daily temperature variance in the northern hemisphere, thus reducing the occurrence of CAO. Besides large-scale circulation, land surface conditions may also play a role in CAOs by providing heat sinks for the polar air masses that migrate with the cold fronts [12,19]. Vavrus [20] found that by removing the seasonal snow cover in a GCM, extreme CAOs were essentially gone as temperature increases both locally over areas with snow cover loss and in remote Arctic areas. As numerous studies have projected significant decline in snowpack under global warming, the potential influence of snow cover loss on CAOs should be considered in conjunction with large-scale circulation changes that may ensue to alter CAOs in the future.
A limited number of studies have investigated future changes in CAO using GCMs, but results to date are mostly based on a small number of GCMs. Recently, Westby et al [21] analyzed the historical simulations from 15 GCMs in the Coupled Model Intercomparison Project Phase 5 (CMIP5) [22] and found no significant historical trends in CAOs from 1949 to 2011 over the US. This study uses a larger multi-model ensemble to evaluate their simulations of historical CAOs, and examine their projected changes of CAO in the near future (2026-2045) and towards the end of the 21st century (2081-2100). We focus on North America where historically CAOs occurred quite often. Notably, CAO may change as a result of changes in many aspects including the mean, variance, skewness, and higher order moments of the surface air temperature frequency distribution. We analyze the contributions of these factors to CAO changes in the future and explore the dynamical and thermodynamical modulations of these factors to gain an understanding of processes linked to CAO changes. This study further analyzes regional climate simulations that have realistic depiction of topographic influence on snowpack to study the potential impacts of snow cover losses on CAOs in the future. Lastly, we investigate changes in both temperature and wind speed during cold extreme events from the CMIP5 ensemble to investigate changes in wind chill temperature (WCT) that have important implications to cold related mortality [23,24].

Model description and analysis methods
A total of 26 CMIP5 models with 39 ensemble members are used in this study, together with regional climate simulations from the Weather Research and Forecasting (WRF) model [25] that downscaled the Community Climate System Model (CCSM4) in the CMIP5 archive for one member of the present climate and two members of future climate scenarios. More details of the WRF simulations are described by Gao et al [26]. The CMIP5 models used in this study are summarized in table S1 in the supporting information, available at stacks.iop.org/ERL/10/044001/mmedia. All global model outputs are interpolated to 2.0 by 2.0°l atitude-longitude grids. The WRF simulations are performed at a spatial resolution of 20 km by 20 km. Both low-to-medium (RCP 4.5) and fossil fuel intensive (RCP 8.5) scenarios [27,28] were used in this study.
There is no universal definition of CAOs, so we adopt the definition used by Vavrus et al [17]. A CAO event is defined as two or more continuous days with daily mean near surface air temperature (T) two standard deviations (stds) (σ) below the winter mean temperature (T ) (i.e., T < T − 2σ). To calculate the mean and std, we calculate the std of each day over multiple years (i.e., 20 years), and then calculate the average of the std of each day in winter (i.e., 90 days) as the mean std. The winter CAO duration is calculated using the total number of CAO days divided by the number of years (20 in this case). The historical period is defined as 1981-2000, while two future periods are used: 2026-2045 and 2081-2100. The CMIP5 archive only includes three hourly data from these periods for most models used in the later part of this study. Similar to Vavrus et al [17], the same threshold determined from the present day simulations is used to define CAO in the future. In addition, we also use thresholds (T − 2σ) calculated from the future climate to detect CAO in the future climate simulations. As will be discussed later, this latter approach is useful for analyzing CAO changes related to changes in temperature skewness in areas where significant warming in the future eliminates most CAO defined using the threshold of the present climate. Future CAOs detected by this method are warmer than those detected using the present-day threshold, but arguably they may still have impacts on humans and ecosystems that could adapt to the warming climate.
Under global warming, changes in the mean, std, skewness and higher moments of the temperature frequency distribution can all impact the occurrence and duration of CAO. In order to quantify the respective contributions from these factors to the changes of CAO, we decompose the changes of CAO into changes attributed to these factors by constructing two pseudo warming scenarios based on statistical transformation of the daily temperature in the present-day climate for comparison with the model simulated future scenarios. The first pseudo warming scenario is generated by adding the mean warming of a future climate scenario to the simulated present-day daily temperature. The resultant probability density function (PDF) of the transformed daily temperature has the same shape as the PDF of the present-day climate, but with the mean temperature shifted to the right. Besides mean warming, decrease/increase of the std of temperature should reduce/increase the frequency of CAO. To isolate the contribution of changing std to CAO, we develop the second pseudo warming scenario by rescaling the present-day temperature anomalies (with respect to its winter seasonal mean) by the ratio of the std of the future scenario to the std of the present day. Adding the mean warming of the future scenario to the rescaled temperature anomaly, the transformed daily temperature has the same mean warming and std as the future scenario simulated by the models. The same approach was employed by de Vries et al [29] for evaluating the change of winter cold spells in Europe under global warming.
Comparing the changes of CAO from the two pseudo warming scenarios with the changes of CAO from the future scenario simulated by the models, we investigate the contributions from changes in mean, std, and higher order moments to the changes of CAO. This method has a limitation, however, as under significant mean warming, the first pseudo warming scenario may eliminate all occurrences of CAO defined by the thresholds of the present climate. Hence the second pseudo warming scenario cannot be used to further discern the contributions from the std and higher moment to CAO changes in the future. As will be shown later, this is a particular challenge in the high latitudes and polar region where warming is amplified. Nevertheless, evaluating the different moments of the temperature PDF is meaningful as climate change can influence CAO through profound changes of these properties, as detailed in sections 4.1 and 4.2.
Because of the caveat discussed above, the effect of skewness and higher moments on CAO changes is further elucidated using an alternative approach that defines CAO in the future scenarios using thresholds (T − 2σ) derived from the future climate. We refer the CAO detected using these future thresholds as 'warm world' CAO events, although these events are still very cold in the high latitudes and polar region. Comparing the warm world CAO events detected from the second pseudo warming scenario, which has the same mean temperature and std as the simulated future climate scenario, with the warm world CAO events from the simulated future climate scenario, we can attribute CAO changes to changing skewness and higher moments in the future climates. This method is used in section 4.2 to study the changes of skewness and their contributions to CAO changes in the future.
The mean winter CAO durations are shown in figure 1 and the number of CAO events per winter is shown in figure S1 in the supporting information. Overall, the spatial distributions from the CMIP5 ensemble mean and CCSM4 and WRF are consistent with the observations, showing two large areas with CAO occurrences. A swath of high CAO frequency stretches from southwestern Canada, northwestern US to the mid-western US, with mean CAO duration of four to five days per winter. As will be shown later, these areas with high CAO frequency have a strong correspondence with large negative skewness in the temperature PDF. Since there is good agreement among CMIP5, CCSM4 and WRF, the following sections mainly describe results from CMIP5 to indicate model agreement and robustness of future changes.

Changes of CAO duration and the dynamical and thermodynamical modulation 4.1. Changes of CAO duration
In this section, changes in CAO duration are calculated for each model by comparing CAO durations for the future and present using a fixed temperature threshold (called present threshold) determined from the historical simulation (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and applied to both present and future climate. To evaluate the robustness of CAO changes in the future, both model agreement and statistical significance are used to assess the changes simulated by the 26 CMIP5 models for a total of 39 members. The multi-model results on each grid are considered to have model agreement if at least 75% of the CMIP5 models show the same sign of change as the CMIP5 ensemble mean. For models with multiple ensemble members, the ensemble means are used. Statistical significance is then tested for each model that agrees with the ensemble mean; if at least 75% of these models show statistical significance at 95% level, then the ensemble mean change for that grid is considered statistically significant. If a model has multiple members, statistical significance is achieved only if at least half of the members show significance.
Figures 2(a) and (b) show the percentage change of winter CAO durations for RCP 4.5 2030s (low warming) and RCP 8.5 2090s (high warming) using the present threshold, while the corresponding results for RCP 8.5 2030s and RCP 4.5 2090s with intermediate warming are shown in figure S2. One of the most robust features is the model agreement on the decrease of winter CAO durations across the entire North America in the future (both RCP 4.5 and RCP 8.5). However, areas from western Canada to the upper Midwestern US that have high CAO duration historically will still experience relatively more CAO days or less fractional reduction in CAO days than other regions. For example, while the percentage decrease of CAO over areas with historically high CAO is 50% or less in the 2030s, the CAO over areas such as northern Alaska, northern and eastern Canada that experience less frequent CAO historically will see a decrease of about 60-80% in the 2030s under RCP 4.5, and CAO will likely disappear by the end of the 21st century under RCP 8.5.
From the decomposition analysis described in the section 2, figures 2(c) and (d) show the CAO changes comparing the first pseudo warming scenario with the present climate. Resemblance between figures 2(c) and (d) and figures 2(a) and (b) suggest that the model simulated CAO duration changes can be largely explained by the shift of the PDF or mean warming captured by the first pseudo warming scenario. Areas with more/less CAO events in the present climate have temperature that corresponds to relatively higher/ lower percentile, so a shift of the PDF by the mean warming without shape change would bring smaller/ greater reduction of CAO fractionally. Figures 2(e) and (f) show the difference in CAO changes between the second and first pseudo warming scenarios to isolate the effects of std changes on CAO in the future. Changes of CAO attributed to std changes are characterized by a reduction over a broad region extending from Alaska to northeastern US and eastern Canada with peak magnitude as large as −20%. The CMIP5 models demonstrate considerable agreement on the sign of the CAO changes over these areas. Adding the effects of PDF shift and std changes can accurately reproduce the total percentage change of CAO between the future and present climate. Similar conclusion was also reached by de Vries et al (2012) for the change of cold spells over western Europe in the CMIP5 future warming simulations.
It is reaffirming to note that the CAO change pattern in figure 2(f) conforms well with the corresponding pattern of temperature std reduction (see figure S3 (b)). In particular, wherever there is a reduction of CAO due to std changes, the std of temperature itself is also reduced in the same areas. The influence of climate warming on the std of surface temperature has been well established and linked to the reduction of the local temperature gradient (e.g., [18,35]). Figure  S3(a) shows the CMIP5 ensemble mean surface temperature warming trend by the end of the 21st century under RCP 8.5-a pattern often referred to as polar amplification emerges. The correspondence between the reduced meridional temperature gradient inferred from figure S3(a) and the decrease of std ( figure S3(b)) is easy to discern. No doubt that a weakened poleward temperature gradient can reduce cold advection and consequently the variance anomalies associated with this mechanism [35]. Coincidentally, studies of the mid-latitude storm track response to global warming (e.g., Chang et al 2012 [36], see their figure 6(a); [37]) showed reduction of storm activity (measured by the variance of the band pass-filtered surface pressure) over exactly the same areas as the CAO reduction induced by the muted std. Thus, it is conceivable that weakening of wind variance (inferred from the geostrophic relation from the change of pressure) might also play a part in the reduction of CAO in the North American continent.

Linkage between CAO and the skewness of temperature
The effect of changes of the higher moments of temperature PDF on CAO can be derived either as the residual after subtracting the effects of changes in mean temperature (figures 2(c) and (d)) and std (figures 2(e) and (f)) of the PDF from the CAO changes in the future (figures 2(a) and (b)), or using the second approach explained in section 2 using future thresholds to detect the 'warm world' CAO. Consistent results from the former and latter approaches are shown in figures 2(h) and 3(d), respectively, for RCP 8.5 in the 2090s. Over the contiguous US, both methods reveal negative contribution to the occurrence of CAO from the higher moments with the impact centered near the western and northeastern US. Their magnitudes are relatively small compared to the effect of mean warming and reduction of std, but still sizable with the residual approach showing contributions up to 10% reduction of CAO and the future threshold approach showing about 20% reduction of 'warm world' CAO. The different magnitudes of contribution should be attributed to the different criteria used to define CAO, but the same spatial pattern that emerged from the two methods argues for some physical basis. Almost the same pattern of CAO changes but with somewhat weaker magnitude is also found for RCP 4.5 in the 2090s (not shown).
Over the Arctic where only the 'warm world' CAO is pertinent, the higher moments effect on CAO is a manifold increase in CAO duration. The cause for this is revealed by the PDFs of temperature examined next. Remarkably, the pattern of the 'warm world' CAO changes ( figure 3(d)) is almost a perfect mirror image of the change of the temperature skewness for the same scenario ( figure 3(c)), suggesting a strong link between the temperature skewness on CAO. It remains to explain why the skewness varies as such from the models. Overall, the skewness change in figure 3(c) compared to the present-day skewness in figure 3(b) can be categorized as (i) a reversal from weak positive to negative, which occurs typically over the Arctic sea ice, and (ii) an increase towards positive skewness in areas with climatological mean negative skewness, which typically occurs over the mid-latitude band. An immediate reason for these skewness changes is the melting of snow and sea ice. This is hinted by the first category of changes occurring over the Arctic with sea ice, and the second occurring near the climatological 0°C isotherm shown in red in figures 3(a) and (c).
We pick two spots representative of the two categories discussed above to investigate the skewness changes. For the first category, we choose an area north of Alaska between the Beaufort Sea and East Siberian Sea in the Arctic, where the present-day temperature PDF is slightly positively skewed ( figure 4(a)). As climate warms, it transitions to a strongly negatively skewed PDF with the mode of the PDF progressively shifted towards 0°C at the end of the 21st century under RCP 8.5. For the RCP 4.5 scenario, the PDF ends up with a bimodal distribution with a new peak formed near 0°C. Thus the increasingly larger negative skewness is partly a result of the dramatic increase in the frequency of daily temperature near 0°C and the limited shift above the melting point in the future. As sea ice is gradually warmed to near 0°C in the course of this century, additional energy that is needed to melt the sea ice dampens the warming rate, causing an accumulation of occurrences of daily temperature or stagnation of the mode near the melting point. This '0°C mode' effect is apparent  Figure 4(b) shows the shifted PDF (red dashed line) from the first pseudo warming scenario, the rescaled PDF (purple dashed line) from the second pseudo warming scenario with mean warming and changes of std, and the simulated future scenario (red solid line) in the 2090s for RCP 8.5. It shows clearly the stark contrast in PDF shape caused by the change in skewness compared to the shifted and rescaled PDFs, contributing to the large increase in 'warm world' CAO events depicted in figure 3(d) over the Arctic.
For the second category of skewness change, we pick a box (blue box in figure 3(c)) over the eastern US where relatively large skewness-induced reduction of CAO is found and the surface is covered by snow during most of the winter. The corresponding PDFs are shown in figures 4(c) and (d). The '0°C mode' effect is also apparent in this region particularly in the near future (2030s) and in RCP 4.5 even in the 2090s. That is, despite the changes in mean temperature towards warming, the mode stagnates near 0°C in the future. Hence its effect is to increase the negative skewness of the PDF in the near future, similar to figure 4(a). However, the PDF for RCP 8.5 2090s and figure 3(c) show a change in skewness towards positive in a much warmer scenario where winter snow cover is abated. This is more apparent in figure 4(d) that compares the shifted and rescaled present-day PDF against the future PDF. As warming continues in areas that are snow free, while snow-covered areas that remain anchor the mode near 0°C, the future PDF has a broadened peak, and positive skewness develops. Hence, similar to sea ice, melting of snow has a thermodynamical modulating effect on the shape of the PDF that results in changes in temperature skewness, but the changes can be positive or negative depending on the rate of warming and the extent of sea ice and snow cover simulated by the models. Interestingly, the 0°C isotherm of the climatological mean temperature nicely demarcates the area where the '0°C mode' effect operates, and it meanders right through the center of the reduction of negative skewness over the US and central Europe in figure 3(c). Similar good correspondence with the 0°C isotherm also holds for CAO changes in figure 3(d).
It is well known that extreme cold events during winter are related to the strong advection by circulation regimes [38], such as blocking. It is conceivable that some of the changes in the CAO, especially those accompanying the skewness change, are related to the changing statistics of blockings. To test this, we plot the present-day frequency of blockings and its change under RCP 8.5 by the end of the 21st century (defined and detected in Masato et al 2013 [39]) as contours overlaying the corresponding patterns of present-day skewness and its change in figures 3(b) and (c), respectively. It is encouraging to note some interesting agreements between the two fields in the present-day climatology, such as the northeast Asia-north Pacific blockings and the negative skewness extending from the northeast Pacific to the west coast of the US and Canada (weighted towards the land mass though), and correspondence between the western European blockings and the skewness over western and central Europe. Although the exact geographical relationship between the location of blocking and cold temperature extremes is not yet well understood [40], the colocation of the dipolar change in the frequency of blocking in the northeast Asia-northern north Pacific sector (an increase over Alaska and the surrounding seas and reduction to the south over the Pacific) with the pattern of the skewness change in the RCP 8.5 2090s scenario (figure 3(c)) is unlikely fortuitous, but suggests a dynamical modulation of skewness and CAO changes in the future. Admittedly, the exact temperature fingerprints of the possible impact of wintertime blockings on the skewness of wintertime temperature over different regions should be investigated to consolidate and understand the apparent links between blockings and skewness discussed above, a task beyond the scope of the current investigation.

Effect of snow water equivalent (SWE) on CAO duration from WRF-CLM results
The above analysis indicates that large-scale circulation, i.e., atmospheric blocking is closely related to the change of CAO. Furthermore, the '0°C mode' effect appears to provide strong thermodynamical modulations of temperature skewness and CAO changes in areas with sea ice and snow cover. Snow cover has a strong cooling effect on surface air temperature due to its high albedo [41] and can potentially impact atmospheric circulation [42]. CAOs are found to correlate highly with snow cover [20]. Vavrus et al [20] conducted a sensitivity experiment with all land snow cover eliminated in a GCM. This resulted in 8-10°C increase in winter surface air temperature over northern North America, which eliminated all occurrences of CAO defined based on the threshold from the control run. To investigate the impacts of snow cover on CAO, we take advantage of the high resolution WRF [26] simulations that better resolve mountain snowpack than the CMIP5 models. In addition, WRF outputs include 6-hourly snow cover while CMIP5 outputs only include monthly snow cover, which does not have sufficient temporal resolution to determine changes in snow cover during CAO events.
Observed SWE is downloaded from the National Aeronautics and Space Administration (NASA) Earth Observations (NEO: http://neo.sci.gsfc.nasa.gov) and used to compare with the WRF model outputs ( figure  S4 in the supporting information). WRF generally well captured the spatial distributions of SWE over mountainous areas such as the Sierra-Nevada mountain range, Cascade mountain range, Rocky Mountains, and in eastern Canada. However, the model underestimates SWE over a band stretching from mid-western US to western Canada.
Winter mean SWE shows a predominant decreasing trend in North America under warming, except for small increases near Hudson Bay in northern Canada (not shown). However, figure 5 shows that the existing SWE one day before the onset of CAO events in the future is higher than that in the present climate, particularly for RCP 4.5 in the 2030s ( figure 5(e)). Using SWE two days before the onset of CAO events for the analysis results in patterns similar to figures 5(e)-(h). Apparently, from figure 5(e), higher SWE in the Rocky Mountains, southwestern US including western Texas, New Mexico, Colorado, intermountain West such as Utah, and near the Great Lakes may contribute to the occurrence of CAO in the near future (i.e., 2030s in RCP 4.5; figure 5(a)) by increasing surface albedo and reducing surface temperature during CAO, despite an overall reduction in winter mean SWE with increasing temperature.
Note that from figure 5(a), the CAO in RCP 4.5 in the 2030s simulated by WRF-CLM increases over the Rocky Mountains and southwestern US, although the multi-model CMIP5 shows a dominate decreasing trend of CAOs (figure 2(a)). Analysis of CAO changes from CCSM4 member 6 that provided initial and boundary conditions for the WRF simulations produces similar patterns as the regional simulations (not shown here). Noting that the CAO changes simulated by WRF and CCSM4 show important differences compared to the CMIP5 ensemble mean, analysis of SWE from the multi-model CMIP5 ensemble may provide further insights on the relationship between SWE and CAO. However, daily SWE is not currently available from the CMIP5 archive. As SWE reductions are more significant in RCP 8.5 or in the 2090s, the likelihood of snow cover contributing to CAO events becomes negligible (figures 5(b) and (f): 2090s in RCP 4.5 and figures 5(d) and (h): 2090s in RCP 8.5) except near the Great Lakes including Ohio, Indiana and Illinois where increases of SWE are still apparent before the onset of CAO events even in RCP 8.5 and near the end of the century. These results demonstrate that snow cover may play an even more important role in CAO events in the near future when snow in the mountains is still dominant in the cold season and its impacts on surface temperature may become critical for cold air masses to reach a temperature below the thresholds of CAO. However, these analyses do not delineate the relative contributions of snow cover and circulation anomaly to CAO. Arguably future CAO events are likely associated with both larger large-scale circulation anomaly and existing snow cover to overcome the effects of global warming.

Future occurrence of the top five historical coldest events
Since climate models projected continued occurrence of CAO events even in a warmer climate, an important question is whether the top five historical CAO events may still occur in the future. To explore this possibility, CAO events are ranked by the mean temperature of each event and the fifth lowest temperature of the present climate is used as a threshold to determine future occurrence of the historical top five events. Each occurrence of future CAO event with mean temperature lower than the threshold is counted towards the likelihood of the top five historical coldest events happening in the future. Figure 6 shows that about 2-3 top five historical events are projected to occur in the near future (2030s) in both RCP 4.5 and RCP 8.5. These future occurrences of the top five historical events are found in western US where historically strong CAO events occur (figure 1) and CAO durations are projected to decrease by less in future (i.e., 2030 in RCP 4.5 in figure 2(a)). In the 2090s, although with less model agreement, the top five historical events may still occur once or twice. Similar analysis has been performed for the top longest lasting events (figure S5 in the supporting information). At present, the longest CAO event lasts around ten days in particular over NW US, and western Canada. In the future, although the longest event generally becomes shorter, it could still last five days or more in the 2030s in both RCP 4.5 and RCP 8.5, and even in the 2090s in RCP 4.5.

WCT and wind chill warning
Cold extreme events have large impacts on human health. However, such impacts are more closely related to WCT which combines the effects of cold temperatures and strong winds. We adopt the formula of Osczevski and Bluestein [43] shown in equation (1)  We compare historical CAO events using daily mean, minimum and maximum WCT derived from three hourly WCT, and compared with the CAO determined using daily average temperature, shown in figure S6 in the supporting information. The CAOs determined by these four metrics have similar spatial distribution and magnitude. Comparing the future changes in CAOs, the four metrics again provide similar spatial distribution of changes (not shown). Since WCT is closely related to wind speed, the changes of wind speed in the future are also evaluated. The winter mean wind speed tends to decrease in North America ( figure S7). However, focusing on the days of the CAO events, there are larger areas showing increase of wind speed (figure S8). In particular, wind speed near the Great Lakes could increase on average by as much as 1 m s −1 (3.6 km h −1 ), which corresponds to almost 1°C decrease in WCT.
Wind chill warning occurs when WCT is low enough to be life threatening. There is no universal threshold to determine wind chill warning. In this study, we use the thresholds from Environment Canada 3 : when WCT is below −28°C, exposed skin can freeze in 10-30 min. Thus, we use −28°C as the threshold for areas south of 45°N near the US-Canada border, and reduce the threshold by 1°C for each degree of latitude to the south. Figure 7(a) shows less than two days of wind chill warning days per winter in the southern US. In western Canada and NW US with historically higher CAO frequency (figure 1), the wind chill warning days reach 4-10 per winter. In the future, the number of wind chill warning day decreases predominantly across North America with a few scattered exceptions. Similar to CAO duration, a swath across the northwestern to mid-western US shows smaller reduction by 10-20% compared to the surrounding, (figures 7(b)-(e)). Hence wind chill warning may still have important societal impacts particularly over northwestern US to central US.

Conclusions
Using a large ensemble of global climate simulations, we found robust decrease of CAO duration in the future due to climate warming with the magnitudes of decrease differ by locations. In particular, over northwestern US, the decrease is relatively smaller than other areas. In addition to the mean warming, the changes of CAO are highly related to the decrease of std and changes of skewness. By decomposing the changes of the PDF of daily surface temperature into changes due to mean warming, changes in std, and changes in skewness and higher order moments, the contributions of each factor to CAO changes in the future are quantified. Combining the decomposition analysis with analysis of future CAO detected using thresholds defined by the future climate that we called 'warm world' CAO, and analysis of PDF changes in two locations with sea ice and snow cover, we further examined the changes in skewness and the modulations by dynamical and thermodynamical processes.
The decomposition analysis shows that CAO duration changes in the future can be largely explained by the mean warming, but the decrease in temperature std contributes to as much as a 20% reduction of CAO over a broad region extending from Alaska to upper Midwestern and northeastern US and eastern Canada. Our results are consistent with recent studies that attributed the decrease of std of temperature to the Arctic amplification [18] and weakening of storm track extending from Alaska, Canada and across the continental US to western North Atlantic [36], although the exact dynamical modulation has not been completely clear so far [44]. Changes in skewness also have non-trivial impacts on CAO and reduce CAO in the western and northeastern US with climatological snow cover in the cold season by up to 10% and 20% using present and future thresholds, respectively. We identified a thermodynamical modulation of the skewness from additional energy needed to melt the snowpack that we called the '0°C mode' effect that operates prominently along the 0°C isotherm hemispherically. The same effect over the Arctic sea ice produces a manifold increase in CAO events detected using the future threshold. The change of skewness is also modulated dynamically through changes in atmospheric blocking, which affects CAO in regions such as Alaska. In particular, the dipole change pattern (increase north of Alaska and decrease south of Alaska) of atmospheric blocking frequency has a counterpart of dipole change in skewness over the same region. Combining the '0°C mode' and atmospheric blocking effects, 'warm world' CAO duration could increase by 3-10 times over Alaska and the Arctic by the 2090s in RCP 8.5.
Taking advantage of our high resolution regional climate simulation, we further investigated the effect of snow cover on CAO changes in the US. In the future, while the winter mean SWE shows predominate decreasing trends, we found the likelihood for CAO to be enhanced by increased antecedent SWE during the onset of CAOs. This demonstrates that future CAO events are more likely to occur on days with existing SWE in the Rocky Mountains, southwestern US and near the Great Lakes than CAO events in the present climate, as snow increases surface albedo and the resulting reduction in surface temperature could be critical for the warmer future cold air masses to cross the thresholds of CAO.
Using three hourly data, we also evaluated potential changes in wind chill warning. We found over western Canada and northwestern US with smaller percentage decrease of CAO days, the decrease of wind chill warning days is also concomitantly smaller than the surrounding areas. In addition, the top five coldest events simulated in the present climate are projected to occur a few times particularly over western Canada and northwestern United States. Overall, our results from multi-models emphasize that cold extremes do not completely disappear in a warmer climate. Depending on the future emission scenarios and time period, CAOs and the associated societal effects may continue to be an important issue for policy makers.