Spatiotemporal variations and its driving factors of ground surface temperature in China

The ground surface temperature (GST) serves as a crucial indicator for understanding land-atmosphere mass and energy exchange. The shift from manual measurement to automated station for GST in China after 2002 introduced inconsistencies at certain stations, potentially distorting research findings. Here, daily automatedly observed GST from 2003 to 2017 at 615 selected meteorological stations were updated by constructing linear regression model based on manually observed air temperature (AT) and GST from 1960 to 2002. Then, the spatiotemporal variations of GST from 1960 to 2017 and its driving factors were investigated. Results indicated that: (1) the AT-GST linear regression model could effectively mitigate the inconsistency caused by the change of GST observation methods, enhancing data reliability. (2) GST in China showed little change from 1960–1980, but increased significantly across all regions from 1980 to 2000, with the increase rate slowed down except in the Qinghai–Tibet plateau (QTP) and southwest China after 2000. Notable GST increase is concentrated in colder regions, including the QTP, northeast (NEC), and northwest China (NWC). (3) Evapotranspiration (ET) and vapor pressure deficit were the primary drivers of annual GST variations at the regional scale, while their contributions to GST variations exhibited notable seasonal variability. Our findings could offer valuable scientific insights for addressing climate change, enhancing surface environmental models, and safeguarding ecological environments.


Introduction
The earth's surface serves as a vital boundary interfacing the atmosphere, biosphere, hydrosphere, and lithosphere.The ground surface temperature (GST), representing the earth's surface temperature, is an important variable to characterize mass and energy exchange between land and atmosphere, with applications spanning climate, ecology, hydrology, and agriculture (Jenerette et al 2007, Luo et al 2015, Le et al 2021, Șerban et al 2023), Tang et al (2007).GST generally refers to surface temperature below 0-5 cm from the ground Luo et al (2015) (specific depth depends on research purpose) which is a comprehensive reflection of various factors such as soil temperature, grass temperature, bare soil temperature (Luo et al 2018).GST data is essential for land surface models (Kalma et al 2008), permafrost distribution models (Gao et al 2023), and permafrost thermal simulations (Westermann et al 2013), playing a vital role in estimating and predicting permafrost distribution and active layer thickness.
Currently, some researchers investigated GST variations and its driving factors based on in-situ GST observations across different regions, with the goals to understand the land-atmosphere interaction and flux transfer processes.Beltrami et al (2003) reconstructed GST histories (1950reconstructed GST histories ( -2000) ) from temperature versus depth profiles measured at 246 stations distributed across Canada and found significant spatial variability of GST, indicating that the greatest warming occurs in southern areas of Canada. Reiter (2007) discussed the variability of recent GST at four stations over a large area within the Albuquerque basin, central New Mexico from 1971-2000.Șerban et al (2023).The monthly GST series from 1900 to 2014 for five sites of QTP were generated using linear calibration models and validated through three other sites using the same methods, demonstrating high correlations (r > 0.90) in capturing the dynamics of in GST for all earth system models data (Xing et al 2023).Most researches focused on the spatiotemporal changes and driving factors of land surface temperature (LST) based on remote sensing data, while there are few studies discussing the driving factors of measured in-situ GST in China due to lacking the longterm and consistent in-situ GST observations.Tian et al (2022) used the geographical detector method to quantitatively analyze the driving factors such as sunshine hours, precipitation, elevation, and NDVI on LST changes in China from 2001 to 2020.Peng et al (2018a) analyzed the dominant factors of LST changes in different seasons in Shenzhen, China using regression analysis.
Long-term in-situ GST observations with good data quality and consistency is essential for reliably investigating the spatiotemporal variations of GST and its driving factors.In China, manual observed data for GST has been available since the 1950s.After 2002After , especially in 2004After -2005, there was a widespread transition from manual to automatic observations for GST.Automatic stations measure GST via platinum resistance ground temperature sensor, while manual stations use curved tube geothermal meter to measure GST.This change of GST observation methods leads to the discontinuities in the GST observations, particularly in the north of China during winter when there is snow cover (figure S1) (Liu et al 2008, Cui et al 2021).The automatic stations usually measure higher GST compared to manual observations.While disparities between automated and manual observations are expected, the changes of GST observation methods could lead to data inconsistencies.
Given the notable variations in measurement principles between the manual and automated observation methods, it is essential to conduct regular consistency checks and make necessary adjustments when utilizing the observations.Currently, there are three main approaches to reduce the uncertainties caused by the changes of GST observation methods: (1) Station or period exclusion: conduct consistency checks to remove the stations with significant discontinuities or the periods after the change in GST observation methods (Qiao et al 2015, Liao et al 2019).However, reducing stations could significantly reduce the total number of stations, limiting the calculation of averages and the accuracy of spatial interpolation.Furthermore, avoiding the analysis of data after the change in observation methods would result in a lack of understanding of the post-change period.
(2) LST substitution: substitute GST with remotely sensed LST data (Peng et al 2018b, Song et al 2021, Zhou et al 2022).Remotely sensed LST data generally starts from the 2000s, which means past several decades before 2000s cannot be analyzed.Moreover, GST and LST have different physical meanings, and the significance represented by their results is also different when analyzed (Luo et al 2018).
(3) Correction model: The third method is to use statistical analysis or machine learning methods to correct GST observations due to the change of GST observation methods (Du et al 2020).Cui et al (2021) analyzed the relationship between snow surface temperature (SST) and GST to update GST observations removing the influence of GST observation methods.Daily air temperature (AT), sunshine duration, precipitation, wind speed, relative humidity, solar radiation, and diurnal temperature range was used to train the machine learning model to correct the daily GST observations after 2003 by Liu et al (2022).However, some potential confounding variables beyond SST need to be investigated in correcting GST observations, and the effectiveness of machine learning models for GST correction may be constrained by the complexity of surface temperature dynamics and data limitations.
In this study, we aim to obtain accurate, reliable, and continuous long-term GST observations, and to explore the spatiotemporal patterns of GST and its dominant driving factors.Two questions will be answered: (1) how to build the simple but reliable GST correction model to remove the uncertainties from the changes of GST observation methods?(2) what are the spatiotemporal patterns of GST and its dominant driving factors?This study provides a scientific basis for the study of GST evolution and its environmental effects, as well as provides guidance for policy makers to deal with global warming.

Data source
Daily GST and AT from 1960 to 2017 of 615 weather stations were provided by China meteorological administration (CMA).The spatial distribution of selected stations, along with seven geographical divisions (i.e. the Qinghai-Tibet plateau, QTP; northeast China, NEC; northwest China, NWC; north China, NC; the Loess Plateau, LP; southeast China, SEC; southwest China, SWC) of China (figure 1).
Monthly gridded precipitation (Prep) dataset was sourced from ERA5_land reanalysis data.Monthly root soil moisture (SM) and evapotranspiration (ET) were downloaded from Gleam 3.6a.Monthly gridded vapor pressure deficit (VPD), surface solar radiation (Rs) and snow depth equivalent (SWE) were obtained from TerraClimate.Monthly GIMMS leaf area index (GIMMS LAI4g) was used as the vegetation data.Annual 500 m irrigated cropland maps across China were obtained from the hybrid irrigation dataset (Zhang et al 2022).The 30 m annual land cover dataset and its dynamics in China from 1990 to 2019 was provided by Yang and Huang (2021).Here we classified the landcover data into three classes including vegetation cover, urban, and water, and calculate their areas.For simplicity, we used Irri, VA, UA and WA to denote the irrigation areas and the areas of these three types of landcover, respectively.Given the small variation in our LUCC data until the 1990s (Liu et al 2014), the 1985 LUCC data is used to replace the LUCC data from 1981 to 1984, and the 1990 LUCC data is used to replace the LUCC data from 1986 to 1989 to ensure the consistency of the LUCC data with the NEP data in the time series.Similarly, we also consider the irrigation area before 2000s remain unchanged.To maintain unity with the time range of vegetation index, all data were extracted from 1982-2017 for attribution analysis in section 3.3, and all the datasets were resampled to 0.1 • × 0.1 • for our analysis.All datasets used in this study were summarized in table 1.

Correction models of GST
In this paper, 0 cm GST was selected to analyze the spatiotemporal variation of GST.However, after 2002, China established an automatic weather station network to obtain official records of ground meteorological data (Liu et al 2008).This system replaced over 50 years of manual observations and led to significant changes in observation instruments and equipment (Cui et al 2021).The change in observation methods has resulted in the inconsistency in the GST sequences, around 2003 (figure S1).The observed values of AT were consistently continuous, consistent, and undisputed between 1960 and 2017, which were not been affected by the changes in monitoring methods (Du et al 2020).Moreover, we have discovered a strong linear relationship between monthly GST and AT (figure S2).Therefore, a monthly linear regression model was created using raw GST and AT data from 1960 to 2003 (equation ( 1)).First, the monthly linear regression model between raw GST and AT during the period from 1960 to 2002 was established.In this way, the slope (K) and intercept (B) of linear regression equations for GST and AT could be calibrated monthly at each station.After ensuring the reasonable and reliable relationship between AT and GST, the predicted GST during the period from 2003 to 2017 could be obtained based on the calibrated monthly linear regression models and observed AT during the where: GST m,n was the corrected monthly GST in month m at station n; AT m,n represented the monthly AT in month m at station n; K m,n and B m,n were the monthly slopes and intercepts of the linear relationship between GST and AT.Nash efficiency coefficient (NSE), determination coefficient (R 2 ), mean average error (MAE) and root mean square error (RMSE) were selected as the evaluation criteria of the developed linear regression models between GST and AT (Moriasi et al 2007).

Spatiotemporal variation of GST
The updated monthly GST data was used to interpolate the whole study area (0.1 • × 0.1 • ) via the Kriging method (Oliver and Webster 1990), which is often used for spatial interpolation of temperature (Holdaway 1996, Wu andLi 2013).
The climate tendency rate its significance test of each grid point GST from 1960 to 2017 was calculated by the Theil-Sen median method (also called Sen slope estimation) (Sen 1968) and Mann-Kendall trend test method (Mann 1945, Kendall 1948), and the spatial distribution of climate tendency rate was analyzed.The Sen's slope equation for a number of N data sample pairs is shown as follows: where: x i and x j are the data series at time j and i (j > i), respectively.The N values of Q i are ranked from smallest to largest and the median of Sen's slopeis computed (Q med ) as: In terms of interannual aspects, this study divided 58 years from 1960 to 2017 into three stages, namely 1960-1979, 1980-1999, and 2000-2017.We investigated the climate tendency (slope) of GST in each stage of 7 regions (NEC, NWC, QTP, LP, NC, SEC, and SWC) in China, in order to analyze the changes in different regions at different stages.In terms of the seasonal aspects, this study analyzed the changes in GST during different seasons from 1960 to 2017.

Canonical analysis
Canonical analysis (called constrained ordination) was used to assess the driving factor of GST variations based on the relative contribution.Here, we regarded GST as the response variable and other climate variables as mentioned above including Prep, ET, Rs, SM, SWE, VPD, LAI and Irri, VA, UA, and WA as predictors.All these analysis were implemented in the 'rdacca.hp'R package developed by (Lai et al 2022)

Evaluation of corrected GST observations
Given the widespread shift in GST observation methods after 2002 across most weather stations, we opted to utilize GST data spanning from 1960 to 2002 as the training dataset for the AT-GST linear regression models and conducted the cross validation.The calibration phase was established using even years within the period from 1960 to 2002, while the validation phase encompassed odd years from the same temporal span.In figure 2, throughout the calibration period, the MAE consistently remained below 0.8 • C for nearly all stations, with median and mean values both below 0.6 • C. The RMSE was slightly higher than MAE, yet its median and mean values stayed below 0.7 • C. Notably, both the median and mean values of NSE and coefficient of determination (R 2 ) surpassed 0.99, while the minimum NSE values were above 0.96, and the minimum R 2 values exceeded 0.98.
Similar trends were evident in the validation period, where the majority of stations reported MAE and RMSE values ranging from 0.4 to 0.6 • C and 0.4 to 0.8 • C, respectively.The median and mean values of MAE and RMSE remained consistently below 0.8 • C. Furthermore, NSE and R 2 values for most stations exceeded 0.99, with the lowest values not falling below 0.95.Concurrently, a sample station was selected from each of the seven distinct regions to compare the predicted GST sequences with the original GST sequences spanning from 1960 to 2017 (figure S1).The analysis revealed that the association between the predicted GST sequences and AT remained relatively consistent prior to 2002.In contrast, the raw GST sequences displayed a tendency to diverge after 2002, especially in the northern regions such as the NEC, and NWC.These findings confirmed the reliability of the monthly linear relationship between AT and GST, thereby substantiating the viability of utilizing updated GST sequences derived from predicted data spanning from 2003 to 2017.
The 95% confidence interval of the regression coefficient of the GST-AT regression model is used to describe the uncertainty of the regression model (figure S3).The narrow confidence interval demonstrated the reliability of the linear regression model of the GST-AT model.In addition, the parameters had the temporal and spatial heterogeneity.For the temporal heterogeneity, the parameter uncertainties during the growing seasons (May-October) were larger than those during the non-growing seasons due to the impacts of vegetation.For the spatial heterogeneity, the parameter uncertainties in the QTP and SEC were larger than those in other regions, as some stations in the QTP and SEC had some data gaps in the observations.

Spatiotemporal variation of GST in China
Mean annual GST in China during the period from 1960 to 2017 varied significantly from approximately −5 • C to 30 • C, with pronounced spatial disparities (figure 3(a)).The maximum GST was found in the SEC and SWC, where over 80% of the regions exhibited a mean annual GST above 15 • C. Regions in the NC, LP, and NWC experienced GST between 10 and 20 • C. Conversely, the lowest GST was observed in the NEC, QTP, and Altai and Tianshan mountains of NWC, at the majority of locations of where GST values ranging from 0 to 10 • C. Particularly in the northernmost part of the NEC, some regions had an average annual GST below 0 • C. Overall, driven by the influence of solar radiation (Karimi Firozjaei et al 2022), the distribution of GST exhibited a gradual decrease from south to north, corresponding with increasing latitude.Elevation was also found to impact the distribution of GST, with lower elevations correlating with higher GST and plains exhibiting higher GST compared to mountainous and plateau areas (Allevato et al 2019).The QTP boasted high elevation, while the NEC possessed high latitude, and the Altai and Tianshan Mountains of NWC featured both high latitude and high elevation.Consequently, these three regions have the lowest GST levels in China.
Figure 3(b) displays the distribution of climate trend rates (slopes) of China's mean annual GST from 1960 to 2017.According to the results of the MK trend test, approximately 5% of the Chinese regions showed a significant increase in GST, while 5% of the regions showed no significant trend.No region showed a significant decrease in GST.Except for some places in SWC, where the GST decline rate is about 0-0.2 • C/10 yr, GST in most regions of China showed a rising trend.Hotspots were found in the central QTP, as well as in the north part of NEC and NWC, with the highest slopes ranging 0.5-0.8• C/10 yr.Overall, the climate trend rates in China primarily exhibited an increasing pattern from south to north.
In order to examine changes in interannual variation of GST, the study period was divided into three stages for analysis : 1960-1980, 1980-2000, and 2).From 1980 to 2000, the climate trend rates in all regions increased sharply.The NEC, LP, and NWC had the fastest increasing rates of 6.87, 6.16, and 5.18 • C/10 yr, respectively, while SWC had the slowest increasing rate of 2.00 • C/10 yr or above.After 2000, compared to the period of 1980-2000, the temperature increasing trend in most regions slowed down, with the exception of SWC and QTP which experienced an accelerated warming rate.During the 1980s and 1990s, global industrialization and human activities resulted in a surge of greenhouse gas emissions, leading to a rapid rise in global temperatures.However, since 2000, the implementation of emission reduction measures and clean energy policies by various countries has helped to slow down the upward trend in temperature.
In addition to interannual differences, China's GST exhibits significant regional variations across four seasons, including spring, summer, autumn, and winter (figure 5).The warming centers in the QTP, NEC, and NWC were consistent across all seasons, with particularly strong increasing trend in spring, autumn, and winter.In contrast, there was a two-level differentiation in summer, with the NWC experiencing the highest warming rate of 0.9 • C/10 yr, while the transitional zone between the NC, LP, SEC, and SWC had a small-scale cooling trend of −0.3 • C/10 yr, which was statistically insignificant.Overall, most regions of China showed the greatest warming rates in spring and winter, with the QTP, NEC, and the northern NWC being key areas of warming.For the averaged change rates of GST in seven regions during four seasons (table 2), maximum change rates in the NEC, NWC, QTP, LP, and NC occurred in spring (0.46 • C/10 yr), spring (0.45 • C/10 yr), winter (0.37 • C/10 yr), spring (0.40 • C/10 yr), and spring (0.40 • C/10 yr), respectively.In contrast, the SEC and SWC had the lowest (or second lowest) change rates in almost all seasons.Although the LP did not experience extreme cold or hot conditions in any season, GST in most areas of the region showed significant increasing trend, with change rates of 0.40 • C/10 yr in spring and 0.37 • C/10 yr in winter.Overall, the average change rates of GST in all regions were greater than 0 • C/10 yr in all four seasons.Except for the SWC, which had the highest change rate in autumn (0.19 • C/10 yr), slightly higher than in winter by 0.01 • C/10 yr, the change rates in other regions were generally higher in spring and winter than in summer and autumn.

Dominant driving factors for the spatiotemporal variations of GST
Table 3 unveiled that total contributions of eleven environmental variables to annual GST changes were 63.6%, with no significant differences in cumulative contributions from climate-, vegetation-, and LUCCrelated variables.Seasonally, environmental variable contributions varied, peaking at 71.2% in summer and dropping to 23.4% in winter.Vegetationrelated variables had higher contributions in summer and autumn, reflecting the growth and senescence periods.Winter, characterized by lower temperatures and dormant vegetation, had minimal contributions.Notably, VPD had the most substantial influence on GST variations, contributing 15.2% annually and 13.7% in summer.Variables like VA, (ET, and precipitation were key drivers in spring, autumn, and winter, contributing 10.8%, 10.1%, and 9.1%, respectively.These findings emphasize the X Gao et al Figure 4. Temporal variations of GST across three stages (1960-1979, 1980-1999, 2000-2017) in seven regions China.
impact of seasonal climatic conditions on GST variance.Despite VPD and ET dominance in summer and autumn, land cover-related variables, including vegetation area (VA), water area (WA), and urban area (UA), made discernible contributions (figure 6), adding significance to the broader context of GST fluctuations.Urban expansion enhances the heat island effect that regulating the balance of heat and energy on land surface, resulting in changes in diurnal LST and surface-AT, thus extensively contributing to the increase in LST (Bian et al 2017, Moazzam et al 2022).Our findings align with the previous conclusion that changes in vegetation have had a significant impact on China's regional temperature (Yu et al 2020).Figure 6 indicated the dominated factors among eleven environmental variables contributing to seasonal (figures 6(a)-(d)) and annual (figure 6(e)) GST variations in seven subregions (i.e.LP, NC, NEC, NWC, QTP, SEC, and SWC).LP, QTP, SEC and SWC collectively accounted for over 80% -of the overall annual GST variations, driven by distinct factors in each region.Notably, VPD influenced 23.6% in LP, WA contributed 15.2% in QTP, ET accounted for 14.4% in SEC, and SWE held about 17.1% in SWC.These results are highly aligned with previous research that evapotranspiration was the key factor affecting daytime LST in south China (Song et al 2021).In NC, NEC, and NWC, total contributions ranged mainly from 60% to 70%, with VPD being the predominant factor, contributing 9.8%, 25.5%, and 22.5% in NC, NEC, and NWC, respectively.Regional variations in seasonal GST variations consistently maintained levels within the range of 70%-80% across seasons, except for specific regions in autumn (LP, NWC, NC) and winter (NWC) (table 3).Noteworthy was the elevated contribution observed in LP, NC, and NEC in spring, while other regions generally peaked in summer.In LP, NC, and NEC areas, climate-related factors primarily influenced GST variations across different seasons.In the NWC, climate indicators were dominant, except in autumn, where LUCC-related indicators played a role.In the SWC, vegetation-related indicators controlled GST variations in autumn and winter.The QTP and SEC showed distinct seasonal patterns in GST variations, with VPD having a significant influence in spring, and ET being prominent in autumn and winter, especially in QTP, NC, and SWC.Surprisingly, NEC region experienced VPD as the main driver of seasonal GST variations.In regions less affected by ET, UA, or VPD, Rs or SWE had notable contributions to seasonal GST changes.For example, ET contributed only 1.5% to winter GST variations, while Rs accounted for 25.1% in the NWC region.Conversely, changes in Irri, Prep, and SM had negligible effects on GST variations at both national and regional scales.
Figure S4(A) delineated discernible patterns in annual GST from 1982 to 2017 for China as a whole.Notably, GST generally demonstrated significant   (Peng et al 2014).
The interaction between decreased surface albedo and increased ET from vegetated areas engendered significant alternations in GST.This was attributed to the cooling effect of evaporation and the diminished albedo linked to verdant canopies, which absorbed greater quantum of solar radiation, leading to an

Conclusions
This study aims to correct the in-situ GST observations post-2002 due to the changes in GST observation method, utilizing the strong correlation between GST and AT.Monthly linear regression models between AT and GST were proposed and found to be a simple yet effective method for mitigating uncertainties caused by changes in observation methods.Then, we analyzed the temporal and spatial variations of GST (1960GST ( -2017) ) and investigated its driving factors  across China.The main conclusions are summarized as follows: (1) Correcting GST observations using AT-based models proved effective in reducing uncertainties arising from observation method changes after 2002.
(2) GST in China exhibited significant spatiotemporal variations.There was no significant change in GST from 1960 to 1980, while the increase rate accelerated sharply from 1980 to 2000.After 2000, there was a slight decrease in the upward trend in some regions, with the most significant warming observed in QTP, NEC, and NWC, particularly during spring and winter.Some areas experienced a maximum change rate of 1.2-1.6 • C/10 yr.(3) ET and VPD emerged as dominant driving factors for national GST variations with high partial correlations.Largescale afforestation programs influenced GST through changes in evapotranspiration and surface albedo.Climate-related indicators played a significant role in summer and autumn, while their contribution was minimal in winter.LP, NC, and NEC experienced dominant influence from climate-related indicators, while QTP and SEC regions showed distinct seasonal characteristics influenced by various indicators.VPD was identified as the dominant factor in spring and summer GST changes, except for QTP and SEC regions, where ET played a significant role in autumn and winter GST variations, except for LP, NEC, and NWC regions.
This research contributes valuable insights into the spatiotemporal variations of long-term GST in China and can inform climate change mitigation efforts in the region.Due to the limitations of data sample and model structure, there are uncertainties in the new GST sequence, which have a certain influence on the spatiotemporal variation of GST.Future studies will employ additional methods and consider more factors to improve our understanding of GST change mechanisms.

Figure 1 .
Figure 1.The location and distribution of weather stations in China.NEC stands for northeast China; NWC stands for northwest China; SEC stands for southeast China; SWC stands for southwest China; LP represents loess plateau of China; QTP stands for Qinghai-Tibet plateau of China.

Figure 2 .
Figure 2. Box diagram of GST predicted model evaluation indicators (MAE, RMSE, NSE, R 2 ) in calibration and validation periods.The scatter points on the left of each box diagram are the distribution of each indicator.

Figure 3 .
Figure 3. Spatial distribution of mean annual GST over multiple years and its change rate in China from 1960-2017.(a), Mean annual ground surface temperature.(b), Change rate ground surface temperature.The black grid indicates that the slope passes a significance test with a P-value of 0.01.

XFigure 5 .
Figure 5. Change rates of GST across China during the period from 1960 to 2017.(a), Spring.(b), Summer.(c), Autumn.(d), Winter.The black grid indicated that the slope passed the significance test with P-value of 0.01.

Table 1 .
Summary of data used in this study.

Table 2 .
Average change rates of GST in different regions of China from 1960 to 2017 ( • C/10 yr).

Table 3 .
The total explained variance in GST variations.

Table 4 .
(Dai 2013, Grossiord et al 2020)tion of the greatest dominated variable to GST.Grossiord et al 2020), thereby establishing a pathway through which VPD can influence GST via the coupling of moisture and temperature dynamics(Dai 2013, Grossiord et al 2020).Meanwhile, the amplification of VPD substantially impacts vegetation transpiration processes (Yuan et al 2019, Lu et al 2022), which in turn exerts an indirect influence on GST alterations.As a result, VPD played the dominant role in GST changes for the NEC and NWC (table4).