Dynamics of seasonally frozen ground in the Yarlung Zangbo River Basin on the Qinghai-Tibet Plateau: historical trend and future projection

Seasonally frozen ground (SFG) is a critical component of the Earth’s surface that affects energy exchange and the water cycle in cold regions. The estimation of SFG depth has generally required intensive parameterization which has limited estimates in data-scarce regions such as the Qinghai-Tibet Plateau (QTP). We propose a simple yet robust modeling framework employing ground surface temperatures as major model inputs to assess the spatiotemporal patterns of the SFG depth in the Yarlung Zangbo River Basin (YZRB) on the QTP. The model was calibrated using SFG depth measurements throughout the YZRB from 1980 to 2010. Results suggest that the SFG depth in the YZRB has decreased at a rate of 2.50 cm · a−1 from 1980 to 2010. Future projections indicate that the SFG depth in the YZRB will continue to decrease in response to future warming. The present SFG may no longer exist by 2180 under the RCP 8.5 scenario (if not considering the transition of permafrost to SFG). The proposed modeling framework provides an important basis for the evaluation of the hydrological cycles (e.g. surface water-groundwater interactions) in cold regions under changing climatic conditions.


Introduction
Seasonally frozen ground (SFG) is defined as the portion of soil that freezes in winter and thaws in summer in cold regions (Evans and Ge 2017). The presence of SFG affects the seasonal infiltration in the vadose zone and thus the groundwater recharge (Ge et al 2011). SFG is a thermal condition dependent on ground thermal regime. With soil freezing and thawing, SFG influences the thermal and hydraulic conditions of the soil layer and thus impacts the regional hydrological cycle and ecosystem function (Yang et al 2002, Hansson et al 2004, Guo and Wang 2016. Assessing the spatiotemporal pattern of the SFG depth is essential for water resources management and environmental conservation in cold regions.
The Qinghai-Tibet Plateau (QTP) is the source region of several important large rivers (e.g. the Yellow River, Yangtze River, Lancang River, and Yarlung Zangbo River) and is often referred to as the water tower of the southeast Asia (Jin et al 2009). The Yarlung Zangbo River Basin (YZRB) is one of the most densely populated regions on the QTP and seasonally frozen soil broadly occurs here (Zhong et al 2014). Many studies have reported that the dramatic climate warming in recent decades has substantially accelerated the decrease of the SFG depth in this region (Zhao et al 2004, Pang et al 2011. However, long-term observations of SFG depth are lacking due to the high cost of direct measurements and maintenance of monitoring systems in cold regions. Hence, the assessment of the SFG depth over the large spatial scales often rely on computational models. Numerical schemes (e.g. finite element and finite diffrenence methods) that solve heat transport equations are often incorporated in hydrological models to simulate the soil freeze-thaw process (Gouttevin et al 2012, Rawlins et al 2013. For instance, the SUTRA model (Voss and Provost 2010) couples heat transport with groundwater flow and has been used to assess the impacts of SFG on groundwater flow patterns (Evans and Ge 2017). These numerical models solving heat transport equations provide sophisticated simulations of soil freezethaw dynamics, but the high computation cost and the large data demand for model calibration hinder their application to large scale investigations. Analytical or semi-analytical solutions relate the SFG depth to ground surface/air temperature and soil thermal-hydraulic properties. Among these models, the Stefan's solution (Stefan 1891) and the Kudryavtsev method (Kudryavtsev et al 1977) have received the widest application. The Stefan's solution uses air temperature or ground surface temperature (GST) and soil parameters (e.g. soil thermal conductivity, soil bulk density, soil water content) as inputs to estimate frozen ground depth (Jumikis 1977). The Kudryavtsev model is a semiempirical model originally developed for engineering purposes in the cold regions in Russia and considers the effects of soil moisture, vegetation, snow cover, and soil thermal properties. At the expense of reduced physical significance, these analytical models allow the estimation of long-term SFG depth at large spatial scales (e.g. Peng et al 2017, Qin et al 2018. However, both numerical models and analytical solutions described above require intensive parameterization and numerous measurement data for model calibration, which limits their application in data-scarce regions such as the QTP. Hence, the aim of this study is to develop a simple yet robust approach to estimate the long-term spatially resolved patterns of the SFG depth in cold regions. We utilize a semi-empirical ground temperature (GT) model using ground surface (0-5 cm below the surface cover) temperature as inputs to estimate the vertical GT profile and then derive the SFG depth as the maximum soil depth with temperature lower than 0 • C. In the following, we denote this approach of estimating SFG depth from GT as the 'GT-SFG depth' model. The paper is structured as follows: first, we introduce the study area and describe the GT-SFG depth model in section 2; next, we calibrate this GT-SFG depth model using measured SFG depth data at the meteorological stations; then, we apply this approach to the entire YZRB using spatially resolved data assimilation products (GLDAS, Global Land Data Assimilation System); and lastly, we estimate future trends of the SFG depth in the YZRB under projected climate change scenarios from CMIP5 (the 5th phase of the Coupled Model Intercomparison Project, Taylor et al 2012).

Study area and data
The study area, the YZRB (82 • 00 ′ -97 • 07 ′ E, 28 • 00 ′ -31 • 16 ′ N, figure 1), is located on the northern slope of the Himalayas with elevations ranging from 149 to 6963 m (Zhong et al 2014, Zeng et al 2018. The upstream and midstream basins (∼70% of the entire basin area) are mostly covered with sandy loam (Zeng et al 2018), while loam dominates the downstream basin with alluvial soil distributed along the river. In total, these two soil types cover 96% of the study area. Most of the soils are seasonally frozen, accounting for about 70% of the entire study area (Zou et al, 2017, figure 1). Meadow is the dominant vegetation type in the SFG zone (upstream and midstream basins), whereas forest mainly occupy the downstream perennial unfrozen region in the southeast (Zeng et al 2018). The mean annual total rainfall amount (based on the rainfall data from 1980 to 2010) gradually decreases from southeast to northwest. The mean annual total rainfall amount in the southeast can reach as high as 1600 mm, whereas in the northwest, the mean annual total rainfall amount rarely exceeds 300 mm. Due to the Himalayan collision and post-collision, the YZRB is tectonically active with three geological structures: Himalayan orogen, Gangdese arc, and Lhasa terrane. The northeastern basin is interspersed by the Bomi-Chayu batholith (Yin 2006). The complex geological structures, lithology, and thermal conditions contribute to numerous geothermal fields in the east of the basin, which may affect the thermal condition and thus the thickness of frozen ground (including SFG and permafrost) in this region (Geothermal and Geological Institute, 1991). The spatial distributions of the topographic elevation, soil type, vegetation, and mean annual total rainfall amount are presented in supplementary information S1 (available online at https://stacks.iop.org/ERL/15/104081/mmedia).
A total of 10 meteorological stations (figure 1) are located within the YZRB. All stations have recorded daily GST data for the period from 1980 to 2013 (China Meteorological Data Service Center, http://data.cma.cn/en). However, the annual SFG depth data, i.e. the maximum frozen ground depth throughout a year, are only available at five stations (Dangxiong, Jiangzi, Rikaze, Nimu, Lhasa, blue triangle marks in figure 1, Climate Center of Tibet). Hence, we use all datasets at the ten stations to adjust the bias associated with the grid-data products (GLDAS and CMIP5). The five stations with SFG depth data were used to calibrate our GT-SFG depth model (described in section 2.2). Because the annual snowfall amount is marginal compared with rainfall at these stations, the snowfall is not considered in the model calibration. The details of the five stations including location, elevation, soil and vegetation types, mean annual GST, mean annual total rainfall and snowfall amount, calibrated thermal diffusivity (D H ), and simulation time period are listed in table 1.
To assess the spatial pattern of SFG depth across the entire basin in the historical period (from 1980 to 2010), we employ the spatially resolved daily surface skin temperature (the temperature on the top canopy layer, Bense et al 2016, Luo et al 2018) products from GLDAS CLSM (Catchment Land Surface Model, Rodell et al 2003) to drive the GT-SFG depth model. The GLDAS dataset has been successfully applied to the QTP for hydrological simulations (Qi et al 2018). The utilized GLDAS data have a spatial resolution of 0.25 • × 0.25 • and span the time period from 1980 to 2010. To project the SFG depth in the study area under future climate change, we opt to use the near-surface air temperature (air temperature measured at a screen-height of 1.5-2 m) products generated by the CCSM4 model developed by the National Center for Atmospheric Research (NCAR) of the United States from the CMIP5 project as the input for the GT-SFG depth model. The CCSM4 products have a spatial resolution of 1.25 • × 0.94 • . This model is selected because it provides long-term  climate projection up to the year of 2300. We consider two typical climate projection scenarios in this study-'RCP 4.5' and 'RCP 8.5' (RCP stands for Representative Concentration Pathways, Moss et al 2008). Because our GT-SFG depth model (as described later in section 2.2) uses GST as model inputs, we adjusted the GLDAS surface skin temperature and CMIP5 near-surface air temperature to approximate gridded GST based on GST measurement before using them to drive the GT-SFG depth model (for details of the data pre-processing, see supplementary information S2).

The GT-SFG depth model
The GT-SFG depth model estimates SFG depth from the temporal dynamics of GTs at different soil depths. Annual GSTs are controlled by the energy inputs to the soil surface in response to natural periodic changes in solar radiation. We first employ a sinusoidal model (equation (1) where A 0 [ • C] is the amplitude of the sinusoidal curve and and T 0min [ • C] representing the maximum and minimum daily GST; P [day] is the period of the sinusoidal function and is defined as 365 days in this study; φ 0 [day] is the initial phase to run the seasonal temperature changes. The amplitude of the temperature cycle decreases with increasing depth as less energy is available to heat the soil at deeper locations. Because it takes time for heat to travel into and out of the soil, there is a time lag C between the sinusoidal temperature curves at deeper soil depth and at the soil surface, and the time lag becomes more pronounced with increasing depth (Hanks 1992). All these changes in soil temperature amplitude and time lag affected by soil depth z [m] is described by equation (2): The variables A and C are linked to depth (z) by the following formulas: The symbol τ is the time step factor (86 400 s · day −1 ). Integrating equations (3)-(4) to equation (2) yields: This model uses daily GST as inputs for equation (6) and computes the daily GTs at different depths. We then derive the GT isotherm plots for each year and compute the SFG depth as the maximum depth of the 0 • C isotherm for each year. In addition to the temperature inputs, soil thermal diffusivity D H largely determine the SFG depth (a sensitivity analysis of the SFG depth to D H is shown in supplementary information S3). In the following analysis, we first calibrate the GT-SFG depth model for each individual station with measured SFG depth data, then apply this model to the entire YZRB using GLDAS and CMIP5 data as model input. It is noteworthy that this model is not applicable to permafrost due to different heat transport mechanisms in SFG and permafrost: only one-side freezing from surface downwards in SFG but two-side freezing from surface downwards and from permafrost table upwards in permafrost (Thomas et al 2009).

Soil spatial variability
We incorporate the spatial variability of soil thermal diffusivity D H into the GT-SFG depth model by establishing the relationship between D H and soil wetness condition for different soil types based on the calibrated D H and the local conditions at the meteorological stations. For soil wetness, as soil moisture is often considered to be among the most difficult hydrologic component to determine and varies drastically with soil depth, we use the easyto-measure mean annual total rainfall amount as a proxy to soil wetness condition (previous studies have demonstrated that soil moisture and rainfall are closely correlated in this region, e.g . Bi et al 2016). The soil wetness is thus represented by the linearly scaled mean annual rainfall (R) ranging from 0 to 1, with larger scaled R indicating wetter soil and vice versa. Loam and sandy loam are the two dominant soil types covering 96% of the study area (supplementary information S1), hence we simplify the soil types in the basin into these two major categories. Based on the calibrated D H and local conditions at the five meteorological stations, a linear relationship between soil thermal diffusivity and soil wetness is derived for loam and sandy loam, respectively (for details, see supplementary information S4). We then apply the linear models to the entire basin according to the site-specific soil type and scaled rainfall amount (R).
Note that soil moisture is also affected by soil type. Due to different soil water retention characteristics, soil type affects the rainwater infiltration, soil water holding capacity, and the drainage dynamics during and after rainfall. Under similar wetness condition, rainwater infiltrates faster into coarse soil with larger infiltration capacity than fine soil and thus increases the soil water content. But immediately after rainfall, coarse soil also evaporates and drains at higher rates, which reduce soil water content. Infiltration, evaporation and drainage jointly determine the dynamics of soil water content. These processes are affected by internal and external factors such as soil type, topography and precipitation patterns. The consideration of these detailed soil water dynamics tremendously increases the complexity of this study. Hence, we opt to use the (scaled) mean annual rainfall amount as a proxy for soil water content in this study for the purpose of regional scale study. It is also noteworthy that the SFG zone in our study area is dominated by one single vegetation type-meadow which occupies 90% of the total SFG area (supplementary information S1). This allows us to simplify the analysis by neglecting the effects of vegetation variability. But we acknowledge that vegetation type and soil organic carbon (SOC) content are important factors influencing the soil thermal diffusivity and thus the SFG depth in regions with spatially variable vegetation distribution. We provide a protocol to incorporate the impacts of vegetation type and SOC on the SFG depths into the GT-SFG depth model based on the links between soil thermal diffusivity, vegetation type and SOC. The details of this protocol are presented in supplementary information S5.

Calibrating the GT-SFG depth model
This section presents the calibration results of the GT-SFG depth model based on the five meteorological stations within the YZRB from 1980 to 2010. The calibrated D H values for the five meteorological stations are presented in table 1. A comparison of simulated and observed SFG depth at the five meteorological stations from 1980 to 2010 is shown in figure 2(a). As an example, the daily dynamics of GT at different depths (from ground surface to 2 m depth which corresponds to the maximum observed SFG depth) and the GT isotherm plots at Dangxiong station for 1980 and 2010 are shown in figures 2(b) and (c), indicating a SFG depth of ∼1.5 m and ∼0.5 m in the year of 1980 and 2010, respectively. The simulated time series of the SFG depth at the five meteorological stations shows good agreement with observations (supplementary information S6). These results indicate that, despite the simplicity, the proposed model calculates reasonable SFG depths using GST as major model input. It appears that there exists relatively large discrepancy between simulated and observed SFG depths at the Lhasa station. We attribute this discrepancy to the influence of land use: The Lhasa station is located in urban area, while the other four stations are located in natural regions covered by meadow. The urban land use type at the Lhasa station hinders the heat transport through the soil profile, leading to shallower SFG depths with smaller variability. But for the regional scale analysis (section 3.2 and 3.3), we have not specifically considered the effects of urban land use type, as the YZRB is a highly underdeveloped region with the urban land use type only occupying less than 1% of the total basin (Fang and Li 2015).

Historical SFG spatiotemporal patterns
This section investigates the spatiotemporal patterns of SFG depth in the entire YZRB for the historical time period  using (adjusted) GLDAS surface skin temperature data as model inputs. The spatial distribution of the SFG depth in 1980 (figure 3(a)) indicates a clear spatial pattern: the seasonally frozen soil is mostly located in the northwest of the basin with SFG depth ranging from 1.5 m to 4 m, whereas the SFG depth in the southeast is less than 1 m with perennial unfrozen area in the southeast (to the south of the Linzhi station). Figure 3(b) depicts three examples of the differences between the SFG depth in different years. It reflects a greater decrease in seasonally frozen soil depth in the northwest than in the southeast ( figure 3(b), a negative value denoting a decreasing SFG depth). This spatial pattern is largely due to more intense increases in the surface skin temperature in the northwest (supplementary information S7). The spatial variability of soil properties also contributed to the spatial pattern of SFG depth (higher soil thermal diffusivity in the northwest compared with south of this region, supplementary information S8). Despite annual fluctuations, the average SFG depth across the entire basin (excluding the unfrozen zone in the southeast) decreases at a rate of 2.50 cm · a −1 for the time period of 1980-2010 (shown in figure 4 in section 3.3). Our simulated SFG depths are also confirmed by the modeling results based on the semi-analytical Stefan's method that relates the SFG depth to GST and soil thermal properties. Details of the Stefan's method and its comparison to our GT-SFG depth model are shown in supplementary information S9.

Projected SFG dynamics under climate change
We employed the (adjusted) near-surface air temperature data generated by the CCSM4 model in two projected climate change scenarios (RCP 4.5 and RCP 8.5) as model inputs for the GT-SFG depth model to estimate future changes of the SFG depth in the YZRB. The simulated basin average SFG depths show continuing negative trends in both RCP 4.5 and RCP 8.5 climate scenarios (figure 4). The shaded blue and red areas in figure 4 shows the ranges of average SFG depth under RCP 4.5 and RCP 8.5 scenarios for the period 2006 to 2300 accounting for the uncertainties of soil properties as described in supplementary information S4. Particularly, we observed a linear decrease of the average SFG depth (log scale) with time in the RCP 8.5 scenario (log10 (SFG) = 0.23-0.0047 · T, R 2 = 0.85, with T denoting the elapsed years since 1980) and that the average SFG depth tends to level off at ∼0.1 m after 2180 in the RCP 8.5 scenario, implying that most areas in current SFG zone are projected to no longer have SFG in the future due to climate warming. The basin average of the minimum near-surface air temperature throughout a year is projected to exceed 0 • C after 2180 in the RCP 8.5 scenario (see the inset of figure 4). In addition to the CCSM4 model, we also computed the future SFG depth using near-surface air temperature products generated by 8 additional GCMs (General Circulation Models) from CMIP5. The results confirm the negative trend of the SFG depth, despite the large range of the SFG depth computed by different GCMs (supplementary information S10).

Summary and conclusions
We propose a semi-empirical GT-SFG depth model to investigate the spatiotemporal dynamics of the SFG depth in the YZRB on the QTP via the simulation of long-term time series of GT at different soil depths. This model was first calibrated using ground station observations and then applied to the entire basin using data assimilation products of GLDAS (1980-2010) and GCM projections (2006-2300 to assess the spatiotemporal patterns of the SFG depth under climate change. The results suggest that the average SFG depth in YZRB has decreased at a rate of 2.50 cm · a −1 from 1980 to 2010, with more severe decrease in the northwestern upstream regions than in the southeastern downstream regions. The future projections indicate that the current SFG depth in the YZRB will continue to decrease in response to future warming and may disappear by 2180 under the RCP 8.5 scenario (if not considering the transition of permafrost to SFG).
The proposed modeling framework provides a simple yet robust approach to estimate the temporal dynamics of vertical temperature distribution and SFG depth. This model is particularly useful in datascarce regions (such as the QTP in this study) as it requires minimum model input (i.e. the surface GT, which is nearly universally available over the globe) compared with other physically-based models. In addition, the reduced computational burden allows the estimation of the depth at the regional scale (e.g. over 10 000 km 2 ) spanning long time periods (e.g. decades). The temporally evolving spatial distributions of the SFG depth at the regional scale will be useful for constraining more sophisticated physically based numerical models that consider specifically water flow and heat transport in multi-phase media in cold regions. More field data with increased spatial density for a longer time period will further constrain our model, but our estimates provide an important basis for the evaluation of the hydrological cycles (e.g. seasonal groundwater recharge and surface water-groundwater exchanges) in cold regions under changing climatic conditions. the State Environmental Protection Key Laboratory of Integrated Surface Water-Groundwater Pollution Control of China. We are grateful to David N Lerner, Yuqing Feng and Yu Feng for helpful discussion and assistance.