Response of seasonal soil freeze depth to climate change across China

The response of seasonal soil freeze depth to climate change has repercussions for the surface energy and water balance, ecosystems, the carbon cycle, and soil nutrient exchange. Despite its importance, the response of soil freeze depth to climate change is largely unknown. This study employs the Stefan solution and observations from 845 meteorological stations to investigate the response of variations in soil freeze depth to climate change across China. Observations include daily air temperatures, daily soil temperatures at various depths, mean monthly gridded air temperatures, and the normalized difference vegetation index. Results show that soil freeze depth decreased significantly at a rate of −0.18± 0.03 cm yr−1, resulting in a net decrease of 8.05± 1.5 cm over 1967–2012 across China. On the regional scale, soil freeze depth decreases varied between 0.0 and 0.4 cm yr−1 in most parts of China during 1950–2009. By investigating potential climatic and environmental driving factors of soil freeze depth variability, we find that mean annual air temperature and ground surface temperature, air thawing index, ground surface thawing index, and vegetation growth are all negatively associated with soil freeze depth. Changes in snow depth are not correlated with soil freeze depth. Air and ground surface freezing indices are positively correlated with soil freeze depth. Comparing these potential driving factors of soil freeze depth, we find that freezing index and vegetation growth are more strongly correlated with soil freeze depth, while snow depth is not significant. We conclude that air temperature increases are responsible for the decrease in seasonal freeze depth. These results are important for understanding the soil freeze–thaw dynamics and the impacts of soil freeze depth on ecosystem and hydrological process.


Introduction
Combining multiple land and ocean surface temperature datasets, the global mean air temperature increased 0.85 • C over 1880-2012 (Stocker et al., 2014).Given that all of the cryosphere's components are inherently sensitive to air temperature changes on different timescales, cryospheric changes serve as indicators of climate change.Frozen ground is an important component of the cryosphere.Permafrost regions underlay approximately 24 % of the exposed land surface of the Northern Hemisphere (Zhang et al., 1999), and seasonally frozen ground (SFG) regions occupy 57 % (Zhang et al., 2003).China has the third-largest frozen ground extent in the world, with a permafrost area of ∼ 2.20 × 10 6 km 2 , or approximately 23 % of its land area, mainly on the Tibetan Plateau; regions with SFG occupy about 50 % of the land area in China (Zhou et al., 2000;Christiansen et al., 2010).Under warming climate conditions, frozen ground regions are vulnerable to subsidence, especially ice-rich permafrost and relatively warm discontinuous permafrost (Morison et al., 2000;Osterkamp et al., 2000;Stendel and Christensen, 2002).Maximum soil freeze depth (SFD) of SFG and active Published by Copernicus Publications on behalf of the European Geosciences Union.layer depth over permafrost play a significant role in cold environments, and all hydrological, ecological, biological, and pedological activities occur within this layer (Hinzman et al., 1991;Kane et al., 1991;Zhao et al., 2004).Simultaneously, SFD influences the surface and subsurface hydrologic cycle, promotes soil texture changes, and alters the availability of soil nutrients for plant growth.The soil freeze-thaw cycle and SFD variations affect the decomposition of soil organic matter and greenhouse gas exchanges between the land surface and the atmosphere (Shiklomanov and Nelson, 2002;Mu et al., 2015;Jafarov and Schaefer, 2016).Thus, seasonal SFD variability and climate are closely linked.
Due to global climate warming, significant efforts have been devoted to permafrost research, such as permafrost variations on the hemispheric-scale permafrost temperature changes (Wu and Zhang, 2008;Romanovsky et al., 2010;Guglielmin and Cannone, 2012;Streletskiy et al., 2014;Wu et al., 2015), permafrost degradation (Jorgenson et al., 2006;Ravanel et al., 2010;Sannel and Kuhry, 2011;Streletskiy et al., 2015a;Park et al., 2016), hydrological processes in permafrost regions (Hu et al., 2009;Wang et al., 2009;Park et al., 2013;Streletskiy et al., 2015b;Ford and Frauenfeld, 2016), feedbacks to climate change (Schuur et al., 2008;Park et al., 2015;Abbott et al., 2016), and other aspects.The increasing thickness of the active layer has been indicated by many observations in permafrost regions at high latitudes and altitudes (Brown et al., 2000;Frauenfeld et al., 2004;Zhang et al., 2005;Fyodorov-Davydov et al., 2008;Smith et al., 2010;Wu et al., 2010;Zhao et al., 2010;Callaghan et al., 2011;X. Li et al., 2012;Liu et al., 2014a, b;Stocker et al., 2014).Less research has focused on SFG areas (Zhang et al., 2003;Frauenfeld et al., 2004;Frauenfeld and Zhang, 2011;Wang et al., 2015), although the near-surface soil freezethaw status has been investigated using satellite passive microwave remote sensing (Zhang and Armstrong, 2001;Zhang et al., 2003Zhang et al., , 2004;;Li et al., 2008;Jin et al., 2015).Peng et al. (2016b) analyzed the response of soil freeze-thaw states to climate change across China, based on observational data.While Peng et al. (2016b) investigated the area extent changes of different soil freeze-thaw states, here we instead focus on seasonal SFD.Regional-scale SFD can be an important indicator of climate change and frozen ground condition in cold regions.Further, SFG is closely related with human activities, because most populated areas are located on SFG.Shiklomanov (2012) similarly pointed out that SFG has not received much attention despite its vast area extent and importance, mainly due to a lack of long-term observational time series to document changes.Evaluating climatic and environmental changes on SFG requires comprehensive spatial assessments of available soil temperature records (Shiklomanov, 2012).To date, no comprehensive investigation of SFD in relation to climate change has been conducted in China despite the prevalence of SFG in this part of the world.Therefore, using long-term observational data, the goals and unique contributions of this study are (1) to estimate the spa- 2 Data and methods

Data
Several datasets are used including daily air and ground surface temperature, daily soil temperature at 0-320 cm depth, mean monthly gridded air temperature, and daily snow depth.In addition, we incorporate a 1 km resolution digital elevation model (DEM) and normalized difference vegetation index (NDVI) data.All datasets are described in detail below.

Mean daily air and ground surface temperature
Mean daily air temperature and ground surface temperature data are collected from the China Meteorological Administration (CMA) for a total of 839 meteorological stations (Fig. 1) available four times daily at 02:00, 08:00, 14:00, and 20:00 (http://data.cma.cn/en;Wang et al., 2015).These data come already quality controlled, and station observations date back to the 1950s and 1960s.Some stations end during the 1990s, while others are available until 2013.Most stations are located in east central China, with fewer sites in the west and at high elevations, such as on the Qinghai-Tibetan Plateau (Fig. 1).These mean daily air and ground surface temperatures are used to estimate temperature changes and to calculate the freeze-thaw index.

Soil temperature
Daily soil temperature data are available for 845 sites across China (Fig. 1) from the CMA, measured at the depths of 0.00, 0.05, 0.1, 0.15, 0.2, 0.4, 0.5, 0.8, 1.6, and 3.2 m.The temporal record varies for these stations, with some observations dating back to the late 1950s, and some only to the 1970s.Some station records end in the 1990s, while others are available through 2006 (Wang et al., 2015).Soil temperature is used to calculate the soil freeze depth; we combine the potential maximum soil seasonal freeze depth in permafrost regions, and the maximum soil freeze depth in SFG.The number of stations with both daily air temperature and soil temperature observations is 729.

Mean monthly gridded air temperature (MMGAT)
Mean monthly gridded air temperature was used to analyze soil freeze depth at the regional scale across China.
We obtained the University of Delaware's 1900-2014 terrestrial air temperature gridded monthly time series (http: //climate.geog.udel.edu/~climate/),with a 0.5 • × 0.5 • spatial resolution.This dataset was produced by combining many observational station records across the world, using spatial interpolation and cross-validation procedures (Legates and Willmott, 1990;Willmott and Matsuura, 1995;Willmott and Robeson, 1995;Peterson and Vose, 1997;Peterson et al., 1998).MMGAT during 1950-2010 is used for assessing its correspondence with seasonal freeze depth across China.

DEM
Considering the complex terrain across China and the impacts of elevation on air temperature, we also used the global 30 arcsec elevation dataset (GTOPO30; https://lta.cr.usgs.gov/GTOPO30) as the DEM for this study to further improve the MMGAT resolution.GTOPO30 was derived from several raster and vector sources of topographic information.Across China, the elevation ranges from −152 to 8752 m (Fig. 1).Based on this DEM, we spatially interpolate the MMGAT data to the DEM's 30 arcsec (1 km) resolution.

Snow depth
We obtained daily mean snow depth data for 672 sites across China (Che et al., 2008).The period of record at these locations varies, with some stations dating back to the late 1950s and some only to the 1970s.Some station records end around the 1990s while others are available through 2005.The snow depth was used to assess its influence on soil freeze depth.We calculate the annual maximum snow depth (SND) from the daily data for 1 July-30 June and match those snow depth stations with the soil temperature stations.If there are missing data in the spring, autumn, and winter season of one station, these station data will not be used.

NDVI
The NDVI dataset used in this study is produced by the Global Inventory Modeling and Mapping Studies (GIMMS) team, available for 1982-2006.It is derived from NOAA AVHRR data, available at 15-day temporal resolution and an 8 km spatial resolution (Tourre et al., 2008).These data were used to assess the influence of vegetation on soil freeze depth.We extracted the NDVI values corresponding to the stations' latitude and longitude coordinates.

Methods
Missing data often present a potential problem for analyzing and averaging time series.Therefore, if fewer than 5 days were missing in a given month, filling in missing daily air temperatures was based on highly correlated neighboring sites using linear regression.Missing daily mean ground surface temperatures were estimated through linear regression with the daily mean air temperature at the same station.Based on the daily air temperature, we also calculate the mean monthly air temperature and mean annual air temperature (MAAT).The interpolated results are strongly correlated with observations, as indicated by regression coefficients larger than 0.95.
To improve the original 0.5 • × 0.5 • MMGAT data to a 1 km resolution, spatial interpolation was used in conjunction with monthly lapse rates and the 1 km resolution DEM (e.g., Willmott and Matsuura, 1995;Gruber, 2012).The data processing steps are to (1) calculate the average monthly atmospheric lapse rate based on all available meteorological stations across China and their elevations; (2) bring each average monthly gridded air temperature value to a reference level (elevation of 0 m) using the average monthly lapse rate; (3) apply a Kriging interpolation to the reference-level adjusted MMGAT; and (4) bring the gridded reference-level air temperature back to the DEM-gridded height.Based on more than 800 sites, we evaluated the interpolated MMGAT against the observational monthly air temperatures and found that the regression coefficient was almost 1.0 with a minimum of 0.98 in April.
The freezing and thawing index can also be an important indicator to assess the variations in frozen ground (Zhang et al., 1997;Nelson, 2003;Frauenfeld et al., 2007).There are two primary types of freezing and thawing indices: the surface freezing and thawing index, calculated from ground surface temperatures, and the air freezing and thawing index, computed from air temperatures.To calculate the freezing index, we sum all temperatures below 0 • C during the freezing periods (Eq. 1) and similarly calculate the thawing index by summing the above 0 • C temperatures during the warm season (Eq.2; Wu et al., 2011;Luo et al., 2014).We define the freezing period to be July-June, to sum the freezing index over a continuous cold season (Eq.1).The warming period is defined as the calendar year (Wu et al., 2011;Peng et al  2013, 2016a) (Eqs. 1, 2).Thus, freezing and thawing index at the point scale was calculated based on the daily mean air temperatures and ground surface temperatures.For the regional-scale air freezing index, we use the adjusted 1 km gridded terrestrial air temperature data.
where N F is the number of days with temperature below 0 • C, i = 1, 2. ..N F , N T is the number of days with temperatures above 0 • C, and T i represents the temperature on a specific day.
Various methods are available to calculate the soil freeze depth.For example, it can be estimated directly from soil temperature, from physical and statistical models, and based on the Stefan solution.In this study, we use the Stefan solution to estimate soil freeze depth, which is determined using Eq. ( 3): where SFD is soil freeze depth (m), K f is the thermal conductivity of the frozen soil (W m −1 • C), n f is the n factor for the freezing season and corresponds to the ratio between the surface freezing index and the air freezing index (Peng et al., 2016), FI a is the annual air freezing index ( • C d −1 ), P b is the soil bulk density (kg m −3 ), w is the soil water content by weight, and L is the latent heat of fusion (J kg −1 ) (Zhang et al., 2005).In Eq. ( 3), many site-specific factors are required to estimate SFD, which are generally not available, particularly at the regional scale.However, based on the SFD and annual freezing index at each observational site, we can quantify the relationship between these two parameters (Fig. 2).We find a strong and statistically significant correlation of R = 0.87.Thus, the relationship between SFD and the annual freezing index can be simplified (Harlan and Nixon, 1978) as where E is defined (Nelson and Outcalt, 1987) as To estimate the SFD at the regional scale across China, we first calculate SFD for every observational station by interpolating the depth of the 0 • C isotherm throughout the 0.0-3.2m soil profile using the daily mean soil temperature (Frauenfeld et al., 2004).Next, we estimate the FI a based on the calculations in Frauenfeld et al. (2007).To estimate the E value for all stations, we use the SFD, FI a , and Eqs. ( 2) and (3).Then, we interpolate the E value to the regional scale at 1 km resolution using kriging in ArcGIS.The SFD is estimated across China based on Eq. ( 4), the 1 km E value, and FI a .We can then estimate the regional-scale SFD for each year from 1950 to 2009 across China and obtain the mean decadal SFD.Finally, we estimate the SFD trend at the regional scale across China based on regression analysis.
From the 1 km scale E factor values, we can extract every site's E factor based on the sites' latitude and longitude.Then, the air freezing index from the sites is used to calculate the annual soil freeze depth at every site by Eq. (4).To evaluate the result, we compare the observational SFD calculated from the soil temperatures and the simulated SFD derived from the Stefan method (Eq.4) in Fig. 3.The result demonstrates that the mean absolute error and root-mean-square error are 0.08 and 0.14 m, respectively, indicating that there is a good agreement between simulated and observational SFD.
A number of climatic and environmental variables including MAAT, mean annual ground surface temperature (MAGST), freezing index, thawing index, SND, and NDVI are selected to investigate the potential drivers of the observed long-term SFD changes across China.We use Pearson correlations to analyze the association between these variables and SFD and employ a 95 % significance level to assess the statistical significance for all analyses.

Soil freeze depth
Figure 4 shows the spatial variability and trends of SFD at every location.The highest SFD was mainly located in northeastern and northwestern China and the Tibetan Plateau.In contrast, the lowest SFD was found in the south of China.Locations with SFD greater than 0.4 m are found north of the Yellow River.In the northwest of China, locations with SFD less than 0.8 m are found in the Taklimakan Desert, and some sites with SFD greater than 2.0 m are located in the Altai, Tianshan, and Pamir mountains.On the Tibetan Plateau, most sites have a SFD greater than 2.4 m.There is an increase in SFD with increasing latitude and elevation.The significant SFD changes are between −0.4 and less than 0 cm yr −1 .The sites with the strongest decreasing trends of −1.2 cm yr −1 are on Tibetan Plateau and −1.0 cm yr −1 in the north of China.
Figure 5 shows the standard deviation of SFD at each site across China.It varies from 0.00 to 0.27 m.The standard deviation of SFD is generally less than 0.03 m south of 35 • N, except on the Tibetan Plateau.In northeastern China, the standard deviation changes between 0.06 and 0.15 m.In the northwest, it is generally 0.06-0.12m.On the Tibetan Plateau, the standard deviation varies from less than 0.09 m but can be greater than 0.18 m at some sites.
Based on the sites' E factors and FI a , we calculate SFD time series anomalies from 1951 to 2012 (Fig. 6).Although a composite time series of all available stations data can be calculated during 1951-2012, few of 839 stations actually contribute to the mean values before the 1960s (Fig. 6).There are fewer than 200 stations in the early years, which therefore does not represent the SFD across China as a whole.Beginning in 1967, more than 800 stations contribute to each year's mean; therefore long-term SFD trends will only be evaluated from then on.There is a statistically significant change in SFD anomalies of −0.18 ± 0.03 cm yr −1 , corresponding to a net decrease of 8.05 ± 1.5 cm.In addition to the overall long-term decrease, there are also some patterns of interdecadal variability during 1967-2012, including slight positive changes in some periods.SFD exhibited both increases and decreases until 1975, followed by a sharp decrease until 1990.However, SFD has remained constant or may perhaps be increasing slightly during 1990-2012.Therefore, the overall SFD change during 1967-2012 was largely controlled by the decrease during 1975-1990.Similar SFD changes, attributable to variability in the North Atlantic Oscillation, were found in high-latitude Eurasia (Frauenfeld and Zhang, 2011).

Spatial and temporal variability of SFD in China
Based on the 1 km resolution E factor and 1 km FI a calculated from MMGAT, we estimate SFD across China from 1950 to 2009 by the Stefan method.Figure 7 shows the spatial pattern of multi-year mean SFD.SFD increases with latitude and elevation, with SFD greater than 1.5 m in northeastern China, the Mongolia Plateau, Tibetan Plateau, and north of the Xinjiang region.In the east of China, the SFD ranges from 0.0 m to more than 4.0 m, and increases with latitude.In the Yellow River region, the elevation decreases from west to east, while the SFD varies from greater than 2.5 m to less than 0.5 m.The SFD in the Taklimakan Desert is lower than in the surrounding area.Figure 8 demonstrates the spatial variability of SFD anomaly for the decades of the 1950s, 1960s, 1970s, 1980s, 1990s, and 2000s, with respect to the 1950-2009 mean across China.Results show that the spatial changes of SFD anomaly ranges from larger than 0.15 m to less than −0.1 m.From the 1950s to 2000s, the SFD anomaly changes from positive to negative, which means that SFD has a decrease trend during this period.In addition, the smaller variability of SFD anomaly occurs in south of Yellow River and desert regions in Xinjiang.The SFD trend ranges between 0.0 and −0.4 cm yr −1 in most areas.SFD trends less than −0.4 cm yr −1 are found in some areas, such as the Tibetan Plateau and the Pamirs.In the two small areas of increasing SFD, we further investigated  1950s, 1960s, 1970s, 1980s, 1990s, and 2000s, with respect to the 1950-2009 mean across China.
the MAAT trend during 1950-2010 based on the MMGAT dataset.There is similarly a statistically significant decrease of MAAT in these same areas during this period.Therefore, air temperature is possibly one of the important factors that influence SFD in these areas.More detailed discussion is provided in Sects.3.3 and 4.1.
Overall, the spatial variability indicates that SFD changes with latitude and elevation at the regional scale across China.As is expected from climate warming, a statistically significant decreasing trend in SFD is evident across China from 1950 to 2009.

Potential forcing variables
To explore the possible variables leading to the documented changes in SFD, we analyze potentially important factors for soil freeze dynamics: latitude, altitude, MAAT, MAGST, freezing index including FI a and the ground surface freezing index (FI s ), thawing index including the air (TI a ) and ground surface thawing index (TI s ), SND, and NDVI.
To explore the spatial variability of SFD, we classify the meteorological stations as either eastern or western based on 110 • E longitude.Figure 10 represents the correlations between SFD and latitude and altitude in the eastern and western parts.In the east, we find an exponential relationship between SFD and latitude and a linear relationship with altitude, with both being statistically significant.The SFD values range from 0.0 m to less than 3.5 m, varying with latitude more so than with altitude.Thus, SFD is mainly affected by latitude in eastern China.In the west, SFD is near 0.0 m where altitude is higher than 1000 m.Similarly, SFD is related statistically significantly with altitude and latitude, but altitude is the main factor affecting SFD in the west.
Temperature -including MAGST and MAAT -at the 839 station locations exhibits a statistically significant increase over the 1951-2013 period of 0.019 and 0.013 • C yr −1 , or approximately 1.2 • C and 0.78 • C over the 63 years, respectively (Fig. 11a and b).MAGST and MAAT are statistically significantly correlated with SFD at R = −0.56 and R = −0.66,which means that 31 and 44 %, respectively, of the variability in SFD can be accounted for by these temperature measures.Further, the negative correlation demonstrates that increasing temperatures result in SFD decreases at the 839 stations.
Soil freeze usually begins in autumn or winter, with temperatures less than 0 • C reaching their maximum freeze depth toward the end of winter season or spring.Therefore, maximum annual SFD occurs during the cold seasons.Freezing index is an important indicator for accumulated cold season temperatures (Frauenfeld and Zhang, 2011).From 1951 to 2013, FI s and FI a underwent a statistically significant decrease of 3.0 and 1.62 • C days yr −1 , respectively (Fig. 11c  and d), indicating warming, which reduces the cold season's magnitude and/or duration.The correlation between FI s , FI a , and SFD was a statistically significant 0.68 and 0.87, indicating that the FI accounts for 46 and 76 % of SFD variability.
The thawing index is used to assess the accumulated positive degree days during the warm season (Frauenfeld and Zhang, 2011)  to 2013, TI s and TI a show statistically significant increases at a magnitude of 3.73 and 2.77 • C days yr −1 , respectively (Fig. 11e and f).The correlation coefficient between TI s , TI a , and SFD is −0.53 and −0.57, respectively, indicating a weak negative association, such that warm summer conditions correspond to a shallower SFD the following cold season.
Figure 12 shows the correlation between SFD and SND.However, the weak negative correlation between SFD and SND of R = −0.13 is not statistically significant, indicating that there is no relationship.This result is also consistent with the findings of Jafarov and Schaefer (2016).
As suggested by Shiklomanov (2012), environmental factors likely also affect SFD.The surface can be affected directly by climate forcing, while the subsurface effects are more complex.The subsurface soil only indirectly receives a climatic signal, which is furthermore altered by site-specific soil processes (e.g., thermal conductivity and analogous soil properties).Vegetation is a likely environmental factor that influences the soil freeze depth (Shiklomanov, 2012).Thus, we investigate vegetation using NDVI (Peng et al., 2013) and find it is significantly correlated with SFD at −0.80, suggesting that 64 % of the variability in SFD can be accounted for by NDVI.The statistically significant negative correlation demonstrates that when NDVI increases (greening) this corresponds to a decrease in SFD (Fig. 13).

Discussion
Soil freeze-thaw depth changes involve a series of interactions, such as energy exchanges, soil moisture exchanges, and gas exchanges between the atmospheric and terrestrial system.Therefore, variations of soil freeze-thaw most likely have an important effect on geomorphic, hydrological, and biological processes.Similarly, soil freeze-thaw depth changes also have destabilizing effects on engineering structures, such as on improperly constructed infrastructure (Smith and Burgess, 1999;Stendel and Christensen, 2002).The release of additional greenhouse gases to the atmosphere also occurs (Michaelson et al., 1996;Mu et al., 2015).In this paper, we use the Stefan method to calculate SFD, analyze the spatial SFD variability and trends, and quantify the potential driving factors affecting SFD.

Climatic and environmental factors
SFD variability is susceptible to climate warming and environmental change and is affected by variables including air temperature, ground surface temperature, freezing and thawing index, and vegetation.Many examples of permafrost degradation have been reported, such as deeper the active layer thickness, reduced freeze time duration, and shifts in the timing of thawing and freezing in seasonally frozen ground regions (Henry, 2008;Callaghan et al., 2011;Stocker et al., 2014;Wang et al., 2015).Negative correlations are found here between SFD and temperature (including MAAT and MAGST) because of solar radiation heating the ground and energy transfer into the soil, ultimately increasing the soil temperature.Thus, increasing temperature is found to be the main factor influencing SFD variability in China, as in previous work focusing only on the Tibetan Plateau (Zhao et al., 2004).
The freezing and thawing indices represent the accumulated negative and positive degree days in the cold and warm seasons, respectively (Wu et al., 2011).The positive and negative correlation between SFD and FI and TI was statistically significant, consistent with previous results in other regions (Frauenfeld and Zhang, 2011).Due to the maximum soil freeze depth occurring in the cold season and SFD being affected by temperature, the positive correlation between SFD and FI is reasonable.Although TI is the accumulated www.the-cryosphere.net/11/1059/2017/The Cryosphere, 11, 1059-1073, 2017 Snow depth can have an effect on soil temperature, which would affect the active layer thickness and seasonal SFD variability.Numerical modeling studies have shown that snow depth does impact SFD (Zhang and Stamnes, 1998;Ling and Zhang, 2003;Park et al., 2015).Park et al. (2015) indicated that both increasing SND and snow structure (e.g., snow density) changes were favorable to soil warming, resulting in active layer thickness decreasing in northern regions as previously found by Frauenfeld et al. (2004).Snow cover insulates the ground during the cold season (Zhang, 2005).Interestingly, in our study we did not find a relationship between SND and SFD.This could be due to the spatial heterogeneity of snow across China.According previous research, snow depth, snow water equivalent, and snow densities are smallest on the Tibetan Plateau compared to other parts of China (Ma and Qin, 2012).Compared with other regions, multi-year average snow depth in general is low in China, especially on the Tibetan Plateau and the east-central mountain regions of China (Zhong et al., 2014), and may therefore have only limited insulating effects.This could lead to the lack of a relationship between SFD and SND across China and motivates future investigation.
A negative correlation between SFD and vegetation, as quantified by NDVI, is found.Vegetation change has a significant influence on the climate system mostly through changes to the surface radiative energy budget, which can be affected the SFD.Based on previous research, vegetation varies in different land cover types and responds to climate change via different physical mechanisms (Snyder et al., 2004), e.g., changes in the surface albedo (e.g., bare ground versus vegetation cover), vegetation transpiration, and shading effects (Kelley et al., 2004;Snyder et al., 2004;Swann et al., 2010;Chang et al., 2012;Zhang et al., 2012).In the cold season, less/decreased vegetation will be more easily snow covered, thus increasing the albedo considerably.Increasing albedo results in less net radiation at the land surface, as more incoming solar radiation is reflected from the surface.Then, the surface air temperature will decrease considerably due to less energy absorbed at the surface.For the colder land surface, the sensible heat flux is reduced.Further, the vegetation decrease results in reducing evapotranspiration, which decreases the latent heat flux (Snyder et al., 2004).Compared to increased vegetation cover, less vegetation causes a large annual-average increase in the surface albedo with the largest changes in the winter and spring seasons, which reduces the amount of net radiation at the surface, making the surface colder and resulting in SFD increases.Conversely, vegetation increases could lead to decreasing SFD.The vegetation's effect on transpiration is primarily important in summer, while SFD primary occurs in winter and spring (Snyder et al., 2004).
The inverse relationship between NDVI and SFD is in agreement with results from many previous studies that similarly found a vegetation increase, or a greening trend, in different regions during recent decades (Peng et al., 2011;Piao et al., 2011;Zhang et al., 2013;Zhu et al., 2016).Because climate change controls the spatial distribution of vegetation, most studies report vegetation increases as impacted by temperature and precipitation increases (Bao et al., 2015;Huang et al., 2016).Similarly, Fig. 11 shows that rising temperature results in a SFD decrease.The negative relationship between SFD and NDVI indicates the effect of vegetation on SFD as well as their inverse relationship.
SFD is affected by many factors, including the climatic and environmental variables considered in this study.However, SFD changes in different regions are also potentially influenced by many other local environmental variables or large-scale teleconnections (Frauenfeld and Zhang, 2011).Thus, it remains difficult to fully account for the spatial variations of SFD at the regional scale.

Soil freeze depth in different climate zones
Our results indicate significant changes of SFD across China.To address the spatial pattern of SFD changes, we divide the study area into five different zones, including tropical monsoon (TPM), subtropical monsoon (SM), temperate monsoon (TM), temperate continental (TC), and Qinghai-Tibetan alpine (QTA) climate zones, which are categorized by temperature, precipitation, and other parameters.Results indicate that the 30-year (1971-2000) average SFD in the SM, TM, TC, and QTA climate zones are 2.8 ± 0.5 cm, 113.6 ± 7.6 cm, 132.5 ± 7.8 cm, and 165.8 ± 6.7 cm, respectively.Similar changes of SFD are found across the TM (−0.27 ± 0.005 cm yr −1 ), TC (−0.26 ± 0.005 cm yr −1 ), and QTA (−0.22 ± 0.004 cm yr −1 ) zones during 1950-2009, while there are no significant changes in the SM region.This is likely due to the higher temperatures in SM climate zone (Fig. 14).Although this study investigates a number of environmental and climatic driving variables of SFD, the degree to which other potential factors (e.g., soil texture, soil moisture, albedo, cloud cover, teleconnections) could also influence SFD remains unknown due to a lack of reliable data.

Summary and conclusions
In this study, we conducted a comprehensive regional-scale investigation of SFD across China.A significant climate indicator, SFD is influenced by many variables including climatic and environmental factors.These factors are often integrated to affect SFD (Lachenbruch and Marshall, 1986;Brown et al., 2000;Frauenfeld et al., 2004).Our results can be summarized as follows.
The spatial distribution of SFD variability is influenced by latitude and elevation across China.High-latitude and highaltitude sites are characterized by large SFD.In contrast, smaller SFD values are generally observed for lower latitude and lower elevation regions.
Of the total 839 sites, we find that the SFD decreased significantly, at −0.18 ± 0.03 cm yr −1 from 1967 to 2012, equal to a net change of 8.05 ± 1.5 cm.The long-term decrease also exhibits interdecadal variability, including some positive changes in some periods and no change since 1990.On the regional scale, the 1950-2009 spatial variation of SFD ranges between 0.0 and 4.5 m across China, with most areas exhibiting significant decreases between less than 0.0 and −0.4 cm yr −1 .Different climatic and environmental factors were explored as potential driving variables of SFD.A negative relationship is evident between SFD and MAAT, MAGST, TI a , and TI s , with statistically significant correlations of −0.66, −0.56, −0.57, and −0.56, respectively.The climatic factors FI s and FI a were correlated positively with SFD, at 0.87 and 0.68, respectively.There is no correlation between SFD and SND.The environmental factor vegetation (NDVI) is negatively correlated with SFD, indicating that 64 % of the changes in SFD can be accounted for by vegetation.Of the potential drivers of SFD explored here, FI and NDVI are most strongly correlated with SFD in China, while SND did not show a significant association.
Data availability.The dataset of mean daily air, ground surface temperature, soil temperature, and snow depth from the China Meteorological Administration is not available for public use.It can be accessed only by scientific researchers in China through the submission of an application.The website is http://data.cma.cn/en.The mean monthly gridded air temperature dataset is open for public use and can be attained at http://climate.geog.udel.edu/~climate.The digital elevation model can be downloaded from https://lta.cr.usgs.gov/GTOPO30.The normalized differential vegetation index is open for public use and can be attained at http://westdc.westgis.ac.cn/data/1cad1a63-ca8d-431a-b2b2-45d9916d860d.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.The observational station distribution across China, including the 839 stations with air and ground surface temperatures (green symbols), 845 soil temperature stations (red symbols), and elevation.The blue solid lines represent the main rivers.

Figure 2 .
Figure 2. Linear least-squares regression between soil freeze depth and annual freezing index based on observational sites.The black solid line represents the linear regression.

Figure 3 .
Figure 3.Comparison of the simulated and observed SFD for all stations.The black solid line is the 1 : 1 line, while the gray dashed line is regression fit between the simulated and observed values.

Figure 4 .
Figure 4. Spatial distribution and variability of SFD at the observing stations.(a) Multi-year mean SFD at each site; (b) the number of sites contributing to the SFD mean; (c) the magnitude of SFD change at each site; (d) the number of sites with SFD change observations.

Figure 5 .
Figure 5.The standard deviation of SFD at each site across China.

Figure 6 .
Figure 6.1951-2012 SFD anomalies with respect to the 1971-2000 mean (red solid line) based on up to 839 stations across China as depicted in Fig. 1.Included also is the 1 standard deviation range (gray shading), the linear trend from 1967 to 2012 (blue dashed line), and the 7-year smoothing (green line).The inset shows the number of stations contributing to the time series.

Figure 9
Figure9represents the SFD trend across China from 1950 to 2009.The gray region represents areas where the SFD trends are not statistically significant, but they are statistically significant in all other regions.In general, the SFD decreased significantly over northern China, except in two small areas.The SFD trend ranges between 0.0 and −0.4 cm yr −1 in most areas.SFD trends less than −0.4 cm yr −1 are found in some areas, such as the Tibetan Plateau and the Pamirs.In the two small areas of increasing SFD, we further investigated

Figure 9 .
Figure 9. SFD trends across China from 1950 to 2009.The gray regions indicate non-significant SFD changes, while trends in all other regions are statistically significant.

Figure 10 .
Figure 10.The relationship between SFD, latitude, and elevation in the east and west of China.

Figure 11 .Figure 12 .
Figure 11.SFD time series and trend (black) and the potential forcing variables: (a) mean annual ground surface temperature (red), (b) mean annual air temperature (green), (c) surface freezing index (cyan), (d) air freezing index (magenta), (e) surface thawing index (yellow), and (f) air thawing index (orange).All variables are standardized to range from 0 to 1. R is the correlation coefficient, and all are statistically significant.

Figure 13 .
Figure 13.Correlation between SFD and mean annual NDVI.

Figure 14 .
Figure 14.Time series of SFD changes in different climate zones: (a) subtropical monsoon, (b) temperate monsoon, (c) temperate continent, (d) Qinghai-Tibetan Plateau alpine, and (e) tropical.The insets are the SFD changes in the four climate zones; the bold black line is SFD and the bold red line is the trend.