An observational 71-year history of seasonally frozen ground changes in the Eurasian high latitudes

In recent decades, significant changes have occurred in high-latitude areas, particularly to the cryosphere. Sea ice extent and thickness have declined. In land areas, glaciers and ice sheets are experiencing negative mass balance changes, and there is substantial regional snow cover variability. Subsurface changes are also occurring in northern soils. This study focuses on these changes in the soil thermal regime, specifically the seasonally frozen ground region of Eurasia. We use a database of soil temperatures at 423 stations and estimate the maximum annual soil freezing depth at the 387 sites located on seasonally frozen ground. Evaluating seasonal freeze depth at these sites for 1930–2000 reveals a statistically significant trend of −4.5 cm/decade and a net change of −31.9 cm. Interdecadal variability is also evident such that there was no trend until the late 1960s, after which seasonal freeze depths decreased significantly until the early 1990s. From that point forward, likely through at least 2008, no change is evident. These changes in the soil thermal regime are most closely linked with the freezing index, but also mean annual air temperatures and snow depth. Antecedent conditions from the previous warm season do not appear to play a large role in affecting the subsequent cold season’s seasonal freeze depths. The strong decrease in seasonal freeze depths during the 1970s to 1990s was likely the result of strong atmospheric forcing from the North Atlantic Oscillation during that time period.


Introduction
Soil seasonal freezing and thawing processes play a significant role in middle-and high-latitude climate. The ground thermal regime is influenced by a variety of factors, including air and ground surface temperature, terrain (slope, aspect, etc), surface and subsurface hydrology, vegetation, snow cover, geology, and natural and anthropogenic disturbances (Osterkamp 2007). The active layer depth gauges the maximum annual depth of thaw over permafrost and thereby represents an indicator of climate change (Nelson et al 1997, Brown et al 2000, Osterkamp 2007): long-term increases in active layer depth indicate that permafrost is degrading from the top down. Similarly, the maximum annual seasonal freeze depth (SFD) can be used to quantify long-term changes in seasonally frozen ground regions . In a warming climate, less and less land area experiences seasonal freezing, and the depth of the frost penetration also decreases over time. Because the maximum annual depth of thaw occurs at the end of the warm season (late summer to early autumn), active layer thickness over permafrost reflects primarily warm season conditions. However, climate has changed disproportionately during the cold season, i.e. the greatest observed warming has been observed during winter (IPCC 2007). Therefore, SFD, which reaches its maximum at the end of the cold season, is a unique indicator of climate change that integrates not only changes at the surface and above, but also changes at and below the ground surface. Hence, SFD depends not only on air temperature forcing, but also, potentially, on snow cover, warm season precipitation and resulting soil moisture variations (Iijima et al 2010), as well as a variety of site-specific soil properties: soil type, bulk density and the resulting porosity, and thermal and hydraulic conductivities. This paper therefore focuses on long-term changes in SFDs, as estimated from historical observations of soil temperature across the former Soviet Union, i.e. the Eurasian high latitudes. We do not include permafrost regions in this work and instead only focus on areas characterized by seasonal freezing and thawing processes. In previous work based on soil temperature observations at 242 locations, changes in SFD for the 1956-1990 period were evaluated on the basis of observations from 211 sites across Russia (Frauenfeld et al 2004). We expand upon that work with new data for 1990-2000, a period characterized by accelerated climate change that was unprecedented during the preceding observational record (Crowley 2000). In addition to rescuing data for the decade of the 1990s, we also expand our station database by 181 sites. Of these 423 total locations, 387 are located on seasonally frozen ground across the Eurasian high latitudes. Furthermore, additional data were also rescued for the earlier periods, which allows us to now characterize the 71-year period from 1930 to 2000 using 387 sites, as compared to the previous 35-year period, from 1956 to 1990, based on only 211 stations.
The goals and unique contributions of this study are thus to improve, expand and extend previous 1956-1990 analyses of the soil thermal regime in the Eurasian high latitudes (Frauenfeld et al 2004, Lemke et al 2007. Furthermore, we utilize a kriging approach to fill in many of the data gaps that plague raw soil temperature observations. Given the length and spatial coverage of this station database, this analysis represents the most comprehensive observational assessment of soil thermal conditions in the Eurasian high latitudes. We establish the importance of various potential driving variables, including measures of air temperature and precipitation. Furthermore, many of these surface variables are driven and potentially integrated by the synoptic-scale circulation of the atmosphere, and vice versa (Bojariu et al 2008). We therefore establish the relationship between SFD and the major mode of Eurasian atmospheric variability: the North Atlantic Oscillation (NAO).

Soil temperatures
Soil temperature data were obtained from the National Snow and Ice Data Center (NSIDC) and are available for 423 sites across the former Soviet Union (figure 1). The observations represent soil temperatures measured with extraction thermometers under natural surface cover such as grass and undisturbed snow cover during the cold season (Gilichinsky et al 1998) at depths from 0.2 to 3.2 m (0.2, 0.4, 0.6, 0.8, 1.2, 1.6, 2.0, 2.4 and 3.2 m). The period of record at these locations varies, with some stations dating back to the late 1800s , and some only to the 1990s (Frauenfeld et al 2004). Some station records end around 1950; while others are available through 2000 (Chudinova et al 2006, Frauenfeld et al 2004, Gilichinsky et al 1998.
Obtaining a robust average of soil temperatures over time is difficult when the number of stations varies. We therefore make best use of available data by filling in missing data points at stations using kriging, a method based on the historic relationship between stations. In order to compute a valid covariance between any two stations for any given month, we required a minimum of five concurrent data values to be available. The primary assumption of this kriging approach is that the anomaly statistics between stations are considered stationary over the period of record. A maximum of the closest 50 stations was used in any one estimation so as to minimize computational time and avoid issues of spurious correlation or overfitting. The inversion of the covariance matrix was performed via singular value decomposition. We cross-validated the kriging results by estimating the temperature for all stations for all times without using available data from the target station during its particular estimation. Estimates were then compared to observations and indicated maximum root mean squared errors of only 1-2 • C.

The seasonal freeze depth
The 423 stations were classified as either permafrost or seasonally frozen ground stations on the basis of the soil temperature at a depth of 3.2 m. If a particular 3.2 m soil temperature time series remained negative for that station's entire record, the site was considered to be underlain by permafrost and was therefore not included for further analysis. The 3.2 m soil temperature time series were also examined visually for each station to verify this classification.
To obtain freeze depths for the locations classified as seasonally frozen ground stations, we interpolate the depth of the 0 • C isotherm throughout the 0.2-3.2 m temperature profile. It is important to note, however, that while the depth of the 0 • C isotherm can be used as an estimate of the freeze depth it does not necessarily correspond to the 'true' freeze depth. For wet or ice-rich soils, latent heat can play an important role which may lead to errors in estimating the freeze depth on the basis of the propagation of the 0 • C isotherm. However, the robustness of this approach was validated in Frauenfeld et al (2004) where a virtually perfect relationship was verified between this interpolation procedure and actual observed SFDs. We follow the identical methodology here and derive the maximum freeze depth for the cold season for each seasonally frozen ground location. We next create departures with respect to the 1961-1990 mean to allow for comparison with previous results (Frauenfeld et al 2004). A composite of all SFD time series is created to represent average conditions across the whole of the Eurasian high latitudes. Given the heterogeneity and site-specific nature of soil conditions, we do not area-weight the individual time series prior to creating the composite. It is therefore important to note that the results presented here apply only at the 387 SFD station locations.

Potential forcing variables
To determine what may be driving potential changes in SFDs we analyze air temperature, freezing index, thawing index, snow depth and summer precipitation. For air temperature we use the monthly 0.5 • × 0.5 • gridded time series data, available for 1901-2002 from the University of East Anglia's Climatic Research Unit (CRU, Mitchell and Jones 2005). We calculate mean annual air temperatures and extract the temperature time series from the grid cell center closest to each of our station locations. For the freezing and thawing index we employ the methodology of Frauenfeld et al (2007) based on monthly 0.5 • × 0.5 • gridded terrestrial air temperature data using the climatologically aided interpolation (CAI) method of Willmott and Robeson (1995) from the University of Delaware's Center for Climatic Research, available for 1900-2008. The freezing index represents the cumulative degree-days of below 0 • C temperature in the course of a cold season, defined from 1 July of the current year to 30 June of the following year. Similarly, the thawing index represents the cumulative degree-days of above 0 • C temperature in the course of the warm season, calculated for the 1 January-31 December period. As for the mean annual air temperature data, the freezing and thawing index time series from the grid cell center closest to each station are extracted. Snow depth data, available for 1881-1995, are used from the Historical Soviet Daily Snow Depth Version 2 (HSDSD) from NSIDC (Armstrong 2001). This product represents daily snow depths at 284 World Meteorological Organization stations in the former Soviet Union. We calculate mean annual snow depths from the daily data for July-June and match up those snow depth stations that are within a 0.5 • radius of one of the soil temperature stations. Finally, precipitation data are obtained for the Eurasian high latitudes from the Serreze et al (2005) gauge precipitation dataset totaling 1583 stations with data for 1950-2003 across our domain. Because we are interested in warm season precipitation only, we sum the monthly precipitation totals for June, July and August (JJA). As for the SFDs, for each of the potential forcing variables, departures from the 1961 to 1990 mean are created and averaged into a time series representing conditions at the 387 SFD stations across the Eurasian high latitudes. It must be noted that, for all analyses, only those SFD stations that can be matched with potential forcing data are used. The number of stations used for the various analyses therefore varies slightly.
The combined climate variability of the aforementioned air temperature measures and precipitation (including snow depth) is likely related to the synoptic-scale variability of the atmosphere. We therefore also compare SFDs to the NAO. The station-based winter (December-March) index of the NAO, based on the difference of normalized sea level pressure between Lisbon, Portugal and Stykkisholmur/Reykjavik, Iceland, is obtained from the Climate Analysis Section of the National Center for Atmospheric Research (NCAR, Hurrell 1995). This index spans the period from 1864 to the present. We choose a 95% significance level for all statistical analyses presented here. These statistical procedures entail simple regression, correlations, and stepwise multiple regression (in addition to the procedures outlined in section 2.1).

The seasonal freeze depth
Of the total 423 sites with soil temperature observations, 387 were classified as seasonally frozen ground stations on the basis of the 3.2 m depth soil temperature (figure 1, red symbols). We calculate SFD at these 387 locations and create a composite time series to represent the temporal SFD variability across the Eurasian high latitudes. Prior to 1930, only approximately 200 stations contribute to the composite time series (not shown) while from 1930 onwards, with the exception of 1948, it is possible to calculate SFD on the basis of 250 or more of the 387 total stations (figure 2, green bars). We therefore focus our subsequent discussion on the 1930-2000 period.
Before increasing our soil temperature database from 242 stations to 423 (211 seasonally frozen ground stations to 387), it was only possible to reliably construct a composite SFD time series dating back to 1956, from which point forward 100 or more stations contributed to the mean (figure 2, blue bars). Therefore, the additional data and kriging efforts now allow us to focus on a 71-year time frame with almost double the number of stations, compared to the earlier results using a 35-year period from only 1956to 1990. Over that 1956-1990 period, the net change in SFD was found to be −34 cm (Frauenfeld et al 2004;their figure 6(b)). This indicated that, on average across the Eurasian high latitudes, the seasonally frozen ground during the cold season froze to a depth 34 cm less in 1990 than it did in 1956. Less robust because of the far smaller dataset, the net change over the 1930-1990 period was −27 cm. Accelerated changes were evident from the early 1970s onward, further suggesting a potential climate change signal (Frauenfeld et al 2004).
Our updated and improved time series now allows for a more robust trend estimate. As shown in figure 2, up to 320 (of the total 387) stations now contribute to each year's composite SFD departure time series. Comparing the 1930-1990 trend based on these new and improved data to the previous estimate, we now find a trend of −3.2 cm/decade, or a net change of −19.5 cm for the 61-year period. The 1956-1990 trend is −8.3 cm/decade, or a net decrease of 28.9 cm. These new estimates represent statistically significant decreases of similar, albeit 14-27% smaller, magnitude, compared with previous estimates. The likely reason that the previous assessment illustrated greater trends is because the composite time series (especially during the pre-1956 period) was based on fewer stations. This early period therefore exhibited much greater interannual variability and was more susceptible to outliers and leverage points (figure 3, gray line).
In addition to having more data and hence more robust trends going back to 1930, we now also have data for the decade of the 1990s. This addition of 10 years results in an overall statistically significant change for the 71-year period from 1930 to 2000 of −4.5 cm/decade and a net change of −31.9 cm (figure 3). In addition to this overall long-term decrease the time series also indicates some interesting patterns of interdecadal variability, including periods of slight positive changes. Freeze depths appear to have increased until ∼1970, and this was followed by a sharp decrease until ∼1990. From 1990 to 2000, however, freeze depth indicates little change and may actually be increasing slightly. The overall 1930-2000 change is therefore largely driven by the decrease during the ∼1970-1990 period.

Drivers of the seasonal freeze depth
We establish the correspondence between SFD and air temperature, freezing index, and snow cover. Also examined were warm season precipitation and the thawing index. Air temperature itself at the 387 station locations exhibits a statistically significant increase over the 1930-2000 period of 0.15 • C/decade, or approximately 1.1 • C over the 71 years ( figure 4(a)). From the late 1980s onward, air temperature  anomalies have been exclusively positive with respect to the 1961-1990 mean. These mean annual air temperatures are significantly correlated with the SFD time series at R = −0.71, indicating that 51% of the variance in SFD can be accounted for by changes in air temperatures. The negative correlation illustrates that, as air temperatures increased at the 387 stations, SFDs decreased from 1930 to 2000.
Because maximum annual SFDs usually occur toward the end of the cold seasons, the freezing index should be a better indicator of cold season conditions than mean annual temperatures. Over the 71-year period, the freezing index has decreased at a rate of 39 • C days/decade ( figure 4(b)), indicating that as the climate has been warming, the magnitude and duration of cold conditions has decreased. The freezing index exhibits some interesting interdecadal variability. No trend was evident at the 387 stations until approximately the 1970s. From around 1980 to 1990, a decrease occurred, followed by another period with no change during the 1990s. The overall long-term trend in freezing index is therefore driven primarily by the decrease during the 1980s. The correlation between freezing index and SFD is a statistically significant and strong R = 0.81, indicating that 65% of SFD variability can be accounted for by the freezing index. This strikingly strong correspondence between these two completely independent datasets is readily apparent in figure 4(b).
The thawing index represents the positive degree-days during the warm season. Its variability at the 387 seasonally frozen ground stations indicates a period of no changes from 1930 to approximately the late 1960s, and a positive trend thereafter ( figure 4(c)). The overall trend, driven primarily by the post-1970 period, is a statistically significant +14.1 • C days/decade. This trend is less than half the magnitude of the cold season changes in freezing index, illustrating that the magnitude of warming in the Eurasian high latitudes is much greater during the cold season than during the warm season. The correlation coefficient between the thawing index and SFDs is R = −0.48, indicating a weak (albeit statistically significant) negative association: a greater thawing index in summer corresponds to a lower SFD the following winter. There is thus potential evidence of warm season preconditioning in the soil thermal regime.
Precipitation changes in the Eurasian high latitudes likely also affect the soil thermal regime and hence SFD. Average snow depth observations, only available at a subset of up to 129 of the 387 SFD stations, indicate a high degree of interannual variability, especially prior to approximately 1970. This is likely an artifact of data coverage. In the early part of the record, only 30-50 of the HSDSD stations contribute to the composite snow depth time series. There is no statistically significant trend in average snow depth for the 1930-1995 period. However, snow depth is significantly but weakly correlated with SFD at R = −0.37 over 1930-1995. From 1950to 1995, when data coverage is better, this correlation increases to R = −0.51. Snow cover insulates the ground during the cold season (Zhang 2005), and the negative correlation indicates and verifies that when snow depth is greater, SFD is shallower.
Much like the temperatures during the previous summer (i.e. thawing index) being able to precondition the soil thermal regime, warm season precipitation might similarly have an effect. Increased rainfall in summer could impact soil moisture and hence the thermal conductivity of soils (Iijima et al 2010). Because the availability of station-based precipitation data is good for our domain, we are able to match virtually all of our SFD stations with precipitation data back to 1950. JJA precipitation at the 387 SFD stations exhibits both a high degree of interannual variability, as well as a small but statistically significant positive long-term trend of 0.72 mm/decade (figure 4(e)). Summer precipitation is also only weakly correlated with SFD (R = −0.27). While this negative correlation is statistically significant, its magnitude is so small that it is likely not a physically meaningful association.
While the above correlations for individual components of the climate system are useful for understanding their direct association with SFDs, in the real world these variables all jointly impact the soil thermal regime. We therefore employ a stepwise multiple-regression model approach to test the combined effect of air temperature, freezing and thawing index, snow depth, and JJA precipitation on SFD. The resulting model's R 2 indicates that it accounts for 86% of SFD variability (R = 0.93; adjusted R 2 = 85%). The significant variables in the model are air temperature, freezing index and snow depth. Thawing index and precipitation are not statistically significant. On the basis of the standardized coefficients, we can therefore summarize the model on the basis of the equation SFD = −0.28 air temperature + 0.55 freezing index − 0.36 snow depth. (1) This model indicates similar relationships to the correlation analysis. The same inverse relationship between SFD and both air temperature and snow depth is evident, as is the direct relationship with freezing index. Freezing index is again associated most strongly with SFD, as indicated by the magnitude of the standardized coefficient. Interestingly, even though snow depth was only weakly associated with SFD, it factors into the equation more strongly than air temperature. Considering all of the potential forcing variables jointly, we can thus conclude that freezing index is still the dominant variable driving SFD. However, the insulating effect of snow cover becomes more important than mean annual air temperature. Antecedent conditions such as thawing index and warm season precipitation are not significant drivers of SFD. It should be noted that a multiple-regression model that uses all possible variables, and models using different combinations of the five independent variables, always produce very similar results in terms of the sign and magnitude of the coefficients. We are therefore confident that the stepwise model presented here is robust and not susceptible to overfitting or multicollinearity issues, despite the use of multiple measures of air temperature. This is also illustrated by the close correspondence between the model's R 2 and adjusted R 2 (86% versus 85%, respectively).

The North Atlantic Oscillation
Another way to test the combined effect of many climate variables on the soil thermal regime is to relate SFDs to the NAO. The NAO provides a synoptic-scale atmospheric forcing that subsumes the individual forcings of variables such as air temperature, snow cover, and precipitation because it initially drives these individual variables. Correlating the NAO index with the SFD time series over the 1930-2000 period produces a statistically significant R of −0.58. The NAO therefore accounts for 34% of the variability in SFD such that during its positive phase, which is associated with mild conditions across Europe and Eurasia (Hurrell 1995), SFD depths decrease as expected. Interestingly, though, when the NAO is trending strongly positively, as was the case from approximately 1970 to the mid-1990s, its correlation with SFD increases to −0.89. In fact, over this approximate 30-year period, the NAO and SFDs almost perfectly mirror each other (figure 5).
To ensure that this strong correlation is not an artifact of both SFDs and NAO having a strong trend over this period, we detrend both variables. The correlation between the detrended NAO and SFDs is still a strong −0.81, suggesting that the association is not driven by the trend components. It is therefore likely that the period of accelerated SFD decreases observed from approximately 1970 to the late 1990s was driven by the NAO. The NAO is considered to be primarily a mode of stochastic internal atmospheric variability (Hurrell et al 2003). However, it is also possible that, if the NAO is subject to external forcing, this same forcing was driving SFD changes during the ∼1970-1990s period.

Changes since 2000
Given that SFDs appeared to be deepening again in the 1990s, a period characterized by increased global and high-latitude temperatures, an important question is: how may SFDs have changed during the 2000s, which was the warmest decade on record going back to at least 1880 (Hansen et al 2010)? The 2000s were also characterized by many record breaking cryospheric anomalies, including minima in sea ice extent (Stroeve et al 2007) and record Greenland melt extent (Tedesco et al 2011). Soil temperature observations for the period from 2000 to the present are not yet available from the NSIDC database used here. However, given the close correspondence between the freezing index and SFDs In fact, no statistically significant trend has been evident at all since approximately 1990. It is important to note, however, that even though SFDs did not continue to decrease, they did not return to pre-1970 levels: they have remained at the decreased 1990 levels.

Summary and conclusions
This paper updates and improves on earlier work characterizing changes in the soil thermal regime across the Eurasian high latitudes (Frauenfeld et al 2004). In that previous assessment, robust estimates over only the 1956-1990 period based on a total of 242 soil temperature stations were possible. Our updates presented here entail (1) the additional data rescue of 181 stations to a total of 423 (figure 1); (2) updates through the year 2000 for all stations, where updates were available; and (3) kriging efforts to fill in data gaps. We are now able to characterize SFD changes over 71 years, for 1930-2000, on the basis of 387 seasonally frozen ground stations. This represents 83% more stations and more than double the length of the original time series, and therefore provides the most comprehensive assessment of observed soil thermal conditions for the Eurasian high latitudes.
We find that SFD has decreased −4.5 cm/decade for 1930-2000, resulting in a net change of −31.9 cm. This overall decrease is primarily driven by changes in the 1970s and 1980s. SFD is very strongly linked with the freezing index. Cold season soil conditions also respond to mean annual air temperatures. Weak but statistically significant correlations with snow depth exist. Antecedent conditions from the previous warm season, as captured by the thawing index and JJA precipitation, are also somewhat evident. When these potential forcing variables are considered jointly using a stepwise multiple-regression approach, the freezing index, snow cover and mean annual air temperature account for 86% of SFD variability. Even though the effect of snow cover was relatively weak when individually correlated with SFD, in the multiple-regression model it emerged as the second-most important driver of SFD, after the freezing index. Antecedent conditions from the previous warm season (thawing index and JJA precipitation) did not factor into the model. Previous work had indicated the beginning of a strong decrease in SFDs from 1970 onward, suggesting observational evidence of a strong climate change signal in the soil thermal regime of the Eurasian high latitudes (Frauenfeld et al 2004). This updated analysis that now includes data through 2000 demonstrates that this decreasing trend did not continue in the 1990s (figure 3), and likely also has not continued through at least 2008 (figure 6). Furthermore, the strong decrease evident from ∼1970 to 1990 was likely driven by the NAO. Once the NAO's strong positive mode ended, so did the decrease in SFDs (figure 5). While the soil thermal regime tends to be highly site-specific given the heterogeneity of soils, this provides observational evidence for the importance of synoptic-scale climate forcing on annual to interannual time scales in high-latitude land areas.