This study used air temperatures from a suite of regional climate models participating in the North American Climate Change Assessment Program (NARCCAP) together with two atmospheric reanalysis datasets to investigate changes in freezing days (defined as days with daily average temperature below freezing) likely to occur between 30-yr baseline (1971–2000) and midcentury (2041–70) periods across most of North America. Changes in NARCCAP ensemble mean winter temperature show a strong gradient with latitude, with warming of over 4°C near Hudson Bay. The decline in freezing days ranges from less than 10 days across north-central Canada to nearly 90 days in the warmest areas of the continent that currently undergo seasonally freezing conditions. The area experiencing freezing days contracts by 0.9–1.0 × 106 km2 (5.7%–6.4% of the total area). Areas with mean annual temperature between 2° and 6°C and a relatively low rate of change in climatological daily temperatures (<0.2°C day−) near the time of spring thaw will encounter the greatest decreases in freezing days. Advances in the timing of spring thaw will exceed the delay in fall freeze across much of the United States, with the reverse pattern likely over most of Canada.
The seasonal occurrence of freezing conditions is an integral element of a region’s ecosystem processes, recreational activities, and economy. For example, freezing air temperatures T play a key role in the overwinter survival of many insects that must rely on external sources to provide their heat (Lee 1989; Bale and Hayward 2010). The primary driver of the interannual variability in seasonal transitions and length is regional seasonal T variability (McCabe et al. 2015). Across much of North America a decrease in the number of days when minimum daily T falls below 0°C, known as “frost” days, has been noted for the period 1951–2003 (Alexander et al. 2006). McCabe et al. (2015) examined daily minimum temperature data for the conterminous United States (CONUS) for the period 1920–2012 from the Global Historical Climatology Network and noted earlier spring final frost dates after about 1983, with a change to later fall freeze most noticeable after about 1993.
Other changes in seasonal metrics over North America have occurred. Kunkel et al. (2004) analyzed daily temperature observations from 1980 to 2000 across the contiguous United States and found increases in frost-free season length of approximately 1 week. Contraction of the frozen season has been concurrent with a shift toward earlier spring thaw across much of North America (Schwartz et al. 2006; Wang et al. 2011). Declines in frozen season length and earlier spring thaw are known responses to warming, which across the United States has been strongest and most extensive in spring (Mutiibwa et al. 2015). Advances in the timing of spring thaw have exceeded the delay in autumn freeze-up across most of North America, possibly owing to feedbacks involving losses in snow (Kim et al. 2015). Snow-cover decreases across North America from 1972 to 2011 were greatest in spring and summer, with no detectable or consistent trend in fall snow cover. The duration of the snow season (first autumn snowfall to last spring snowfall) declined by 5.3 days decade−1 from 1972/73 to 2007/08, driven primarily by an earlier season end in spring (Choi et al. 2010). Snowfall has declined sharply across the western United States (Kunkel et al. 2009), likely associated with recent spring season warming. Advances in the timing of spring thaw, however, are not ubiquitous. Ault et al. (2015) reported that across the United States in recent decades, spring onset is not advancing uniformly and that later spring thaw dates have occurred over the northwestern Cascades from 1979 to 2013. They further suggested that interannual to decadal variations appear to pace regional trends.
The trends are expected to continue. Meehl et al. (2004) examined possible future regional changes in frost days in a global coupled model and found that patterns of relative changes of frost days are indicative of regional-scale atmospheric circulation changes that affect nighttime minimum temperatures. An analysis of statistically downscaled model outputs (with bias correction) between historical (1950–2005) and late twenty-first-century (2071–2100) periods suggested that leaf out (an indicator of spring onset) is projected to shift earlier by 22.3 days across the conterminous United States by the end of the century (Allstadt et al. 2015). Simulations with coupled atmosphere–ocean general circulation models (AOGCMs) indicate that temperature increases over high latitudes will exceed those at low latitudes (Holland and Bitz 2003), with this arctic amplification influenced by declines in sea ice (Serreze et al. 2009) and snow cover (Vavrus 2007).
In this study, we use 2-m T drawn from outputs of the North American Regional Climate Change Assessment Program (NARCCAP; Mearns et al. 2007, 2009) and from atmospheric reanalysis to examine likely changes in freezing days ΔFD, defined here as days with daily average temperature at or below 0°C, between baseline (1971–2000) and midcentury (2041–70) periods across North America. We also examine influences on ΔFD and its spatial patterns along with associated changes in the timing of spring thaw and autumn freeze across the region.
2. Data and methods
a. RCM and reanalysis datasets
Daily gridded T estimates for 30-yr time slices in recent past and mid-twenty-first-century periods were drawn from regional climate model (RCM) outputs. In the NARCCAP program boundary conditions for six RCMs were provided by four AOGCMs for 30 yr of baseline climate (1971–2000) and 30 yr of a future climate (2041–70) forced by the IPCC Special Report on Emissions Scenarios (SRES) A2 (differentiated world) emissions scenario. Under this scenario, CO2 levels are projected to be 490 ppm in 2040, increasing to 635 ppm by 2070 (CMIP5 2016). A suite of 11 GCM–RCM pairings was available for analysis. NARCCAP outputs have been used to examine projections for extreme precipitation events (Gutowski et al. 2010) and likely future changes in T and precipitation across the northeastern United States (Rawlins et al. 2012; Fan et al. 2015). The Hadley Centre GCM uses a 360-day calendar. For the HadCM3–MM5 and HadCM3–HRM3 pairs we linearly interpolated values for missing days to a 365-day data series. The six RCMs participating in NARCCAP have different native grids. Analysis of RCM outputs frequently requires spatial interpolation of the data from the model’s native grid to a common one. We interpolated daily T for each GCM–RCM pair for each year to a 0.5° grid using an inverse-distance-weighted spatial interpolation method described by Willmott et al. (1985). The method is based on the original two-dimensional algorithm of Shepard (1968). For each GCM–RCM pair we then generated daily climatological T for both the 30-yr NARCCAP “baseline” period and the 30-yr midcentury “future” period on the 0.5° grid. An 11-member ensemble mean daily climatological T was then produced by averaging T values over the GCM–RCM pairs. Hereafter, and refer to these climatological daily T, averaged over the 11 ensemble NARCCAP members. Those terms are also used here for climatological daily T derived from two reanalysis datasets, as described below.
Biases in baseline T estimates will negatively affect the derived change estimates that are based on exceedance of defined thresholds such as freezing days FD. The NARCCAP models exhibit a slight cool bias across North America, though two of the six RCMs have a warm bias, particularly over Canada (Mearns et al. 2012). While observed meteorological station density is relatively high across the United States, coverage is low enough over areas of northern Canada to limit confidence in spatial and seasonal variations in T, and thus FD if only station data are used. A combination of observations and model assimilation known as atmospheric reanalysis can be useful in estimating daily long-term climatological T over broad regions such as North America.
To assess robustness in the spatial pattern and magnitude in ΔFD we used two reanalysis T datasets, from which respective and, in turn, were derived. One dataset is the North American Regional Reanalysis (NARR; Mesinger et al. 2006). The NARR assimilation scheme ingests near-surface observations hourly, with atmospheric profiles of temperature, winds, and moisture assimilated every 3 h. Daily data are available on a 32-km grid. When compared to observations from surface stations, T biases in the NARR have been found to be generally smaller and have less of a diurnal cycle (Mesinger et al. 2006) than biases seen in the earlier NCEP–NCAR Global Reanalysis 1 (Kalnay et al. 1996). NARR data include leap years. We removed 29 February values each leap year and subsequently interpolated the daily data to the 0.5° grid. NARR data are available beginning in 1979, so we generated gridded daily over the 22-yr period 1979–2000. We also leveraged the Water and Global Change (WATCH) forcing dataset (WFD; Weedon et al. 2011), which is based on the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40; Uppala et al. 2005). For the WFD, daily Ts on a global 0.5° grid were produced (by the creators) by extracting and interpolating ERA-Interim to 0.5° spatial resolution with a bias correction from gridded observations. From this product we derived daily for the same 1971–2000 period as NARCCAP, after removing leap days. The calendar-day adjustments are inconsequential for the derivation of 30-yr T climatologies. A fourth field of was produced by imposing a spatially uniform warming of 2.5°C onto calculated from WFD. A warming of 2.5°C is approximately a lower bound on projected midcentury temperature change averaged over North America (Karl et al. 2009).
b. Derived statistical metrics
Our analysis centers on ΔFD obtained between the baseline and future periods. For each NARCCAP GCM–RCM pair, FD was determined from daily climatological T for both the baseline and the future periods. Gridded estimates of for each of the NARRCAP models and for the NARCCAP 11-model ensemble mean T were derived simply from the model T outputs interpolated to the common 0.5° grid. For example, gridded ΔFDs from the NARCCAP ensemble were calculated as the difference in FD between and periods, with a negative ΔFD indicating a decrease, as expected, with time. Estimates of for NARR and WFD were obtained by applying the change in daily T as calculated from the NARCCAP ensemble T over the baseline and future periods () to the respective reanalysis field. This method of bias correction is referred to as the “delta method” (Hawkins et al. 2013). Application of this method can minimize unwanted effects of biases in present-day temperatures simulated by climate models.
We also determined likely future changes in the timing of spring thaw ST and autumn freeze AF from each set of and . The day of ST was determined as the first day when 12 out of 15 consecutive days (i.e., 80% rule), from January through June, had a climatological T > 0°C. The day of AF was selected as the first day when 12 out of 15 consecutive days from September through December were at or below 0°C. The 15-day moving window size and 80% rule was determined to provide consistent detection of these events across the study region, relative to in situ station observations (Zhang et al. 2011).
Characteristics of the seasonal T cycle that influence ΔFD, ΔST, and ΔAF were also examined. Baseline annual temperatures (gridded , for each dataset) were calculated as the simple arithmetic mean of the 365-day climatological daily T. To further understand influences on the freezing day and seasonal-transition metrics we also estimated the rate of change in climatological daily T (, °C day−1) over the two weeks before and two weeks after spring thaw. Last, climatological T normals (1971–2000) from several stations in the northeastern United States (National Oceanic and Atmospheric Administration 2014) were used to help explain spatial variations in ΔFD as air temperatures vary across a roughly north–south transect of stations.
3. Projected decreases in frequency of freezing days
a. Spatial patterns
Seasonal changes in 2-m T captured by the NARCCAP ensemble average and are shown in Fig. 1. Over North America winter warming ranges from approximately 2°C across the southern United States to over 4°C near Hudson Bay (Fig. 1a). Strong autumn warming is evident across far-northern Canada (Fig. 1d). The south–north winter warming gradient is consistent with expectations and aligns with the greater warming observed across high-northern-latitude land areas over recent decades (Serreze et al. 2000; Hinzman et al. 2013), driven in part by losses in snow cover (Choi et al. 2010). Interestingly, over the entire region, autumn is projected to warm more than spring (Figs. 1b,d).
Maps of ΔFD calculated from and for each individual NARCCAP model (Fig. 2) exhibit spatial patterns characterized by the greatest decreases (most negative ΔFD) along the southern boundary of the part of the continent that experiences seasonally freezing conditions and across far-western Canada. The declines are smallest over north-central Canada. To some extent they also reflect an influence of the driving GCM. The ΔFD s derived from daily and calculated from ensemble mean T show this spatial pattern more distinctly (Fig. 3a). Locally, declines in excess of 60 days are found across extreme western Canada and through the northwestern and central United States.
The spatial pattern in ΔFD obtained using NARCCAP applied to NARR and WFD (Figs. 3b,c) should reflect more accurate baseline fields of climatological daily T. For both NARR and WFD the maps point to greater decreases in freezing days over the western United States as compared to the results using NARCCAP data alone. For WFD the area of the continent that will experience freezing conditions declines by approximately 1 × 106 km2 (6.3%), from 15.7 to 14.7 × 106 km2 (domain area 20.6 × 106 km2), between the baseline and future periods (Fig. 4). The decline obtained from NARR is similar, approximately 5.7% (from 15.8 to 14.9 × 106 km2). Latitudinal transects through the region (Fig. 5) show ΔFDs declining most, in the eastern and central United States, near 40°N. In the west, there is a broad area of ΔFD exceeding −50 days from 40° to 50°N. Above 62°N, ΔFDs are above −20 days. In other words, more modest declines will happen that far north. While continental averages are similar (−22, −23, and −22 days for NARCCAP, NARR, and WFD), ΔFD values derived from NARCCAP T alone show a higher percentage of grids with more modest decreases (Figs. 6a–c). It is clear that across the central United States the zone of greater ΔFD is 1°–2° in latitude farther north when the reanalysis data were used to estimate and, in turn, , consistent with their slightly warmer temperatures compared to the NARCCAP models. Applying a spatially uniform warming of 2.5°C to the gridded from WFD results in a ΔFD pattern (Fig. 3d) similar to the WFD pattern shown in Fig. 3c.
The pattern in Fig. 3a follows from the averaging of climatological T over the available NARCCAP model GCM–RCM pairs. It represents the transient mean response to GHG concentration increases present in the A2 emissions scenario used in the simulations. The pattern in Figs. 3b,c estimated using the two reanalysis datasets reflects the changes arising through use of an ensemble of models and the improved initial fields of from NARR and WFD. Areas of the Intermountain West are biased slightly cold in the NARCCAP simulations (figure not shown). The changes derived from NARR and WFD thus point to a higher risk for large FD decreases in that region compared to the results using NARCCAP model T outputs alone. Areas with baseline mean annual T between approximately 2° and 6°C (Figs. 6a–c) are likely to experience the greatest freezing day declines.
b. Influences on ΔFD spatial variations
While seasonal warming clearly plays a role, other aspects of regional climate further explain the spatial patterning in ΔFDs. Spatial variations shown in Figs. 3a–d are related to the rate of change in climatological daily T near the time of spring thaw calculated from from each of the three datasets. Results show only modest future decreases in FD in the coldest areas with relatively higher (Figs. 6d–f). For both the NARR and WFD, ΔFDs exceeding −60 days are noted across parts of California, Oregon, Washington, and British Columbia where baseline is between 2° and 6°C and < 0.3°C day−1.
The spatial pattern and magnitude in ΔST and ΔAF provide additional insights into the anticipated future changes. Shifts to earlier spring thaw ΔST are found across much of the central and western United States (Figs. 7a,b), a result consistent with but larger in magnitude than recent observed changes for the Northern Hemisphere (McCabe et al. 2015). Differences between ΔST and ΔAF (Fig. 8) suggest that ΔST will exceed ΔAF over parts of the central and western United States. However, over much of northern Canada that tendency is reversed, particularly west of Hudson Bay. Across the conterminous United States the T trend (1895–2015) is 40% greater in spring compared to autumn. The ΔST and ΔAF rates are of similar magnitude to recent declines in snow season length, which averaged 5.3 days decade−1 between the winters of 1972/73 and 2007/08 (Choi et al. 2010).
Climatological T plots of and for four select grid cells (Fig. 9) provide further insights into future changes in ΔFD, ΔST, and ΔAF. The grids were selected along similar parallels and meridians and show well the thaw–freeze timing shifts against the climatological T variations. As previously mentioned, drawn from the NARCCAP ensemble mean T tends to be cooler than NARR and WFD. This bias and a lower spatial variance in the model simulated fields leads to a smoother pattern in ΔFD, with more modest reductions projected over the Intermountain West (Fig. 3a). Comparing daily T variations for a representative grid cell in Oregon (Fig. 9c) with those for grid cells in Nebraska, North Dakota, and Alberta (Figs. 9d, 9b, and 9a, respectively) shows how the rate of change in climatological daily T (expressed by ) are much lower over areas of the western United States relative to other areas (Fig. 9).
Air temperature variations near the time of spring thaw are influenced by modifications due to a number of factors including water phase transitions involving, for example, snowmelt and evaporation. The rate of change in sun angle may also play a role. Put simply, areas where temperatures rise relatively slowly in spring will experience a greater advance in the timing of spring thaw for any given amount of future warming. The seasonality of surface feedbacks involving land snow-cover disappearance is expected to cause warming at differing rates. As mentioned above, freezing days will decrease most in areas where mean annual temperature is close to the freezing point. Snow–albedo feedbacks, which should be expressed in the models, are another obvious influence. Northern areas that do not lose much snow cover are far less susceptible to strong feedback-induced warming, and thus large FD declines. In their study using a coupled GCM, Meehl et al. (2004) suggested that patterns of the future change in frost days depended partly on regional changes in circulation. However, we note considerable decreases in FD across the western United States under the uniform 2.5°C warming (Fig. 3d). Interestingly, over much of western Canada warming captured by ensemble mean T is greater in autumn compared to spring (Fig. 1). This differential warming may explain ΔAF exceeding ΔST across much of Canada. Sea ice losses simulated in the driving GCMs represent a potential influence on future autumn warming across the high latitudes (Serreze et al. 2009).
The connection between regional patterns in ΔFD and climatological T variations near the time of spring thaw are further illustrated by applying a uniform 2.5°C warming to observed normal (1971–2000) data for a series of stations in a roughly north–south orientation in the northeastern United States. The daily climatological T for each station is shown in Fig. 10. Applying a uniform warming to each time series leads to ΔST exceeding ΔAF (Fig. 11). This comparison clearly shows how ΔST, ΔAF, and ΔFD are correlated with annual temperature and the rate of change in T near the time of spring thaw. In other words, each of the three metrics increases in magnitude moving from north to south through this region. Similarly, differences in explain the general north–south gradient in ΔFD over the central United States for the uniform warming case shown in Fig. 3d and as seen in Figs. 9b,d. The northeast station comparison shows shifts to earlier ST that exceed the delays in AF (Fig. 11).
4. Summary and conclusions
Daily 2-m T from the NARCCAP models was used to assess likely changes in seasonal temperatures, freezing days, and thaw timings by midcentury across North America. Averaging the daily 30-yr climatological T across the 11 GCM–RCM pairs minimizes the unwanted effects of individual model outliers, allowing us to elucidate the underlying change patterns that arise as a result of climatic warming. We expect that the analysis and results drawn from the NARCCAP, NARR, and WFD datasets used here represent an improvement over direct GCM outputs in two ways. First, dynamical downscaling by RCMs has been demonstrated to be more accurate in simulating climate than in GCMs, particularly in areas of strong topographic relief (Gao et al. 2011; Elguindi and Grundstein 2013). Second, errors in baseline climate model temperatures were ameliorated using the delta method of bias correction. The results of this study reflect the degree of warming projected by the RCMs forced with boundary conditions from the associated GCM. Averaging the model temperatures over the ensemble members helps to reveal the mean temperature response to rising GHGs that largely controls the magnitude and spatial pattern in ΔFD.
Our results suggest that large parts of the United States, particularly western areas, will experience considerable declines in freezing days. The frequency of freezing days will decrease by as much as 90 days across parts of the United States and coastal western Canada. The greatest decreases can be expected near the southern border of the part of the continent currently experiencing seasonal freezing conditions. Changes generally will be greatest in areas of western North America with annual mean T between 2° and 6°C. Locally higher declines can be expected in regions where climatological daily temperatures tend to rise more gradually near the time of spring thaw. Across parts of the western and central United States spring thaw will advance more than autumn freeze will delay, with the opposite trend likely to occur across much of Canada. Climate model projections of temperature changes related to snow processes and associated snow–albedo feedbacks carry large uncertainties. By extension, it is reasonable to assume that confidence in future decreases in freezing day frequency across the higher elevations in western North America is lower there than over other parts of the continent. The frequency of freezing days represents just one of several useful temperature-based metrics that can be derived from climate model simulation outputs—for example, heating and cooling degree days or the number of extreme heat events. Continued development of models capable of accurate simulations of surface climate at higher spatial scales should improve projections of spatial variations and local impacts of future changes in seasonal transitions and freezing conditions.
We wish to thank the North American Regional Climate Change Assessment Program (NARCCAP) for providing the data used in this paper. NARCCAP is funded by the National Science Foundation (NSF), the U.S. Department of Energy (DOE), the National Oceanic and Atmospheric Administration (NOAA), and the U.S. Environmental Protection Agency (EPA) Office of Research and Development. We also thank the three anonymous reviewers for their constructive comments on the earlier version of the manuscript. RSB was supported by the Department of the Interior Northeast Climate Science Center, Grant G12AC00001 from the United States Geological Survey. The content of the paper is solely the responsibility of the authors and does not necessarily represent the views of the Northeast Climate Science Center or the USGS. This manuscript is submitted for publication with the understanding that the United States government is authorized to reproduce and distribute reprints for governmental purposes.