Abstract

The present study examines the effects of historic climate variability on cold-season processes, including soil temperature, frost depth, and the number of frost days and freeze–thaw cycles. Considering the importance of spatial and temporal variability in cold-season processes, the study was conducted in the midwestern United States using both observations and model simulations. Model simulations used the Variable Infiltration Capacity (VIC) land surface model (LSM) to reconstruct and to analyze changes in the long-term (i.e., 1917–2006) means of soil frost variables. The VIC model was calibrated using observed streamflow records and near-surface soil temperatures and then evaluated for streamflow, soil temperature, frost depth, and soil moisture before its application at the regional scale. Soil frost indicators—such as the number of frost days and freeze–thaw cycles—were determined from observed records and were tested for the presence of significant trends. Overall trends in extreme and mean seasonal soil temperature from 1967 onward indicated a warming of soil temperatures at a depth of 10 cm—specifically in northwest Indiana, north-central Illinois, and southeast Minnesota—leading to a reduction in the number of soil frost days. Model simulations indicated that by the late-century period (1977–2006), soil frost duration decreased by as much as 36 days compared to the midcentury period (1947–76). Spatial averages for the study area in warm years indicated shallower frost penetration by 15 cm and greater soil temperatures by about 3°C at 10-cm soil depth than in the cold years.

1. Introduction

Many recent studies of observation records have suggested that climate variability has increased and climate change has accelerated in the recent years. For instance, global and regional surface temperatures have increased in the twentieth century, with the greatest warming occurring in the last three decades (WMO 2006). The frequency of extreme temperature and precipitation events has also increased; specifically, the earth has experienced 11 of its top 12 warmest years in the last 12 yr (1995–2006) based on the instrumental record available since 1850 (Alley et al. 2007). In the Northern Hemisphere, there has been a significant reduction in mean snow cover area (e.g., Dye and Tucker 2003; Mote et al. 2005), while the maximum area covered by seasonally frozen ground has decreased by about 7% since the 1900s, and spring snowmelt has occurred earlier (Lemke et al. 2007; Stewart et al. 2005; Hodgkins and Dudley 2006). Easterling (2002) found that the number of days with subfreezing air temperature has decreased by four days per year in the United States during the period of 1948–99, while in recent decades the frost-free season based on air temperature has started 11 days earlier in the northeastern United States (Cooter and Leduc 1995). The air temperature frost-free season length has increased more in the western United States than in the eastern United States (Easterling 2002; Kunkel et al. 2004). While most of the studies were dedicated to studying the effect of climate variability at larger scales (Easterling 2002; Kunkel et al. 2004; Frauenfeld et al. 2007), fewer studies have concentrated on regional scales (Cooter and Leduc 1995; Kling et al. 2003). Cold-season processes are associated with a high degree of spatial (e.g., northern Minnesota versus southern Indiana) and temporal (e.g., onset–thaw of soil frost) variability in the midwestern United States, which underscores the need for a study dedicated to this region.

Climate variability and climate change may directly influence cold-season variables. For instance, winter soil temperature is directly governed by the increase in both air temperature and precipitation—mostly snowfall (Zhang et al. 2001, 2003). Hardy et al. (2001) suggested that the shifting of regional climate toward lower snow accumulations and shorter durations of snow on the ground will result in deeper soil freezing and more days with soil frost in temperate forests. Recently, an observational study found that decreasing snow cover has led to cooler soil temperatures in northern Indiana (Sinha and Cherkauer 2008). However, increases in frozen ground are short-term effects of limited spatial extent (Zhang 2005), as increasing air temperatures eventually overwhelm changes in insulation due to the snowpack. For instance, although winter air temperature increased by 4°–6°C at Irkutsk, Russia, during the past century, the soil temperature has increased by up to 9°C because of an increase in early winter snowfall and early snowmelt in spring (Zhang et al. 2003). Thus, changes in cold-season soil temperatures and its derived variables, such as the number of soil frost days, are influenced by changes in air temperature and winter precipitation as snow.

Historic observations in the midwestern United States indicate that climate variability has increased during the past century, which is likely to reinforce spatial–temporal variability in regional cold-season processes. The northern regions of the Midwest have warmed by more than 2°C, whereas the southern portions of the Ohio River Valley have cooled by about 0.5°C during the past century (Easterling and Karl 2001). This region has also experienced an increase in annual precipitation by up to 20% with about double the frequency of heavy precipitation events in the latter part of twentieth century than in the early part of the century (Easterling and Karl 2001; Kunkel et al. 1998). In addition, the midwestern United States experiences uncertainties in the frequency and severity of cold-season storms that affect flooding in late winter and early spring, the number of soil frost days and its time of occurrence, and soil temperatures. Therefore, we addressed the following questions in the present study: (i) How significant were the trends–changes associated with cold-season variables (soil frost, number of frost days and freeze–thaw cycles) in the midwestern United States? (ii) How has climate variability influenced regional cold-season variables in the past century? (iii) How sensitive were the responses of cold-season variables during the extreme climatic years (e.g., warm–cold and wet–dry)? We addressed these questions by integrating observations and modeling forced with long-term meteorological observations.

2. Background

To determine how climate variability affects soil frost, cold season descriptive variables were determined from both observed and simulated records. The selected simulated soil frost variables were calibrated and evaluated with observations before being used to extend the time series for the past century and the entire study area using the Variable Infiltration Capacity (VIC) model.

a. VIC model

The VIC model (Liang et al. 1994, 1996; Cherkauer and Lettenmaier 1999) is a large-scale hydrologic model that has been used for large watersheds (>10 000 km2) and has performed well under different climatic regimes (e.g., Liang et al. 1996; Wood et al. 1997; Lohmann et al. 1998a,b; Nijssen et al. 2001a,b; Cherkauer and Lettenmaier 2003; Bowling et al. 2003a,b). The details of the VIC model are described in Liang et al. (1994, 1996) and Cherkauer and Lettenmaier (1999), and only a brief description of the details is provided here. Distinguishing characteristics of the model include its representation of the spatial subgrid variability of infiltration through the Variable Infiltration Capacity curve, and it use of a nonlinear Arno parameterization for baseflow recession (Abdulla et al. 1999). The horizontal land surface is classified into any number of user-defined vegetation classes based upon the fraction of the grid cell it occupies. The model can be driven by daily meteorological variables, and it predicts surface temperature, evapotranspiration, runoff, snow water equivalent (SWE), and other variables for each grid cell. Runoff and baseflow can be routed to the basin outlet from each grid cell using a stand-alone routing model developed by Lohmann et al. (1998a,b).

Updates to the representation of cold-season processes in the VIC model are detailed in Cherkauer and Lettenmaier (1999) and Cherkauer et al. (2003), though a brief description is provided here for clarity. Cherkauer and Lettenmaier (1999) developed a frozen soil algorithm that solves for thermal fluxes at a user-defined number of soil nodes through the soil column and predicts soil water and ice content for each soil layer. In addition, the moisture fluxes are adjusted by the presence of ice. Calculating fluxes using both ice and liquid soil moisture makes the soil appear wet, which reduces infiltration, while drainage is reduced when only liquid soil moisture is used.

The model employs a two-layer snow algorithm, in which a thin surface layer is used to estimate energy exchange with atmosphere and the bottom layer acts as storage to simulate deeper snowpacks (Cherkauer and Lettenmaier 1999). The snow algorithm was originally designed for use within the Distributed Hydrology Soil Vegetation Model (DHSVM; Storck and Lettenmaier 1999). Cherkauer et al. (2003) added explicit representation of snow interception in the canopy to improve representation of snow accumulation and ablation in forested regions. Canopy intercepted snow is removed through sublimation, mass release, and meltwater drip. The model, as used for this study, assumed a spatially uniform snow cover within each vegetation type, which suggests an underestimation of snowmelt rates under thin and discontinuous snow. It also assumes a new snow density of 100 kg m−3. Densification of older snowpacks is based on Snow Thermal Model 89 (SNTHRM.89; Jordan 1991). The snowpack surface albedo curve parameters are based on snow age and season (Bras 1990). Maximum temperature when precipitation can occur as only snow was set to 0.25°C, whereas the minimum temperature when precipitation can fall completely as rain was defined as −0.25°C.

The VIC model version 4.1.0 r3 was used in this study with the finite difference soil thermal solution described by Cherkauer and Lettenmaier (1999) using a constant bottom boundary temperature. The thermal damping depth was set to 10 m, and the model was allowed to spin up for 2 yr.

b. Study area

The study area consists of the following six states in the midwestern United States: Minnesota (MN), Iowa (IA), Wisconsin (WI), Illinois (IL), Michigan (MI), and Indiana (IN; Fig. 1). These regions were selected because the spatial and temporal variability in seasonal soil frost is likely to vary with air temperature, which increases from north to south, and precipitation, which increases from the northwest to the southeast. This area comprises parts of the upper Mississippi River basin, the Ohio River basin, and the Great Lakes drainage basin. The area consists of the rolling forested landscapes of southern Illinois, and Indiana; a central region of the relatively flat farmland that produces mostly corn and soybean; and the evergreen and deciduous forests of northern Wisconsin, Minnesota, and Michigan (Easterling and Karl 2001). In winter, the region is influenced by cold Arctic air masses, which move into the region, resulting in frequent storms and windy conditions. The total annual precipitation varies from as low as 640 mm in the western Minnesota and Iowa to more than 1020 mm along the Ohio River. Annual average temperatures vary from 2.75°C in the north to 15°C in the south.

Fig. 1.

Map showing the locations of 1) selected sites that have been collecting soil temperatures: (a) MOR, (b) WAS, and (c) Lamberton (LAM) in MN; (d) Chatham (CHA) and (e) northwest Michigan (NMI) in MI; (f) Hancock, WI (HAN); (g) DEC, (h) AME, and (i) Burlington (BUR) in IA; (j) Freeport (FRE), (k) Urbana (URB), and (l) Orr (ORR) in IL; (m) WAN, (n) West Lafayette (WLA), and (o) Dubois (DUB) in IN; and 2) calibrated watersheds in the region: ANOKA–Mississippi River at Anoka, MN; MINNR–Minnesota River at Jordan, MN; IOWAR–Iowa River at Wapello, IA; CHIPR; GRAND–Grand River at Grand Rapids, MI; ILLIR–Illinois River at Valley City, IL; and WABAS. Stars represent streamflow gauging stations.

Fig. 1.

Map showing the locations of 1) selected sites that have been collecting soil temperatures: (a) MOR, (b) WAS, and (c) Lamberton (LAM) in MN; (d) Chatham (CHA) and (e) northwest Michigan (NMI) in MI; (f) Hancock, WI (HAN); (g) DEC, (h) AME, and (i) Burlington (BUR) in IA; (j) Freeport (FRE), (k) Urbana (URB), and (l) Orr (ORR) in IL; (m) WAN, (n) West Lafayette (WLA), and (o) Dubois (DUB) in IN; and 2) calibrated watersheds in the region: ANOKA–Mississippi River at Anoka, MN; MINNR–Minnesota River at Jordan, MN; IOWAR–Iowa River at Wapello, IA; CHIPR; GRAND–Grand River at Grand Rapids, MI; ILLIR–Illinois River at Valley City, IL; and WABAS. Stars represent streamflow gauging stations.

3. Methods and data

a. Observational data

1) Daily soil temperatures

Data for soil temperatures were obtained from sites that have collected at least 10 years of consistent soil temperature at a depth of 10 cm. The longest period for which consistent soil temperature data were available was 1966 to 2006. Data since 1982 were obtained from the National Climatic Data Center’s (NCDC) Summary of the Day (SOD) dataset (available online at http://www.ncdc.noaa.gov/oa/ncdc.html) and the Illinois State Water Survey (WARM 2008). Data prior to 1982 were available only as paper records, and they were digitized from the Climatological Data for five of the six states within the study domain (U.S. Department of Commerce 1966–1982). Maximum and minimum soil temperatures at a depth of 10 cm were taken using a Palmer dial-type temperature sensor, which has a precision of 0.5°C (Schaal et al. 1981). All the measurements for these sites were taken under a bare soil surface (refer to the Appendix A). Data for the Illinois sites were collected by the Illinois State Water Survey and represent soil temperatures at 10 cm below a grass surface.

2) Weekly frost depth and snow depth

Weekly frost and snow depths were obtained from the National Snow and Ice Data Center (NSIDC) for selected sites in Minnesota during the two cold seasons of 1975/76 and 1979/80 for comparison of observed and simulated variables. The frost depths were measured using frost tubes (Haugen and King 1998).

3) Monthly stream flow data

Monthly stream flow data were obtained since 1982 for selected basins monitored by the U.S. Geological Survey (USGS) based on criteria described in section 3e. These data were useful for model calibration and evaluation. Several gauged watersheds within the study domain were selected for calibration (Fig. 1).

b. Model input data

1) Soil data and vegetation data

Soil parameters were processed from a multilayer soil characteristics dataset at a 1-km resolution for the conterminous United States (CONUS-SOIL; Miller and White 1998). The CONUS-SOIL data are based on the original State Soil Geographic database (STATSGO). Land use parameters were acquired from the Land Data Assimilation System (LDAS) project (Mitchell et al. 2004; Maurer et al. 2002). Additional parameters required for the frozen soil algorithm were taken from Mao and Cherkauer (2009). The land use map was processed using vegetation classifications based on 1-km resolution satellite data from 1992 to 1993 (Hansen et al. 2000). These data have been gridded to a resolution of ⅛th-degree latitude and longitude for the conterminous United States. Monthly leaf area index (LAI) data for each vegetation type and each ⅛° grid cell have been obtained from Myneni et al. (1997). Vegetation library parameters were obtained from Mao and Cherkauer (2009).

2) Extension of daily gridded meteorological data

Long-term records of daily precipitation and maximum and minimum air temperature from 1915 to 2006 were obtained from the NCDC’s SOD dataset and gridded at a spatial resolution of ⅛° latitude by longitude (about 14 km × 11 km) for the midwestern United States using the methodology of Hamlet and Lettenmaier (2005). This gridding methodology is applicable for long-term trend analysis, as it accounts for changes in meteorological station location and temporal inconsistencies, which otherwise can lead to spurious trends in precipitation and temperature. Hamlet and Lettenmaier (2005) found that despite acceptable streamflow simulations through the VIC model for recent decades, streamflow prior to 1950 was strongly biased in comparison with that after 1950 due to temporal heterogeneities in the driving data. Hamlet and Lettenmaier (2005) suggested a methodology to make temporal corrections based upon monthly U.S. Historical Climatology Network (HCN) records and topographic corrections based upon monthly precipitation maps produced by the Parameter-elevation Regressions on Independent Slopes Model (PRISM). This methodology has been tested and verified for the state of California and has now been applied to the midwestern United States starting in 1915 for this study. Daily wind speeds since 1949 were obtained from the National Center for Atmospheric Research–National Centers for Environment Prediction (NCAR–NCEP) reanalysis project (Kalnay et al. 1996), whereas the daily wind speed prior to 1948 was derived from post-1949 wind climatology as described in Hamlet and Lettenmaier (2005).

c. Data processing

1) Soil temperature

Daily observed soil temperature data from selected stations were subjected to quality-checking procedures to identify missing dates, inconsistencies, and unreasonably high diurnal ranges in soil temperatures (Sinha and Cherkauer 2008). The erroneous data were replaced with “no data” values for making comparisons. Years missing more than 10% of observed values during the cold season were eliminated to avoid any biases in trend analysis. Cold-season statistics were computed for a period from September through May, which incorporates the earliest observed occurrence of soil frost and the latest observed thaw at all the sites in the region. The following cold-season variables were estimated from the daily time series of maximum and minimum soil temperatures (Tsoil) at a depth of 10 cm as described in Sinha and Cherkauer (2008):

  • (i) Seasonal mean maximum Tsoil is determined by calculating the mean of daily maximum Tsoil data series during the cold season.

  • (ii) Seasonal mean minimum Tsoil is determined by calculating the mean of daily minimum Tsoil data series during the cold season.

  • (iii) Extreme minimum temperature is determined by calculating the extreme value of daily minimum Tsoil data series during the cold season.

  • (iv) Annual soil frost days are estimated by counting the number of days when the soils are frozen during the cold season.

  • (v) Annual freeze–thaw cycles are estimated by counting the number of times soil temperature changes between frozen and thawed states during the cold season. Details provided in Sinha and Cherkauer (2008) are summarized for clarity. A freeze event occurs when at least daily minimum Tsoil on the first day is at or below 0°C and both daily maximum and minimum Tsoil on the subsequent day are also at or below 0°C. Similarly, a thaw event occurs following a freeze event when both maximum and at least one minimum Tsoil on two consecutive days are greater than 0°C. A freeze event, followed by a thaw event at a given soil depth is defined as one freeze–thaw cycle. A threshold value of −0.25°C (instead of 0°C in the case of observed data) was used to estimate freeze–thaw cycles from model-simulated daily maximum and minimum soil temperatures. This step is useful to make similar comparisons from frost variables computed from observed soil temperatures. The observed soil temperature sensor has a precision of 0.5°C, which implies that any fluctuation of less than 0.25° around 0°C will not be measured by the sensor.

  • (vi) Onset day of soil frost is determined by computing the first day since 1 September when soil is frozen.

  • (vii) Last thaw day is determined by computing the last day of soil frost since 1 September.

2) Frost depth

Among the stations measuring frost depth in the study region, only data from the Minnesota sites passed quality checks for consistency and continuity and were included in this analysis.

3) Climate variability in the gridded meteorological forcing data

Climate variability for the study area was assessed by estimating extreme temperature and precipitation events during the cold season derived from the daily gridded meteorological forcing data for precipitation and air temperature. Mean air temperatures were calculated by averaging daily maximum and minimum air temperature over the cold season, from September through May, for the entire region. Seasonal anomalies were normalized for both precipitation and average air temperature (Fig. 2) using

 
formula

where Xnormalized is the normalized time series, Xt is the time series for the cold season, Xmean is the mean, and Xmax is the maximum value and Xmin is the minimum value of the time series. A threshold value of 0.3 was selected to define multiple extreme events that occurred during the past century to distinguish into warm and wet years based upon air temperature and precipitation, respectively (Fig. 2). This means that the years with values greater than 0.3 were considered as warm based on average air temperature and considered as wet based upon precipitation. Similarly, a threshold of −0.3 was selected to identify cold and dry years. Only one year—1917—met the criteria of a cold year, whereas other extreme events were spread over multiple years. Seasonal statistics were analyzed to study the effects of climate variability on selected soil frost indicators.

Fig. 2.

Cold-season normalized anomalies of gridded meteorological forcing data from 1917 to 2006 averaged over the entire study area for (a) precipitation (PRCP) and (b) mean air temperature. The continuous line indicates a significant trend; the dashed line indicates a nonsignificant trend. The Mann–Kendall slope of trend line is indicated.

Fig. 2.

Cold-season normalized anomalies of gridded meteorological forcing data from 1917 to 2006 averaged over the entire study area for (a) precipitation (PRCP) and (b) mean air temperature. The continuous line indicates a significant trend; the dashed line indicates a nonsignificant trend. The Mann–Kendall slope of trend line is indicated.

d. Observed data analysis

Data analysis was performed in two stages: 1) data series were tested to identify the presence of any significant autocorrelations, which is the correlation of a variable with itself over successive time intervals, and 2) data series were tested to identify significant trends in soil frost variables using the Mann–Kendall test. The presence of autocorrelation increases the possibility of detecting trends even if they are absent (Rao et al. 2003). In the case of significant autocorrelation in a data series, the effective number of independent observations decreases and thus accounting for it reduces the chances of detecting spurious trends.

Data series with more than 10 years of observations remaining, after data quality checks, were tested for significant trends using the Mann–Kendall test, as it is applicable for 10 or more observations. This test has been widely used for trend analysis in hydrologic studies and was applied in a similar form used by Hirsch et al. (1992). The details of the test were described in Sinha and Cherkauer (2008). In cases of significant autocorrelation, a modified Mann–Kendall test (Hamed and Rao 1998) was used that calculates the autocorrelation between the ranks of the data after removing the apparent trend, as described in Sinha and Cherkauer (2008). The trend analysis was performed for the soil frost variables described in section 3c for the following periods to maximize the number of stations measuring consistent soil temperature at a depth of 10 cm: (i) 1967–2006, (ii) 1974–2006, (iii) 1983–2006, and (iv) 1991–2006.

e. Modeling Strategy

The VIC model was calibrated using observed streamflow from 1982 to 1992 and seasonal average daily soil temperature at a depth of 10 cm. The model was evaluated for streamflow from 1993 to 2003 as well as average seasonal soil temperature, soil moisture (three sites in Illinois), and frost depth (three sites in Minnesota) based upon data availability. Direct comparisons of trends in model-derived soil frost variables (such as soil frost days and extreme seasonal daily Tsoil) have not been made with those estimated from the observed data at selected sites, as trends are difficult to identify in model-derived variables. However, the calibrated model was used to compare the general patterns of selected soil frost variables using Spearman rank correlation coefficient r before it was used to quantify long-term changes in cold-season parameters for the study domain.

1) Calibration and evaluation of streamflow

Seven major watersheds were selected that were minimally influenced by major dams and reservoirs (Fig. 1). Since land use was based upon 1992/93 satellite data, the 10 years prior to this period (1982–91) were selected as the calibration period, whereas the subsequent 10 years (1994–2003) were selected as the evaluation period. This minimized the effects on streamflow due to land use change. The VIC model was calibrated on a watershed basis using observed streamflow measurements. Calibration was performed to improve the agreement of observed and simulated discharge volumes and the shape of hydrographs. A single set of model parameters that control infiltration, runoff, and baseflow were obtained for the entire study domain, to minimize the effect of parameters change between watersheds on any spatial analysis. Figure 3 shows the mean monthly hydrographs for the seven watersheds for the entire length of the calibration and validation periods.

Fig. 3.

Comparison of monthly observed (Obs) and simulated (Sim) streamflow from 1982 to 2003 at the seven calibrated watersheds. Shown are E (%) and d (%) between observed and simulated streamflow over the entire period.

Fig. 3.

Comparison of monthly observed (Obs) and simulated (Sim) streamflow from 1982 to 2003 at the seven calibrated watersheds. Shown are E (%) and d (%) between observed and simulated streamflow over the entire period.

The following accuracy indicators were used to evaluate the model performance:

The E and the d were applied to monthly streamflows for the period 1982–2003 (Fig. 3). Values of E varied from 0.27 for the Chippewa River at Durand, WI (CHIPR), to 0.84 for the Wabash River at Mt. Carmel, IN (WABAS), and the values of d varied between 0.81 and 0.96 for the CHIPR and WABAS, respectively.

2) Calibration and Evaluation of soil temperatures

Fifteen sites, with soil temperature observations at 10 cm soil depth, were selected to calibrate and evaluate model performance. The key soil parameters—such as Ds (which is the fraction of maximum velocity of baseflow Dsmax from the lowest soil layer where nonlinear baseflow begins), parameter describing the variation of saturated hydraulic conductivity with soil moisture in the three soil layers, and average soil temperature used as the bottom boundary for soil heat flux solutions—were adjusted during the calibration process. The Spearman rank correlation coefficient r was selected to compare relationships between simulated and observed soil frost variables, because it does not require any assumption about the frequency distribution of the variables; instead, it uses the ranks of the data to estimate the correlation between the variables. It is computed using

 
formula

where R(x) and R(y) are the ranks of a pair of variables (x and y), with each containing n observations.

4. Results and discussion

a. Observational data analysis for soil frost variables

Trends in observed soil temperature records were assessed for four periods to maximize the length of record, as well as the number of active stations, since the majority of stations came online in the 1980s and 1990s. Trend analysis indicated that a higher percentage of stations show statistically significant trends in soil frost variables for the period of 1967 to 2006 than for any other period (Table 1). However, this period also had the fewest available stations. The percentage of stations showing significant trends in soil frost variables varied over the different periods considered. Most of the significant trends were observed in extreme minimum Tsoil and seasonal mean maximum Tsoil for all the periods under consideration. Interestingly, among significant trends, the magnitude of the Mann–Kendall slope for extreme minimum Tsoil was higher than the corresponding slopes in mean seasonal Tsoil for all the periods (Tables 2 –5), indicating greater changes in extreme minimum Tsoil. The greatest increase in extreme minimum Tsoil of 0.68°C since 1983 as well as the largest decrease in mean maximum Tsoil (−0.52°C since 1991) were observed at Morris, MN (MOR; Tables 4 and 5), indicating that variability in cold-season soil temperatures has decreased. This site also experienced the greatest rate of increase in annual number of soil frost days from 1991 to 2006.

Table 1.

Percentage of stations with statistically significant trends in soil frost variables for different periods at a significance level of 5%. The number of stations available for each period is provided in parentheses.

Percentage of stations with statistically significant trends in soil frost variables for different periods at a significance level of 5%. The number of stations available for each period is provided in parentheses.
Percentage of stations with statistically significant trends in soil frost variables for different periods at a significance level of 5%. The number of stations available for each period is provided in parentheses.
Table 2.

Stations showing significant trends in soil frost variables from 1967 to 2006. The values indicate the Mann–Kendall slope. Numbers in parentheses indicate the p value tested for significance; p < 0.05 is statistically significant and N indicates the absence of significant trend at a 5% significance level.

Stations showing significant trends in soil frost variables from 1967 to 2006. The values indicate the Mann–Kendall slope. Numbers in parentheses indicate the p value tested for significance; p < 0.05 is statistically significant and N indicates the absence of significant trend at a 5% significance level.
Stations showing significant trends in soil frost variables from 1967 to 2006. The values indicate the Mann–Kendall slope. Numbers in parentheses indicate the p value tested for significance; p < 0.05 is statistically significant and N indicates the absence of significant trend at a 5% significance level.
Table 5.

Same as Table 3 but for 1991–2006.

Same as Table 3 but for 1991–2006.
Same as Table 3 but for 1991–2006.
Table 4.

Same as Table 3, but for 1983–2006.

Same as Table 3, but for 1983–2006.
Same as Table 3, but for 1983–2006.
Table 3.

Stations showing significant trends in soil frost variables from 1974 to 2006. The values indicate the Mann–Kendall slope, and the number in parentheses indicate the p value tested for significance; N same as in Table 2.

Stations showing significant trends in soil frost variables from 1974 to 2006. The values indicate the Mann–Kendall slope, and the number in parentheses indicate the p value tested for significance; N same as in Table 2.
Stations showing significant trends in soil frost variables from 1974 to 2006. The values indicate the Mann–Kendall slope, and the number in parentheses indicate the p value tested for significance; N same as in Table 2.

During the periods since 1967 and 1974, trends that were statistically significant indicated a warming in extreme and mean seasonal Tsoil—specifically in northwest Indiana, north-central Illinois, and southeast Minnesota—leading to a reduction in the number of soil frost days [Figs. 4(1a)–4(4b)]. However, periods since 1983 and 1991 indicated a decrease in soil temperatures at a limited number of sites, such as Morris, where soil frost days increased with time. Significant negative trends in soil temperatures were identified only for the shorter duration records. For instance, Chatham, MI (CHA), and Farmland, IN (FRM), have experienced a reduction in extreme minimum Tsoil, with Farmland, IN, also experiencing a reduction in seasonal mean Tsoil from 1983 to 2006. The mixture of statistically significant positive and negative trends was indicative of the dominance of localized controls over climate signals on seasonal soil frost variables for shorter periods (i.e., since 1983 and since 1991, respectively).

Fig. 4.

Observed trends in the cold-season variables showing (1) mean maximum Tsoil (°C), (2) mean minimum Tsoil (°C), (3) extreme minimum Tsoil (°C), (4) number of frost days, and (5) number of freeze–thaw (F/T) cycles for the following periods: (a) 1967–2006, (b) 1974–2006, (c) 1983–2006, and (d) 1991–2006. Symbols represent station locations. Filled triangles indicate that trends are significant at the 5% significance level and empty triangles indicate nonsignificant trends. Triangles show the direction of trends.

Fig. 4.

Observed trends in the cold-season variables showing (1) mean maximum Tsoil (°C), (2) mean minimum Tsoil (°C), (3) extreme minimum Tsoil (°C), (4) number of frost days, and (5) number of freeze–thaw (F/T) cycles for the following periods: (a) 1967–2006, (b) 1974–2006, (c) 1983–2006, and (d) 1991–2006. Symbols represent station locations. Filled triangles indicate that trends are significant at the 5% significance level and empty triangles indicate nonsignificant trends. Triangles show the direction of trends.

Higher slopes with significant trends since 1991 indicate that for this shorter period, soil frost variables have generally experienced more rapid change than in the past (Table 5). In a few cases, when the trends were not significant, the direction of change was sensitive to the period in consideration. For example, at Waseca, MN (WAS), the direction of change reversed for the number of soil frost days when computed since 1974 to that computed since 1983 [Figs. 4(4b) and 4(4c)]. The 1970s were characterized by cold winters in the midwestern United States, which may have caused greater soil frost days during early 1970s. Despite the increase in soil frost days since 1983, the overall direction of change since 1974 indicated a decrease at WAS, where the changes were statistically insignificant.

b. Comparison of observations with simulated variables

Calibration and evaluation of the VIC model typically focuses on streamflow statistics; however, to use the model for analysis of soil frost processes, additional evaluations of soil temperature and frost depth were performed. Fifteen sites were selected for evaluation of observed and simulated cold-season variables (Fig. 1). For the evaluations, all sites, except for the three in Iowa, made use of daily maximum and minimum Tsoil data to compute seasonal statistics. The Iowa sites recorded only a single daily temperature measurement taken at the time of observation; hence, these sites were unable to capture the daily extreme values observed by the other sites. This implies that the seasonal extremes and frost dynamics were underrepresented at those sites.

Cold-season mean daily Tsoil comparisons indicated that the VIC model generally captured overall patterns of the observed daily Tsoil at most of the selected sites (Fig. 5); however, there were differences in absolute values. Here r between the observed and simulated mean daily Tsoil ranged from 0.75 for Dubois, IN (DUB), to −0.27 for MOR, with 7 out of the 12 sites having r ≥ 0.5 (excluding the three IA sites because of different observed Tsoil format). The general patterns in soil temperature were captured by the model at most sites, except at MOR and Decorah, IA (DEC), where the Spearman rank correlations were negative. Since the Iowa sites had a different format for the observed Tsoil, the simulated mean daily Tsoil were expected to be colder than observed Tsoil recorded at a specific time of a day. The Morris, MN, site, on the other hand, experienced the highest variability in observed Tsoil among all the sites, as described earlier, resulting in larger differences in simulated and observed Tsoil. The discrepancy in absolute values between observed and simulated Tsoil may be due to differences in meteorological data and soil data at the selected sites. The regional simulations were conducted using the ⅛th-degree gridded meteorological forcing data, which were spatially interpolated to the center of the grid cell using neighboring meteorological stations and may differ from the actual observations of daily precipitation and air temperatures at selected point locations. The model further disaggregated daily meteorological data to subdaily scale, which affected snowfall timings. In addition, the model assumes a maximum (minimum) temperature of 0.25° (−0.25°C) when precipitation falls as snow (completely as rain). This may cause differences between observed and simulated snow depth, resulting in changes in cold-season daily average Tsoil. Furthermore, the VIC model simulations make use of average soil conditions for the grid cell, which may not capture site-specific conditions despite forcing the VIC grid cell with the same vegetation type as the observational site. Overall, the VIC-simulated daily Tsoil were deemed acceptable, as the observed patterns were simulated at more than half of the selected sites and thus the model was used as an analysis tool to estimate long-term changes in soil frost variables.

Fig. 5.

Cold-season time series of observed and simulated mean daily soil temperatures at a depth of 10 cm for the following sites: (a) MOR, (b) LAM, (c) WAS, (d) CHA, (e) NMI, (f) HAN, (g) DEC, (h) AME, (i) BUR, (j) FRE, (k) URB, (l) ORR, (m) WAN, (n) WLA, and (o) DUB. Observations at the three sites in IA (DEC, AME, and BUR), unlike the other states, are for a specific time during a day and do not provide information about daily extremes.

Fig. 5.

Cold-season time series of observed and simulated mean daily soil temperatures at a depth of 10 cm for the following sites: (a) MOR, (b) LAM, (c) WAS, (d) CHA, (e) NMI, (f) HAN, (g) DEC, (h) AME, (i) BUR, (j) FRE, (k) URB, (l) ORR, (m) WAN, (n) WLA, and (o) DUB. Observations at the three sites in IA (DEC, AME, and BUR), unlike the other states, are for a specific time during a day and do not provide information about daily extremes.

Simulated snow depths and soil frost depths were evaluated with observations at three sites in Minnesota using observational data for the cold seasons of 1975/76 and 1979/80 (Fig. 6). Snow depth and frost depth were well simulated for 1979/80 at the Morris and Waseca sites, whereas the underestimation of snow depth at Lamberton resulted in deeper soil frost penetration [Figs. 6(1c)–6(3d)]. Overestimation of snow depth at Waseca for the 1975/76 season led to an underestimation of frost depth, because of an increase in the insulation effect of snow at the ground surface. A smaller overestimation of snow cover at the Lamberton site during 1975/76 did not have a significant effect on soil frost penetration, which was slightly deeper than observed [Fig. 6(2b)]. Since the distribution of precipitation between rain and snow relative to air temperature is dependent on geographical locations, snow depth is difficult to simulate (Fassnacht and Soulis 2002). Additionally, the VIC model is reconstructing subdaily precipitation and air temperature from daily values averaged for an ⅛th-degree grid cell, which further affects the occurrence of snow. With these caveats in place, the simulation of snow and soil frost depths by the calibrated VIC model were deemed satisfactory.

Fig. 6.

Daily time series of observed (red) and simulated (black) variables (cm) snow (snow depth) and frost depth (F. Depth) and observed thaw depth (blue, in cm) at the following sites in MN: (1) MOR, (2) LAM, and (3) WAS for the two cold seasons of 1975/76 and 1979/80.

Fig. 6.

Daily time series of observed (red) and simulated (black) variables (cm) snow (snow depth) and frost depth (F. Depth) and observed thaw depth (blue, in cm) at the following sites in MN: (1) MOR, (2) LAM, and (3) WAS for the two cold seasons of 1975/76 and 1979/80.

As a final evaluation, the calibrated model was used to estimate the same soil frost indicators computed from observed soil temperature records. Spearman correlation coefficients between the annual numbers of soil frost days estimated from observations and simulations ranged from −0.05 at Wanatah, IN (WAN), to 0.86 at Ames, IA (AME), with 5 out of 15 sites indicating r ≥ 0.50 (Fig. 7). Three sites had r ≤ 0, indicating significant differences between observed and simulated soil frost days, most likely related to the non-site-specific meteorology and calibration mentioned previously. Generally, the model predicted a higher number of soil frost days than observations. In contrast, there were fewer fluctuations in freeze–thaw events estimated from simulated Tsoil compared to those estimated from observed Tsoil. This resulted in r ≤ 0 for 4 out of 12 sites (excluding IA sites) for the number of freeze–thaw cycles (Fig. 8). It is harder to predict freeze–thaw cycles than the number of soil frost days, as soil temperature sensors may not observe small fluctuations around 0°C because of their limited precision and unknown accuracy, especially around 0°C.

Fig. 7.

Same as Fig. 5, but for annual time series of number of soil frost days estimated from observed and simulated soil temperatures at a depth of 10 cm.

Fig. 7.

Same as Fig. 5, but for annual time series of number of soil frost days estimated from observed and simulated soil temperatures at a depth of 10 cm.

Fig. 8.

Same as Fig. 5, but for annual time series of number of F/T cycles estimated from observed and simulated soil temperatures at a depth of 10 cm.

Fig. 8.

Same as Fig. 5, but for annual time series of number of F/T cycles estimated from observed and simulated soil temperatures at a depth of 10 cm.

To study the effects of climate variability on temporal and spatial variability in cold-season processes, simulations for the entire period (since 1917) as well as observed daily soil temperature were sorted into cold and warm and wet and dry years as described in section 3c. Because only one year (1917) met the criteria of cold year, and it did not overlap with observations, years 1978 and 1979 were used to represent cold years for comparisons. The differences between observed and simulated Tsoil were within 15% for at least half the sites for cold, wet, and dry climatic divisions (Table 6), whereas the model did not performed as well in warm years. Differences for warm years were still within 20% for five out of eight sites; therefore, the model performance was deemed acceptable to capture the long-term means in daily Tsoil for all climatic divisions and was used as an analysis tool to quantify climate variability effects.

Table 6.

Obs and Sim average cold-season daily soil temperature during extreme climatic years (warm, cold, wet, and dry) for sites measuring soil temperatures (°C) since 1970. Values in bold indicate that the difference between Obs and Sim soil temperatures were within 15% of each other.

Obs and Sim average cold-season daily soil temperature during extreme climatic years (warm, cold, wet, and dry) for sites measuring soil temperatures (°C) since 1970. Values in bold indicate that the difference between Obs and Sim soil temperatures were within 15% of each other.
Obs and Sim average cold-season daily soil temperature during extreme climatic years (warm, cold, wet, and dry) for sites measuring soil temperatures (°C) since 1970. Values in bold indicate that the difference between Obs and Sim soil temperatures were within 15% of each other.

c. Cold-season sensitivity to snow cover and air temperature under extreme conditions

To understand the role of snow cover and air temperature on soil temperature and soil frost development, observations of air temperature and precipitation since 1917 were used to sort simulated variables into warm, cold, wet, and dry years as described in section 3c. Time series of a representative year for each climate type were compared with air temperature for four locations across the study domain (Fig. 9). For all sites and climate conditions, once the soil was frozen, daily changes in soil temperature were smaller than for air temperature; when snow exceeded a threshold depth [∼(5–10 cm)], changes in daily near-surface soil temperature became very small. This is in agreement with the findings of Nikol’skii et al. (2002), who found that the top 10 cm of snow cover had a strong insulation effect on soil surface temperature. The minimum value of soil temperature is strongly correlated to the timing of snow accumulation. When snow accumulates soon after the occurrence of freezing air temperatures [e.g., Figs. 9(2c) and 9(2d)] and remains throughout the winter, minimum soil temperatures are close to 0°C. If snow accumulation occurs later, then soil temperatures will be colder even under warm [e.g., Fig. 9(2a)] or wet [Fig. 9(3c)] climate conditions. When little to no snow is present during the winter, air temperature is the dominant control on soil temperature and freezing [e.g., Fig. 9(1d)]; however, snow cover clearly is the dominant controlling factor, as the colder temperatures during a wet winter [Fig. 9(1c)] do not result in colder soils than a slightly warmer winter without snow [Fig. 9(1d)].

Fig. 9.

Daily values of observed air temperature (°C), simulated soil temperature (°C) and simulated snow depth (cm) in (1) northwest MN (48.8125°N, 96.1875°W), (2) northern WI (46.0625°N, 91.0625°W), (3) southern IA (41.5625°N, 94.0625°W), and (4) southern IN (38.8125°N, 86.6875°W), and the following cold seasons: (a) warm, (b) cold, (c) wet, and (d) dry. Simulated soil temperatures and snow depth were used because of the lack of observations at the selected sites during extreme climatic years.

Fig. 9.

Daily values of observed air temperature (°C), simulated soil temperature (°C) and simulated snow depth (cm) in (1) northwest MN (48.8125°N, 96.1875°W), (2) northern WI (46.0625°N, 91.0625°W), (3) southern IA (41.5625°N, 94.0625°W), and (4) southern IN (38.8125°N, 86.6875°W), and the following cold seasons: (a) warm, (b) cold, (c) wet, and (d) dry. Simulated soil temperatures and snow depth were used because of the lack of observations at the selected sites during extreme climatic years.

d. Spatial and temporal analysis

Subsequent to calibration and evaluation, the VIC model was used to extend the spatial and temporal analysis of soil frost. Model simulations starting in 1917 indicated that the annual number of soil frost days decreased significantly from north to south in the region [Fig. 10(1a)], with the average number of soil frost days ranging from 160 in northern Minnesota to about 4 days in southern Illinois and Indiana. The average number of freeze–thaw cycles at a depth of 10 cm was highest in northern Indiana and lower Michigan, the central parts of Wisconsin and Illinois, and in southern Iowa [Fig. 10(2a)]. The central regions of Wisconsin and lower Michigan were under agricultural or mixed use conditions based upon the land use data used for this study. Agricultural lands are more susceptible to soil frost than forests because without a canopy, agricultural fields are more exposed to cold air and surface residues. The lack of a canopy does lead to greater accumulation of snow in agricultural fields; however, that snow is more likely to be redistributed by wind, and it is fully exposed to the sun so that it melts more quickly. Thus agricultural fields were more exposed and left with less insulation for longer periods of time than forested environments, which resulted in higher numbers of freeze–thaw cycles. In contrast, western Minnesota and northern Wisconsin were mostly forested, which reduced the depth of ground snow but maintained it over longer periods of time. Snow cover insulates the soil surface from changes in air temperature and thus reduces the frequency of freeze–thaw cycles. Furthermore, in the northern regions of the study domain, the presence of more persistent below-freezing temperatures resulted in fewer freeze–thaw cycles. The number of freeze–thaw cycles was also limited in southern Illinois and Indiana because the warmer air temperatures resulted in fewer days with soil frost and thus fewer opportunities to complete freeze–thaw cycles.

Fig. 10.

Spatial plot indicating (a) annual average values since 1917, (b) late-century (1977–2006) minus early-century conditions (1917–46), and (c) late-century minus midcentury conditions (1947–76) for the following variables: (1) frost days (number), (2) F/T cycles (number), (3) onset day of soil frost, and (4) last thaw day.

Fig. 10.

Spatial plot indicating (a) annual average values since 1917, (b) late-century (1977–2006) minus early-century conditions (1917–46), and (c) late-century minus midcentury conditions (1947–76) for the following variables: (1) frost days (number), (2) F/T cycles (number), (3) onset day of soil frost, and (4) last thaw day.

The timing of the onset of soil frost and last spring thaw both varied across the study area. The earliest occurrence of soil frost in the region, on an average, took place around early November in northern Minnesota, whereas in the southern regions frost occurred around early to mid-December [Fig. 10(3a)]. The last thaw in spring occurred around early February in the southern regions and around the end of April in the northern regions [Fig. 10(4a)].

To study temporal variability, changes in soil frost variables were calculated between 30-yr groups representing early- (1917–46), mid- (1947–76) and late-century (1977–2006) periods. By studying the differences during early-century and midcentury periods with respect to the late-century period, systematic biases—such as the effects of annual variations that are not generally well represented in the model-simulated variables—were removed. Changes in the number of soil frost days in the late century with respect to the middle and early century indicated a decrease in soil frost days in most of the northern and central regions while an increase was observed in southern Illinois [Fig. 10(1b)] and the upper peninsula of Michigan [Fig. 10(1c)]. The decrease in soil frost was by as much as 24 days in lower Michigan and 18 days in southern Minnesota, southern Wisconsin, and the northern regions of Illinois and Indiana [Fig. 10(1c)]. Most of these regions are agricultural lands, which are more exposed to changes in air temperature. Despite these reductions, the number of freeze–thaw cycles was unaffected in most of the region in the late century as compared to the early-century period [Fig. 10(2b)]. In comparison, by midcentury, regions in the north indicated an increase or no change in freeze–thaw cycles while regions in the south indicated a decrease [Fig. 10(2c)]. Increases in soil temperatures may have resulted in more fluctuations around 0°C in the north while providing fewer opportunities for soil frost development in the south.

In the late-century period, fewer sites experienced an earlier onset of frost in central and southern Illinois, north-central Minnesota, and in the upper peninsula of Michigan [Fig. 10(3c)]. During this period, both the first occurrence of frost and the last thaw changed by as much as 18 days [Fig. 10(4c)]. Interestingly, the area that displayed a delay in the onset day of soil frost between the late- and early-century periods was larger than the region that experienced the same change between the late- and midcentury periods [Figs. 10(3b) and 10(3c)]. These sites were generally those where the number of soil frost days decreased [Figs. 10(1b) and 10(1c)]. On the other hand, the last day of thawing in the late-century period occurred earlier in most of the region [Figs. 10(4b) and 10(4c)], by as much as 27 days in lower Michigan and southern Illinois. This indicates a reduction in soil frost duration between the mid- and late-century periods by as much as 36 days—specifically in southeast Minnesota, northeast Iowa, and north-central Indiana.

Simulations from 1917 indicate that average winter soil temperatures at a depth of 10 cm for the months of December, January, and February ranged from −8°C in the north to +2°C in the south [Fig. 11(1a)]. In spring (March–May), the monthly average temperature varied from about 1° to about 12°C [Fig. 11(2a)]. Frost depths were greatest in the northern regions and increased to their maximum values in the spring [Figs. 11(3a) and 11(4a)]. This may be due to the lower-than-average SWE during spring in comparison to the winter season [Figs. 11(5a) and 11(6a)]. In winter, higher average monthly SWE was indicative of more days with snow cover and deeper snow, which further insulated the soil surface and decreased soil frost penetration. Furthermore, the development of soil frost continued from winter into spring, leading to increased frost depths.

Fig. 11.

Spatial plots indicating (a) seasonal average since 1917, (b) late-century (1977–2006) minus early-century conditions (1917–46), and (c) late-century minus midcentury conditions (1947–76) in Tsoil (°C), frost depth (cm), and SWE (cm). The odd numbers represent winter (December–February) and even numbers represent spring (March–May).

Fig. 11.

Spatial plots indicating (a) seasonal average since 1917, (b) late-century (1977–2006) minus early-century conditions (1917–46), and (c) late-century minus midcentury conditions (1947–76) in Tsoil (°C), frost depth (cm), and SWE (cm). The odd numbers represent winter (December–February) and even numbers represent spring (March–May).

Winter soil temperatures increased in the northern part of the study domain with respect to the early-century period, while the central part of the study domain experienced warming between the middle and late century [Figs. 11(1b) and 11(1c)]. In contrast, northeastern Minnesota and the upper peninsula of Michigan underwent a cooling trend in winter soil temperatures during the latter part of the twentieth century. This may be due to a decrease in average monthly winter SWE, which was indicative of reduced snow cover in these regions of Minnesota and Michigan [Fig. 11(5c)]. This also led to increased frost depth during the winters of the late-century period with respect to midcentury [Fig. 11(3c)]. During spring, most of the northern regions experienced an increase in soil temperature, leading to reductions in frost depth in northern Minnesota and western Wisconsin [Figs. 11(4b) and 11(4c)]. Winter SWE has increased in southern-central Illinois and Indiana, southern Minnesota and northeast Iowa by up to 9% in the late century in comparison to the midcentury while spring SWE has decreased in the northern region, southeastern Wisconsin, and northern Michigan by 9% [Figs. 11(5c) and 11(6c)]. The decrease in spring SWE coincided with the regions that were classified as agricultural and mixed use. Agricultural regions have higher snow accumulations on the ground surface than forested sites because trees intercept snowfall, which reduces snow accumulation. Similarly, the regions where SWE has increased during winter were mostly under agricultural areas.

Cold-season soil temperatures in warm years varied from 0°C in the north to 12°C in the south, while in the cold years temperatures ranged from −5°C in the north up to 9°C in the south, an average decrease of 3°C (Fig. 12). Warm years had shallower frost penetration, which may even be shallower by up to 100 cm at a few sites in northern Minnesota in comparison to the cold years [Figs. 12(2a) and 12(2b)]. Cold-year soil frost was deeper, specifically in most of Minnesota and northern Wisconsin, with average depths for the entire region exceeding by 15 cm than those simulated for warm years. Cold years also distinguished themselves from warm years through higher snow accumulations. The variability in SWE during cold years was higher, with a greater number of sites experiencing extreme values of SWE, than in the warm years (Fig. 13). However, total precipitation in the cold years was lower than in warm years.

Fig. 12.

Cold seasonal averages for the following selected soil frost variables: 1) Tsoil (°C), 2) frost depth (cm), and 3) SWE (mm), 4) PRCP (mm), and 5) monthly air temperature Tair (°C) with the following years: (a) warm, (b) cold, (c) wet, and (d) dry.

Fig. 12.

Cold seasonal averages for the following selected soil frost variables: 1) Tsoil (°C), 2) frost depth (cm), and 3) SWE (mm), 4) PRCP (mm), and 5) monthly air temperature Tair (°C) with the following years: (a) warm, (b) cold, (c) wet, and (d) dry.

Fig. 13.

Box plots showing (top) frost depth, (middle) SWE, and (bottom) PRCP for cold seasons representing: warm, cold, wet, and dry years. The tops and bottoms of the boxes indicate the 25th and 75th quartiles, respectively; the lines in the middle represent the median value, the whiskers indicate the maximum and minimum values in the distribution, and crosses indicate values that fell outside of 1.5 × the interquartile range and were therefore classified as outliers.

Fig. 13.

Box plots showing (top) frost depth, (middle) SWE, and (bottom) PRCP for cold seasons representing: warm, cold, wet, and dry years. The tops and bottoms of the boxes indicate the 25th and 75th quartiles, respectively; the lines in the middle represent the median value, the whiskers indicate the maximum and minimum values in the distribution, and crosses indicate values that fell outside of 1.5 × the interquartile range and were therefore classified as outliers.

In wet years, soil temperatures were warmer than in the dry years (Figs. 12c and 12d), but the spatial average of soil temperatures over the entire study domain for both wet and dry years were similar (5.06° and 4.98°C, respectively). Wet years experienced higher SWE than dry years, resulting in reduced depth of frost penetration, as expected. The median precipitation in wet years was about 30% higher than that for dry years (Fig. 13). Although soil temperature changes generally followed patterns similar to those of air temperature during wet and dry years, regions with higher SWE, such as the upper peninsula of Michigan, experienced warmer soil temperatures than the corresponding air temperatures because of the insulation effect of snow cover. The frequency and severity of extreme temperature and precipitation events are likely to be enhanced even more in the future than has occurred in the past, resulting in more dramatic changes in cold-season processes.

5. Conclusions

This study focused on the effects of historic climate variability on cold-season variables in the midwestern United States, using both observations and model simulations. Soil frost indicators—such as the number of frost days and freeze–thaw cycles—were determined from observed records for different periods based upon the availability of data and were tested for the presence of significant trends. We used the VIC model to reconstruct long-term historic time series of cold-season variables and then analyzed trends and associated spatial and temporal variability. We also analyzed the hydrologic response to extreme warm, wet, cold, and dry years to understand the effect of climate variability on seasonal soil frost. The primary findings of this study are summarized as follows:

  • An analysis of observed data indicated that the percentage of stations showing significant trends and the direction of those trends varied for the different observational periods considered. A higher percentage of stations experienced statistically significant trends in soil frost variables since 1967 in comparison to the shorter periods considered in this study. Overall trends in extreme and mean seasonal Tsoil from 1967 and from 1974 to present have generally indicated a warming in soil temperatures, specifically at northwest Indiana, north-central Illinois, and southeast Minnesota, leading to a reduction in the number of soil frost days. The record is mixed on shorter time scales (since 1983 and since 1991), indicating localized effects dominate the regional climate signal in this region.

  • Analysis of daily time series of simulated soil temperature and snow depth at four sites indicated that snow is the dominant control on soil temperatures. When significant snow (depth > 5 cm) is present, soil temperatures are relatively constant and do not reach the extreme minimums of a soil directly exposed to cold air temperatures. Actual soil temperatures under snow are controlled by the timing of snow accumulation, such that soil temperatures can be substantially colder in years with late accumulation or discontinuous snow cover even when winter air temperatures are higher. With limited to no snow cover, soil temperatures are highly correlated to air temperatures.

  • The duration of soil frost has decreased by as much as 36 days in the late-twentieth-century period (1977–2006) as compared to the midcentury (1947–76) period in southeast Minnesota, northeast Iowa and north-central Indiana, increasing the length of growing season. In most of the central and southern regions, the length of the growing season increased by about two weeks in the late twentieth century as compared to the midcentury period. These changes in soil frost duration suggest a shift in the time of sowing and fertilizer applications earlier in the spring by two weeks. In contrast, soil frost days have increased by two weeks between the middle and end of the twentieth century in the upper peninsula of Michigan, where the onset of frost also now occurs two weeks earlier. The central regions of Wisconsin and lower Michigan have experienced increased freeze–thaw cycles in the late twentieth century relative to the early- and midcentury periods. Therefore, these regions may experience increased potential for soil erosion due to enhanced freezing and thawing over winter months.

  • Historically, soil temperatures in warm years were greater by up to 3°C than those in cold years with shallower frost penetration by 15 cm on average. Cold years were typically influenced by deeper snowpacks (higher SWE totals), which led to shallower penetration of soil frost, specifically in Minnesota and northern Wisconsin. This indicates that snow accumulation played a key role in soil frost formation and seasonal dynamics. During the cold season, the extreme wet years experienced especially high SWE, resulting in further reduced penetration of soil frost and warmer soil temperatures compared to dry years.

Long-term model simulations provided useful information on spatial and temporal patterns of cold-season variables for the entire six-state study domain. However, there are some limitations of this study. This study did not account for the feedback between vegetation and climate, as vegetation parameters were based on 1992/93 land cover data. This implies that the model results are not as adaptive to real-world changes in agricultural and natural vegetation. The model was also not configured to represent the effects of urbanization or the lakes and wetlands that are common in the northern parts of the study domain. Representation of the earlier-mentioned processes may improve estimations of various water and energy balance components in the model, reducing uncertainties in estimation of freeze–thaw cycles and cold-season processes.

Acknowledgments

We are thankful to NASA for providing funding for this research through the Grant NNG04GP13P. This manuscript is Purdue Climate Change Research Center (PCCRC) paper 0901.

REFERENCES

REFERENCES
Abdulla
,
F. A.
,
D. P.
Lettenmaier
, and
X.
Liang
,
1999
:
Estimation of the ARNO model baseflow parameters using daily streamflow data.
J. Hydrol.
,
222
,
37
54
.
Alley
,
R. B.
, and
Coauthors
,
2007
:
Summary for policymakers.
Climate Change 2007: The Physical Science Basis, S. Solomon et al., Eds., Cambridge University Press, 18 pp
.
Bowling
,
L. C.
, and
Coauthors
,
2003a
:
Simulation of high-latitude hydrological processes in the Torne-Kalix basin: PILPS Phase 2(e) 1: Experiment description and summary intercomparisons.
Global Planet. Change
,
38
,
1
30
.
Bowling
,
L. C.
, and
Coauthors
,
2003b
:
Simulation of high-latitude hydrological processes in the Torne-Kalix basin: PILPS Phase 2(e) 3: Equivalent model representation and sensitivity experiments.
Global Planet. Change
,
38
,
55
71
.
Bras
,
R. A.
,
1990
:
An Introduction to Hydrologic Science.
Addison-Wesley, 643 pp
.
Cherkauer
,
K. A.
, and
D. P.
Lettenmaier
,
1999
:
Hydrologic effects of frozen soils in the upper Mississippi River basin.
J. Geophys. Res.
,
104
,
19599
19610
.
Cherkauer
,
K. A.
, and
D. P.
Lettenmaier
,
2003
:
Simulation of spatial variability in snow and frozen soil.
J. Geophys. Res.
,
108
,
8858
.
doi:10.1029/2003JD003575
.
Cherkauer
,
K. A.
,
L. C.
Bowling
, and
D. P.
Lettenmaier
,
2003
:
Variable infiltration capacity cold land process model updates.
Global Planet. Change
,
38
,
151
159
.
Cooter
,
E. J.
, and
S. K.
Leduc
,
1995
:
Recent frost date trends in the north-eastern USA.
Int. J. Climatol.
,
15
,
65
75
.
Dye
,
D. G.
, and
C. J.
Tucker
,
2003
:
Seasonality and trends of snow-cover, vegetation index, and temperature in northern Eurasia.
Geophys. Res. Lett.
,
30
,
1405
.
doi:10.1029/2002GL016384
.
Easterling
,
D. R.
,
2002
:
Recent changes in frost days and the frost-free season in the United States.
Bull. Amer. Meteor. Soc.
,
83
,
1327
1332
.
Easterling
,
D. R.
, and
T. R.
Karl
,
2001
:
Potential consequences of climate variability and change for the midwestern United States.
Climate change impacts on the United States: The potential consequences of climate variability and change, National Assessment Synthesis Team Foundation Rep., 167–188
.
Fassnacht
,
S. R.
, and
E. D.
Soulis
,
2002
:
Implications during transitional periods of improvements to the snow processes in the Land Surface Scheme—Hydrological model WATCLASS.
Atmos.–Ocean
,
40
,
389
403
.
Frauenfeld
,
O. W.
,
T.
Zhang
, and
J. L.
McCreight
,
2007
:
Northern hemisphere freezing/thawing index variations over the twentieth century.
Int. J. Climatol.
,
27
,
47
63
.
Gaile
,
G. L.
, and
C. J.
Willmott
,
1984
:
On the evaluation of model performance in physical geography.
Spatial Statistics and Models, Kluwer, 443–460
.
Hamed
,
K. H.
, and
A. R.
Rao
,
1998
:
A modified Mann-Kendall trend test for autocorrelated data.
J. Hydrol.
,
204
,
182
196
.
Hamlet
,
A. F.
, and
D. P.
Lettenmaier
,
2005
:
Production of temporally consistent gridded precipitation and temperature fields for the continental United States.
J. Hydrometeor.
,
6
,
330
336
.
Hansen
,
M. C.
,
R. S.
Defries
,
J. R. G.
Townsheed
, and
R.
Sohlberg
,
2000
:
Global land cover classification at 1 km spatial resolution using a classification tree approach.
Int. J. Remote Sens.
,
21
,
1331
1364
.
Hardy
,
J. P.
, and
Coauthors
,
2001
:
Snow depth manipulation and its influence on soil frost and water dynamics in a northern hardwood forest.
Biogeochemistry
,
56
,
151
174
.
Haugen
,
R.
, and
G.
King
,
1998
:
Seasonal frost depths, midwestern USA. Circumpolar Active-Layer Permafrost System (CAPS), version 1.0.
National Snow and Ice Data Center, Boulder, CO, digital media. [Available online at http://nsidc.org/data/ggd498.html]
.
Hirsch
,
R. M.
,
D. R.
Helsel
,
T. A.
Cohn
, and
E. J.
Gilroy
,
1992
:
Statistical treatment of hydrologic data.
Handbook of Hydrology, D. R. Maidment, Ed., McGraw-Hill, 17.1–17.55
.
Hodgkins
,
G. A.
, and
R. W.
Dudley
,
2006
:
Changes in late-winter snowpack depth, water equivalent, and density in Maine, 1926-2004.
Hydrol. Processes
,
20
,
741
751
.
Jordan
,
R.
,
1991
:
A one-dimensional temperature model for a snow cover: Technical documentation for SNTHERM.89.
U.S. Army Corps of Engineers, CRREL Special Rep. 91-16, 61 pp
.
Kalnay
,
E.
, and
Coauthors
,
1996
:
The NCEP/NCAR 40-Year Reanalysis Project.
Bull. Amer. Meteor. Soc.
,
77
,
437
471
.
Kling
,
G. W.
, and
Coauthors
,
2003
:
Confronting Climate Change in the Great Lakes Regions: Impacts on our Communities and Ecosystems.
The Union of Concerned Scientists and the Ecological Society of America Rep., 92 pp
.
Krause
,
P.
,
D. P.
Boyle
, and
F.
Base
,
2005
:
Comparison of different efficiency criteria for hydrologic model assessment.
Adv. Geosci.
,
5
,
89
97
.
Kunkel
,
K. E.
, and
Coauthors
,
1998
:
An expanded digital daily database for climatic resources applications in the midwestern United States.
Bull. Amer. Meteor. Soc.
,
79
,
1357
1366
.
Kunkel
,
K. E.
,
D. R.
Easterling
,
K.
Hubbard
, and
K.
Redmond
,
2004
:
Temporal variations in frost-free season in the United States: 1895–2000.
Geophys. Res. Lett.
,
31
,
L03201
.
doi:10.1029/2003GL018624
.
Legates
,
D. R.
, and
G. J.
McCabe
Jr.
,
1999
:
Evaluating the use of “goodness-of-fit” measures in hydrologic and hydroclimatic model validation.
Water Resour. Res.
,
35
,
233
241
.
Lemke
,
P.
, and
Coauthors
,
2007
:
Observations: Changes in snow, ice and frozen ground.
Climate Change 2007: The Physical Science Basis, S. Solomon et al. Eds., Cambridge University Press, 337–384
.
Liang
,
X.
,
D. P.
Lettenmaier
,
E. F.
Wood
, and
S. J.
Burges
,
1994
:
A simple hydrologically based model of land surface water and energy fluxes for general circulation models.
J. Geophys. Res.
,
99
,
14415
14428
.
Liang
,
X.
,
D. P.
Lettenmaier
, and
E. F.
Wood
,
1996
:
One-dimensional statistical dynamic representation of subgrid spatial variability of precipitation in the two-layer variable infiltration capacity model.
J. Geophys. Res.
,
101
,
21403
21422
.
Lohmann
,
D.
,
E.
Raschke
,
B.
Nijssen
, and
D. P.
Lettenmaier
,
1998a
:
Regional scale hydrology: I. Formulation of the VIC-2L model coupled to a routing model.
Hydrol. Sci. J.
,
43
,
131
142
.
Lohmann
,
D.
,
E.
Raschke
,
B.
Nijssen
, and
D. P.
Lettenmaier
,
1998b
:
Regional scale hydrology: II. Application of the VIC-2L model to the Weser River, Germany.
Hydrol. Sci. J.
,
43
,
143
158
.
Mao
,
D.
, and
K. A.
Cherkauer
,
2009
:
Impacts of land-use change on hydrologic responses in the Great Lakes region.
J. Hydrol.
,
374
,
71
82
.
Maurer
,
E. P.
,
A. W.
Wood
,
J. C.
Adam
,
D. P.
Lettenmaier
, and
B.
Nijssen
,
2002
:
A long-term hydrologically based dataset of land surface fluxes and states for the conterminous Unites States.
J. Climate
,
15
,
3237
3251
.
Miller
,
D. A.
, and
R. A.
White
,
1998
:
A conterminous United States multilayer soil characteristics dataset for regional climate and hydrology modeling.
Earth Interactions
,
2
.
[Available online at http://EarthInteractions.org]
.
Mitchell
,
K. E.
, and
Coauthors
,
2004
:
The multi-institution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrologic modeling system.
J. Geophys. Res.
,
109
,
1
32
.
Mote
,
P. W.
,
A. F.
Hamlet
,
M. P.
Clark
, and
D. P.
Lettenmaier
,
2005
:
Declining mountain snowpack in western North America.
Bull. Amer. Meteor. Soc.
,
86
,
39
49
.
Myneni
,
R. B.
,
R. R.
Nemani
, and
S. W.
Running
,
1997
:
Estimation of global Leaf Area Index and absorbed par using radiative transfer models.
IEEE Trans. Geosci. Remote Sens.
,
35
,
1380
1393
.
Nash
,
J. E.
, and
J. V.
Sutcliffe
,
1970
:
River flow forecasting through conceptual models part I—A discussion of principles.
J. Hydrol.
,
10
,
282
290
.
Nijssen
,
B.
,
G. M.
O’Donnell
,
D. P.
Lettenmaier
,
D.
Lohmann
, and
E. F.
Wood
,
2001a
:
Predicting the discharge of global rivers.
J. Climate
,
14
,
3307
3323
.
Nijssen
,
B.
,
R.
Schnur
, and
D. P.
Lettenmaier
,
2001b
:
Global retrospective estimation of soil moisture using the variable infiltration capacity land surface model, 1980–93.
J. Climate
,
14
,
1790
1808
.
Nikol’skii
,
A. A.
,
E. E.
Roshchina
, and
O. V.
Soroka
,
2002
:
Snow cover as a factor of winter ecology of small mammals in the steppes zone.
Dokl. Biol. Sci.
,
383
,
158
160
.
Rao
,
A. R.
,
K. H.
Hamed
, and
H-L.
Chen
,
2003
:
Time domain analysis.
Nonstationarities in Hydrologic and Environmental Time Series, Kluwer Academic Publishers, 27–54
.
Schaal
,
L. A.
,
J. E.
Newman
, and
K. A.
Scheeringa
,
1981
:
Climatology of soil temperatures in Indiana.
Department of Agronomy, Agricultural Experiment Station, Purdue University, Station Bulletin 307, 87 pp
.
Sinha
,
T.
, and
K. A.
Cherkauer
,
2008
:
Time series analysis of freeze and thaw processes in Indiana.
J. Hydrometeor.
,
9
,
936
950
.
Stewart
,
I. T.
,
D. R.
Cayan
, and
M. D.
Dettinger
,
2005
:
Changes toward earlier streamflow timing across western North America.
J. Climate
,
18
,
1136
1155
.
Storck
,
P.
, and
D. P.
Lettenmaier
,
1999
:
Predicting the effect of a forest canopy on ground snow pack accumulation and ablation in maritime climates.
Proc. 67th Western Snow Conf., South Lake Tahoe, CA, Western Snow Conference, 1–12
.
U.S. Department of Commerce
,
1966–1982
:
Indiana, Michigan, Illinois, Iowa, Minnesota, Wisconsin: Climatological Data, Vols. 71–87, Environmental Data and Information Service, National Climatic Data Center
.
WARM
,
cited
.
2008
:
Illinois State Water Survey.
.
Willmott
,
C. J.
,
1981
:
On the validation of models.
Phys. Geogr.
,
2
,
184
194
.
WMO
,
2006
:
Statement on the status of the global climate.
WMO Rep. 998, 12 pp. [Available online at http://www.wmo.ch/index-en.html]
.
Wood
,
E. F.
,
D. P.
Lettenmaier
,
X.
Liang
,
B.
Nijssen
, and
S. W.
Wetzel
,
1997
:
Hydrologic modeling of continental-scale basins.
Annu. Rev. Earth Planet. Sci.
,
25
,
279
300
.
Zhang
,
T.
,
2005
:
Influence of the seasonal snow cover on the ground thermal regime: An overview.
Rev. Geophys.
,
43
,
RG4002
.
doi:10.1029/2004RG000157
.
Zhang
,
T.
,
R. G.
Barry
,
D.
Gilichinsky
,
S. S.
Bykhovets
,
V. A.
Sorokovikov
, and
J.
Ye
,
2001
:
An amplified signal of climatic change in soil temperatures during the last century at Irkutsk, Russia.
Climatic Change
,
49
,
41
76
.
Zhang
,
T.
,
M.
Serreze
,
R. G.
Barry
,
D.
Gilichinsky
, and
A.
Etringer
,
2003
:
Climate change: Evidence from Russian historical soil temperature measurements.
Geophysical Research Abstracts, Vol. 5, Abstract 1485. [Available online at http://www.cosis.net/abstracts/EAE03/01485/EAE03-J-01485.pdf]
.

APPENDIX A

Location and Soil Types for Sites Collecting Soil Temperatures at 10-cm Depth

Footnotes

Corresponding author address: Tushar Sinha, School of Life Sciences, Arizona State University, P.O. Box 874501, Tempe, AZ 85287-4501. Email: tsinha4@asu.edu