Abstract

Land surface temperature Ts and near-surface air temperature Ta are two main metrics that reflect climate change. Recently, based on in situ observations, several studies found that Ts warmed much faster than Ta in China, especially after 2000. However, we found abnormal jumps in the Ts time series during 2003–05, mainly caused by the transformation from manual to automatic measurements due to snow cover. We explore the physical mechanism of the differences between automatic and manual observations and develop a model to correct the automatic observations on snowy days in the observed records of Ts. Furthermore, the nonclimatic shifts in the observed Ts were detected and corrected using the RHtest method. After corrections, the warming rates for Ts-max, Ts-min, and Ts-mean were 0.21°, 0.34°, and 0.25°C decade−1, respectively, during the 1960–2014 period. The abnormal jump in the difference between Ts and Ta over China after 2003, which was mentioned in existing studies, was mainly caused by inhomogeneities rather than climate change. Through a combined analysis using reanalyses and CMIP5 models, we found that Ts was consistent with Ta both in terms of interannual variability and long-term trends over China during 1960–2014. The Ts minus Ta (TsTa) trend is from −0.004° to 0.009°C decade−1, accounting for from −3.19% to 5.93% (from −3.09% to 6.39%) of the absolute warming trend of Ts (Ta).

1. Introduction

Global mean surface temperature (GMST) change is an important metric used to measure climate change (Hansen et al. 2010; Wang and Dickinson 2013). However, different types of measurements are typically used over land and ocean. For most datasets used to evaluate the GMST (see Table S1 in the online supplemental material), the surface temperature over the ocean is represented by sea surface temperature (SST) rather than marine air temperature, while the surface temperature over land is represented by near-surface air temperature Ta rather than land surface temperature Ts (Cramer et al. 2013; Jones et al. 1999; Mann et al. 2017). The response of Ts to climate change and its difference from Ta are critical factors for assessing the effects of climate change, as well as promoting the evaluation and improvement of climate model simulations (Cowtan et al. 2015).

As a key variable at the land–atmosphere interface, Ts drives the exchange of longwave radiation and turbulent heat fluxes. Meanwhile, this parameter is also important in the physical processes of the land surface and water cycle from the global to regional scale (Anderson et al. 2008; Karnieli et al. 2010). Therefore, Ts can directly characterize spatiotemporal variations in the surface equilibrium state and is widely used to validate climate models (Koch et al. 2016; Singh et al. 2016), estimate drought conditions (Mühlbauer et al. 2016; Singh et al. 2016), and depict urban heat islands (Jin 2012; Schwarz et al. 2011; K. Wang et al. 2007). Due to the larger heterogeneity of land surface properties, Ts shows a larger spatial difference and thus in situ observations are less representative of large areas (Liu et al. 2006). Therefore, unlike Ta, Ts is not a routinely measured variable at meteorological stations in most regions (Wang et al. 2014).

However, Ts is one of the most important variables of the Earth system. Many specific measurement campaigns have recognized Ts as a high-priority parameter [e.g., the First ISLSCP Field Experiment (FIFE) and the Coordinated Enhanced Observing Period (CEOP); Catherinot et al. 2011; Jin et al. 1997; Koike 2004]. In those field experiments, Ts was estimated by an infrared radiometer over a surface with known emissivity (Jiménez et al. 2012). The in situ observations of Ts have been used widely to evaluate and validate the algorithm for satellite retrievals and the surface scheme for model simulations (Bosilovich 2006; Jiménez et al. 2012; Wang et al. 2005). However, most of those observations are only available at a very limited number of stations, and the lengths of most records are only a few months. The sparse and short observations for Ts make it difficult to analyze the interannual variability and long-term trend of Ts with climate change.

As a routine observation at metrological stations, the measurement of Ts began as early as 1951 in China, and Ts has been continuously recorded at approximately 2400 stations (see details in section 2a). The Ts observations at weather stations have long been questioned for their weak spatial representations due to large land surface heterogeneities (Fuchs and Tanner 1968), but recent studies have shown that the observed Ts has good stability and spatial consistency on the yearly and monthly scales (Good et al. 2017).

Wang et al. (2017) analyzed the interannual variability and warming pattern of Ts based on the station observations from China Meteorological Administration and noted that the variability in Ts is similar to that in Ta from 1965 to 2004. Liu et al. (2017) analyzed the variabilities in Ta and Ts over the temperate steppe region of North China during the 1965–2007 period and suggested that long-term changes in Ts and Ta have significant differences in spatial patterns and seasonal contrasts. Zhang et al. (2016) argued that the soil surface temperature increased 31% more than the air temperature across China from 1962–2011, which may result in more soil carbon release to the atmosphere than predicted.

However, it is important to note that most previous studies have found that the warming rate for Ts became much faster than that for Ta after 2003, especially in the northern regions during winter (Liu et al. 2017; Wang et al. 2017; Zhang et al. 2016). For instance, Wang et al. (2018) determined that the difference in the warming rate between the surface and near-surface air daily mean temperatures was 0.05°C decade−1 (10.9% of the absolute warming trend of Ts during the same period) from 1980 to 1998 and increased to 0.92°C decade−1 from 1999 to 2011. Various explanations have been proposed for this sudden increase in Ts after 2003, such as the enhancement of snowfall and changes in total cloud amount. However, these explanations are merely speculation and do not provide convincing evidence or quantitative attribution.

In this study, we compared the annual mean difference between observations of the Ts and land–atmosphere temperature differences (TsTa) before and after 2003, evaluated the inhomogeneous observation records, and found that the transition from manual to automatic measurements during the 2003–05 period in China is the main reason for the abnormal warming in the observed Ts. During the manual observation period, Ts is a measure of the snow surface temperature; however, for automatic observation, Ts is the temperature under the snow cover when snow is present. Owing to the heat preservation effect of snow, the Ts observed under snow cover is much warmer than that observed on the snow surface, especially during nighttime.

This paper focuses on the correction methods of Ts observations under snow cover conditions during the automatic observation period. Based on the robust relationship among Ts and Ta, snow depth (SD), and surface solar radiation Rs, a new empirical method is presented and tested in various contrast experiments. Using this method, we corrected the observed Ts and built a homogenized dataset with the help of RHtestV4 software (Wang and Feng 2013). Combined with reanalysis and model simulations, we explored the variability in the difference between Ts and Ta during the 1960–2014 period and further trends during the 2015–2100 period under different representative concentration pathways (RCPs) for greenhouse gas and aerosol concentrations.

2. Data and method

a. Observation dataset

The daily observations from approximately 2400 meteorological stations across China are used in this study, which are provided by the China Meteorological Administration (CMA; http://data.cma.cn/). The raw data from the dataset have been passed through strict quality controls by the CMA, including identification of outliers, assessment of spatiotemporal consistency, and correction of suspicious and erroneous data.

As shown in Fig. S1b, the number of observation stations rapidly increased before 1960 and then stabilized. As the station observations of SD are only available before 2015, this study only focuses on the 1960–2014 period. We selected 2134 stations to be used in this paper and all of them have more than 40 years of valid observations during the 1960–2014 period (see Fig. S1a). The remaining 266 stations have too many missing values and not enough valid observations to be used in this study.

As reference variables, Ta, SD, and Rs are used to correct and homogenize the observations of Ts in this study. Therefore, it is critical that these variables are homogeneous in time. The SD is different from the other two input variables, which have always been measured manually and will not be affected by automatic transformation. However, the instruments used to measure Ta and SD changed after the transformation to automated instruments. For example, the observation equipment for Ta transformed from glass thermometers to platinum resistance thermometers. However, the physical attributes of those variables are not changed after automatic transformation, and the observations made during the manual and automatic periods appear to be consistent (Fig. S2).

However, subtle inhomogeneities may be present in these records due to changes in instrument and observation frequencies. As shown in Fig. S3, there is no significant difference in the probability density distribution for observed Ta and Rs between manual and automatic observations. Therefore, any possible heterogeneity problems can be corrected by the RHtest method, which has been widely used in previous studies to homogenize climate variables. Hence, we homogenized the observations of Ta, SD, and Rs based on the RHtest software used in previous studies (Sun et al. 2016; Xu et al. 2013; Zhou et al. 2018). A detailed description of the homogenization process is provided in the supplementary material.

Surface temperature Ts is a routine observation variable at meteorological stations in China. The observation field for Ts is bare soil at meteorological stations, and the surfaces must be kept loose, grassless, and flat. Before 2003, Ts was measured by glass thermometers, including a common thermometer, a maximum thermometer, and a minimum thermometer. The readings of the thermometers were recorded manually. The thermometers were placed horizontally on the observational field; half of each thermometer was buried in the soil, and the other half was exposed to the air; thermometers were checked by observers daily. The buried part of the thermometers had to be in close contact with the soil without gaps, and the exposed part had to be kept clean from dust and dew. When the surface was covered by snow, the thermometers were required to be removed from under the snow and placed on the snow surface in the same way as on the soil surface. Therefore, during the manual observation time period, Ts represents the snow surface temperature when the surface is covered by snow. However, during the 2003–05 period, the meteorological stations were gradually transformed from manual observations to automatic observations (Ying et al. 2006).

In automatic observations, Ts is observed by one platinum metallic thermometer, which is placed in the same way as glass thermometers. However, when there is snow cover on the land surface, an automatic weather station does not require the platinum metallic thermometer to be removed; therefore, Ts measured by an automatic weather station is the temperature of the soil surface beneath the snow cover. Snow cover has a much higher albedo than soil, so that less surface incident solar radiation is absorbed by the surface. At same time, the thermal conductivity of the snow cover is approximately 5–20 times lower than that of soils, which prevents the soil surface under snow from cooling (Zhang 2005). Therefore, the transformation in the method of measurement may lead to a warming effect on the observed record of Ts-min but an unknown effect on the observed record of Ts-max.

b. Reanalysis products

In addition to observations, reanalysis products are one of the most reliable datasets reproducing the state of Ts. The Ts of reanalyses can provide continuous Ts in both clear conditions and cloudy conditions. Compared with the model simulations, the reanalysis system assimilates the observations of surface radiation fluxes and atmospheric conditions from satellites and radiosondes. The assimilation of relevant variables can provide a better boundary condition for numerical forecast models, which is helpful to improve the reliability of the Ts output based on physical consistency. Existing studies have suggested that there is good consistency between the reanalysis and observations, especially for large and long-term scales (Trigo et al. 2015; Tsuang et al. 2008).

Therefore, we chose several recent reanalysis products to explore the variability in Ts and its warming contrast with Ta in China. (Detailed introductions to the reanalysis products are shown in Table 1.) First, we interpolated all the reanalyses on a 1° × 1° grid based on bilinear interpolation and then selected all the grids with valid observed Ts to calculate the national average based on the area-weighted mean. We selected the monthly Ts-mean and Ta-mean from a reanalysis, and the overlap period of all reanalyses was selected as the study period (1980–2014). Then, based on the monthly anomalies, we compare the interannual variability in Ts and the long-term trend of Ts-meanTa-mean between the observations and reanalyses.

Table 1.

Detailed information on the five reanalyses used in this paper. In the assimilation column, 3D-Var (4D-Var) indicates the three-dimensional (four-dimensional) variational assimilation scheme; and 3D-FGAT indicates the three-dimensional variational assimilation scheme with the first guess at the appropriate time.

Detailed information on the five reanalyses used in this paper. In the assimilation column, 3D-Var (4D-Var) indicates the three-dimensional (four-dimensional) variational assimilation scheme; and 3D-FGAT indicates the three-dimensional variational assimilation scheme with the first guess at the appropriate time.
Detailed information on the five reanalyses used in this paper. In the assimilation column, 3D-Var (4D-Var) indicates the three-dimensional (four-dimensional) variational assimilation scheme; and 3D-FGAT indicates the three-dimensional variational assimilation scheme with the first guess at the appropriate time.

c. Model simulations

Although a reanalysis is generally considered to be quite accurate, most of the reanalysis data are available only after 1979 due to the limitations of assimilated data, which are relatively short and unable to provide predictions for future scenarios. Model simulations not only have a long time series, but also simulate the change characteristics of climate factors under different future scenarios (Knutti and Sedláček 2013; Taylor et al. 2012). Therefore, simulations are widely used in the detection, attribution, and prediction of climate change (Rogelj et al. 2012; Screen and Williamson 2017; Wuebbles et al. 2014). The CMIP5 dataset (phase 5 of the Coupled Model Intercomparison Project) constitutes a set of model simulations of the surface and atmospheric state provided by many institutions using different integrated Earth system models and is used widely in climate science (Taylor et al. 2012).

In this paper, we selected 44 simulations from 19 models (see Table 2). They all provide monthly simulations of Ts and Ta during the historical period and four future scenarios (including RCP2.6, RCP4.0, RCP6.0, and RCP8.5) (Taylor et al. 2012). Using an approach similar to that used for the reanalysis datasets, the monthly grid dataset provided by all models was converted into the annual anomalies of the national average. We combined simulations in the historical period with those in different future scenarios to form a long time series of Ts-mean and Ta-mean from 1960 to 2100.

Table 2.

Detailed information on the CMIP5 models used in this paper. The ensembles are the model running in different configuratiosn of physical parameterization schemes. The grid number indicates the number of rows (latitude) and the number of columns (longitude).

Detailed information on the CMIP5 models used in this paper. The ensembles are the model running in different configuratiosn of physical parameterization schemes. The grid number indicates the number of rows (latitude) and the number of columns (longitude).
Detailed information on the CMIP5 models used in this paper. The ensembles are the model running in different configuratiosn of physical parameterization schemes. The grid number indicates the number of rows (latitude) and the number of columns (longitude).

d. Statistical analysis method

The analysis of interannual variabilities and long-term trends were based on the monthly anomalies. The monthly anomalies relative to the 1980–2010 climatology were calculated based on a monthly mean value of the daily values, and when a month was missing more than seven daily values, that month was classified as a missing value. For the annual anomalies, the monthly anomalies were averaged for the entire year. The anomalies in the warm seasons were the averages of the monthly anomalies from May to October, and the anomalies in the cold seasons were the averages of the monthly anomalies from November to April.

The Pearson’s correlation coefficients were calculated by the detrended anomalies that removed the impact of trends to estimate the linear relationship between variables. The long-term trends of variables were evaluated by linear regression. The Student’s t test was performed for the statistical significance of trends, and 95% confidence intervals were provided for trends. Additionally, owing to the uneven distribution of weather stations across China, the national means for all indexes (e.g., correlation coefficients and trends) were calculated on 1° × 1° grids transformed from station observations (see Fig. 1b). Additionally, all trend differences and their confidence intervals between variables were calculated by the time series of differences between variables, rather than the trend and their confidence level of every variable.

Fig. 1.

The large changes in the station observations of land surface temperatures Ts as the manual measurements transform to automatic measurements. (a) Daily values of the observed Ts-mean from 1960 to 2014 from one station (Mohe station; 53°29′N, 122°21′E). (b)–(d) The differences in the multiyear mean between Ts and Ta during the manual observation period (1960–2003). (e)–(g) The differences in the multiyear mean between Ts and Ta during the automatic observation period (2006–14).

Fig. 1.

The large changes in the station observations of land surface temperatures Ts as the manual measurements transform to automatic measurements. (a) Daily values of the observed Ts-mean from 1960 to 2014 from one station (Mohe station; 53°29′N, 122°21′E). (b)–(d) The differences in the multiyear mean between Ts and Ta during the manual observation period (1960–2003). (e)–(g) The differences in the multiyear mean between Ts and Ta during the automatic observation period (2006–14).

3. Detection and attribution of the Ts inhomogeneity

a. Detection of the Ts inhomogeneity

As shown in Fig. 1a, a notable jump occurred in the daily series of Ts when the manual measurement was transformed to an automatic measurement, especially during cold seasons. During the manual observation period, there are no strong spatial patterns in the differences between Ts and Ta (see Figs. 1b–d). In the automatic observation period, the difference between Ts-max and Ta-max is different from that observed during the manual observation period and exhibits no coherent spatial pattern (see Figs. 1b,e). For Ts-min and Ts-mean, the differences between Ts and Ta in the automatic observations have clear spatial patterns and the differences are more positive in North China (see Figs. 1f,g).

The differences between the automatic and manual measurements have substantial effects on the interannual variability in the observed Ts. As shown in Fig. 2a, compared with homogeneous Ta, Ts is consistent in terms of interannual variability and long-term trends from 1960 to 2003. However, Ts-min and Ts-mean sharply increased from 2003 to 2005 and remained stable from 2006 to 2014, which was reported in previous studies (Liu et al. 2017; Wang et al. 2017, 2018). In spatial patterns, the warming difference between Ts-max and Ta-max was largest in Northwest China and decreased to Southeast China (see Fig. 2b), which may be caused by the trend pattern of Rs during this period (Du et al. 2017). The difference in the warming rate between Ts-min (Ts-mean) and Ta-min (Ta-mean) had no significant spatial difference during 1960–2003 (see Figs. 2c,d), but the warming rate of Ts-min (Ts-mean) was much faster than that of Ta-min (Ta-mean) during 2003–14 (see Table S2), especially in North China (Figs. 2f,g). Both patterns of the difference in the multiyear mean and warming rate between Ts and Ta are similar to the spatial pattern for the multiyear mean SD rather than the difference in SD between the automatic and manual periods (see Figs. 1e,f, 2e,f, and Figs. S1c,d).

Fig. 2.

The effect of the automatic transformation of meteorological stations on the interannual variabilities and warming rates of land surface temperatures Ts with near-surface air temperatures Ta as a reference. (a) The yearly anomalies for land surface and near-surface air temperatures from 1960 to 2014. (b)–(d) The differences in trends between land surface and near-surface air temperatures during the manual observation period (1960–2002). (e)–(g) The differences in trends between the land surface and near-surface air temperatures during the automatic observation period (2003–14). (h)–(j) The differences in trends between the land surface and near-surface air temperatures during the entire study period (1960–2014).

Fig. 2.

The effect of the automatic transformation of meteorological stations on the interannual variabilities and warming rates of land surface temperatures Ts with near-surface air temperatures Ta as a reference. (a) The yearly anomalies for land surface and near-surface air temperatures from 1960 to 2014. (b)–(d) The differences in trends between land surface and near-surface air temperatures during the manual observation period (1960–2002). (e)–(g) The differences in trends between the land surface and near-surface air temperatures during the automatic observation period (2003–14). (h)–(j) The differences in trends between the land surface and near-surface air temperatures during the entire study period (1960–2014).

b. Attribution of the Ts inhomogeneity

To prove that the abrupt increase in Ts is caused by transformation to automatic stations rather than climate change, we divided the observed records of Ts during 1965–2014 into five equal length segments (1968–76, 1977–85, 1986–94, 1995–2003, and 2006–14) without the transitional period (2004–05). The first four segments belong to the manual period (1960–2003) and the last segment belongs to the automatic period (2006–14). In addition, all segments are divided into snowy days (SD ≥ 1 cm) and snowless days (SD < 1 cm).

As shown in Figs. 3a–c, for the land–atmosphere temperature difference (TsTa), no significant difference appears among the four segments in the manual observation period, and the difference between the automatic and manual periods is also not obvious on snowless days. However, significant differences emerge between the automatic and manual periods on snowy days. On snowy days, the difference between Ts-min and Ta-min is negative in the manual period and becomes positive in the automatic period (see Fig. 3b), and the difference between Ts-max and Ta-max is smaller in the automatic period than in the manual period (see Fig. 3a). Consequently, the difference between Ts-mean and Ta-mean is much larger in the automatic period than in the manual period (see Fig. 3c). These results indicate that the differences in observed Ts between the automatic and manual periods are mainly caused by the measurement method rather than by climate change or instrument replacement.

Fig. 3.

The abnormal change in TsTa after the observation method changed from manual to automatic and its relationship with SD, for (left) Ts-maxTa-max, (center) Ts-minTa-min, and (right) Ts-meanTa-mean. (a)–(c) The multiyear mean TsTa on snow days and snowless days during the manual observation period (the first four segments) and the automatic observation period (the last segment). (d)–(f) Probability density functions (PDFs) of TsTa on snow days during the five segments mentioned above. (g)–(i) Change characteristics of TsTa with SD. All whiskers denote the standard errors of the corresponding interval.

Fig. 3.

The abnormal change in TsTa after the observation method changed from manual to automatic and its relationship with SD, for (left) Ts-maxTa-max, (center) Ts-minTa-min, and (right) Ts-meanTa-mean. (a)–(c) The multiyear mean TsTa on snow days and snowless days during the manual observation period (the first four segments) and the automatic observation period (the last segment). (d)–(f) Probability density functions (PDFs) of TsTa on snow days during the five segments mentioned above. (g)–(i) Change characteristics of TsTa with SD. All whiskers denote the standard errors of the corresponding interval.

To explore how snow cover affects the Ts observations, we plotted the probability density functions (PDFs) of TsTa on snowy days during the five segments. As shown in Figs. 3d–f, all PDFs during the manual segments are nearly identical and close to the normal distribution. However, there are significant differences in PDFs between the automatic and manual periods. The PDFs during the automatic period are positively skewed and wider than those in the manual period rather than simply moving left or right. The results show that the biases in Ts for the automatic observations are not caused by system errors and are sensitive to changes in other factors, such as SD.

As shown in Fig. 3g, with increasing SD, the value of Ts-maxTa-max sharply decreases and then stabilizes after the SD exceeds 8 cm, which indicates that the impact of snow cover on the observed Ts-max mostly depends on the much higher albedo to the surface solar radiation Rs of snow than soil (Barry 1996). Notably, Ts-min decreases with SD on the snow surface (manual measurement) due to the higher longwave radiation emissivity of snow compared with soil (Dewey 1977; Serreze et al. 1992), whereas Ts-min under a snow layer (automatic measurement) increases with snow depth due to the heat preservation effect of snow cover (see Fig. 3h) (Goodrich 1982; Sokratov and Barry 2002). Therefore, the difference between the automatic and manual measurements is much larger for Ts-min than for Ts-max. As shown in Fig. 3i, Ts-mean observed by automatic measurement is controlled by a higher albedo when the snow layer is thinner and is dominated by thermal insulation for SD exceeding 8 cm. Overall, the heat preservation effect of snow cover on the land surface is much larger than the reflective cooling effect.

Recently, Xu et al. (2019) found a homogeneity problem in the observed Ts, which was caused by the automatic transformation, and the monthly observed Ts was homogenized by the RHtest method. In addition, they removed the records with obvious jumps rather than correcting those errors caused by snow cover, as was done in this study. However, the method used by Xu et al. (2019) has two defects. First, it reduces the spatial coverage of the dataset, resulting in no valid stations in northeast China and northern Xinjiang. Second, as long as there is snow cover, this observation error will exist, although the snowy days for some stations are relatively few and make observation errors not obvious on monthly scales. Importantly, those observation errors will increase with the number of snowy days and cannot be corrected by statistical homogenization methods. Therefore, we used an effective model to correct the observed errors of Ts caused by snow cover during the automatic observation period and then carried out homogenization corrections. The new method can more effectively homogenize the observed records of Ts, and the dataset will have more valid stations and better spatial coverage.

4. Correction of inhomogeneities in the observed Ts

a. Correcting the automatic Ts measurement on snowy days

As shown in Figs. 4a–f, Ts had a significant and stable linear relationship with Ta during the manual period. The national mean correlation coefficients between Ts and Ta are 0.773, 0.859, and 0.821 for the maximum temperatures, mean temperature, and minimum temperature during 1960–2003 (all p < 0.05), respectively. Additionally, we found notable differences in the correlation coefficients between Ts and Ta on snowy days and snowless days.

Fig. 4.

The scatter between Ts and Ta on (a)–(c) snow days and (d)–(f) snowless days during the five segments, respectively (SL is the slope of Ta to Ts; R2 is the determination coefficient of the regression model), showing scatterplots for (left) daily maximum temperatures (Ts-max and Ta-max), (center) daily minimum temperatures (Ts-min and Ta-min), and (right) daily mean temperatures (Ts-mean and Ta-mean). (g) Daily value of Ts after adjusting for abnormal values during 1960–2014 from one station (Mohe station; 53°29′N, 122°21′E). (h) Boxplot of the multiyear mean differences (multiyear mean diff.) between the observed and derived Ts from three methods [method I is Eq. (1); method II is Eq. (2); method III is Eq. (3)] for 1995–2003. (i) The national mean correlation coefficient between the observed and derived Ts from the three methods for 1995–2003. (j) Trend differences (trend diff.) between the observed and derived Ts from the three methods for 1960–2003, and the whiskers denote the 95% confidence intervals.

Fig. 4.

The scatter between Ts and Ta on (a)–(c) snow days and (d)–(f) snowless days during the five segments, respectively (SL is the slope of Ta to Ts; R2 is the determination coefficient of the regression model), showing scatterplots for (left) daily maximum temperatures (Ts-max and Ta-max), (center) daily minimum temperatures (Ts-min and Ta-min), and (right) daily mean temperatures (Ts-mean and Ta-mean). (g) Daily value of Ts after adjusting for abnormal values during 1960–2014 from one station (Mohe station; 53°29′N, 122°21′E). (h) Boxplot of the multiyear mean differences (multiyear mean diff.) between the observed and derived Ts from three methods [method I is Eq. (1); method II is Eq. (2); method III is Eq. (3)] for 1995–2003. (i) The national mean correlation coefficient between the observed and derived Ts from the three methods for 1995–2003. (j) Trend differences (trend diff.) between the observed and derived Ts from the three methods for 1960–2003, and the whiskers denote the 95% confidence intervals.

First, the correlation between Ts-max and Ta-max on snowy days was higher than that on snowless days in the manual observations (see Figs. 4a,d), which is mainly because the surface property (albedo and moisture) of snow cover is more stable than that of soil and thus less affected by flux changes (e.g., evaporation and surface solar radiation). Second, the differences between the correlation coefficients for the daily minimum and maximum temperatures are smaller on snowy days than on snowless days. The correlation coefficient between Ts-min and Ta-min (r = 0.974) is much larger than that between Ts-max and Ta-max (r = 0.829) on snowless days, while there was no significant difference between the correlation coefficients for the daily minimum temperatures (r = 0.958) and maximum temperatures (r = 0.945) on snow days. Overall, the relationship between Ts and Ta is closer and more stable on snow days, especially for the daily maximum temperatures, and the Ta observation are less affected by the automatic transformation (see Fig S2a). Therefore, it is a good choice to use the observed Ta to correct the observation of Ts on snowy days during the automatic observation period.

Figure S4 shows the flowchart of correcting and homogenizing the observed Ts by multiple regression models and the RHtest method with proxy variables. Next, we introduce the process for building regression models between Ts and proxy variables, which provide Ts over the snow surface during the automatic period and the reference series used to homogenize Ts after correction.

Based on the robust linear relationship between Ts and Ta on snowy days, we can derive Ts on snowy days during the automatic observation period by Ta based on a regression model:

 
Ts=a1Ta+c+ε,
(1)

where Ts and Ta denote the daily temperatures of the land surface and near-surface air, respectively; a1 denotes the regression coefficient, c denotes the constant term, and ε is the residual. First, we build the regression model based on the observations of Ts and Ta from 1960–2003. Then, we derive Ts from Ta for 2003–14 based on this model.

According to the analysis in section 3, the SD has a significant impact on the relationship between Ts and Ta. Therefore, we add the SD as a potential predictor for the regression model:

 
Ts=a1Ta+a2SD+c+ε,
(2)

where SD denotes the daily snow depth and a2 is the corresponding regression coefficient for SD. Other variables and coefficients are the same as those in Eq. (1). Based on the analysis in section 4, TsTa has a nonlinear relationship with SD; for example, Ts-maxTa-max would first decrease and then increase slowly with increasing SD (Fig. 3g).

In addition, surface radiation fluxes are essential elements in land–atmosphere interactions. There is a significant difference in the TsTa relationship between clear conditions and cloudy conditions (see Fig. S5). Thus, the changes in the effective cloud cover (aerosol and cloud) have a significant effect on the relationship between Ts and Ta. Previous studies have revealed that the ratio between the surface solar radiation Rs and theoretical maximum surface solar radiation Rc is a metric indicating variabilities in the surface radiation fluxes caused by clouds and aerosols (Feng and Wang 2018). Therefore, based on Eq. (2), we add Rs/Rc as a potential predictor for the regression model:

 
Ts=a1Ta+a2SD+a3(Rs/Rc)+c+ε,
(3)

where Rc and Rs denote the theoretical maximum surface solar radiation on clear days and the surface solar radiation derived from daily observations of sunshine duration, respectively, which were calculated with a hybrid model provided by Yang et al. (2006) and are widely used in existing studies (Du et al. 2017; He et al. 2018; Wang et al. 2015).

Based on the relationship between TsTa and snow depth, we divided the snow depth into four intervals, including SD = 0 cm (snowless), 0 cm < SD ≤ 2 cm, 2 cm < SD ≤ 8 cm, and SD > 8 cm. The coefficients for the regression models are shown in Fig. S6 (SD = 0 cm), Fig. S7 (0 < SD ≤ 2 cm), Fig. S8 (2 < SD ≤ 8 cm), and Fig. S9 (SD > 8 cm). As shown in Figs. S6–S9, the sensitivity of Ts to SD is highest when it is in the interval (0, 2], decreases when SD becomes larger than 2 cm, and tends to stabilize when SD becomes larger than 8 cm. The change coefficients of Ts for SD in different intervals are consistent with the relationship between TsTa and snow depth (see Figs. 3g–i).

Based on the above methods, we built the regression models for each station separately, using observations acquired at each station during the 1960–2002 period. Then, we derived the Ts values on snowy days during 2003–14 and replaced the corresponding observations, and one corrected time series of Ts is shown in Fig. 4g. Compared with Fig. 1a, the anomalous values in the observation series of Ts almost disappeared.

To further verify the accuracy of the derived Ts, we use the observations during 1960–94 to build the regression models and derive Ts for 1995–2003 based on the observations of Ta, SD, and Rs. The comparison between the observed Ts and corrected Ts is used to verify whether the latter can well reflect the mean value, interannual variability, and long-term trend of Ts.

As shown in Fig. 4h, the medians of the differences between the derived and observed Ts based on the three methods are close to 0, and the intervals of the differences decrease when more variables are considered in the regression model. The correlation coefficients between the derived and observed Ts are all greater than 0.91 and notably increase after considering the effects of SD and Rs (see Fig. 4i). Additionally, we analyzed the differences in trends between the observed Ts and derived Ts during 1960–2003. As shown in Fig. 4j, the differences in trends for Ts derived by Eq. (1) are 0.021°C decade−1 (10.9% of the trend of observed Ts), 0.008°C decade−1 (2.3%), and 0.016°C decade−1 (6.7%) for Ts-max, Ts-min, and Ts-mean, respectively. After considering SD and Rs, the differences in trends between the observed and derived Ts decrease to 0.006°C decade−1 (2.1%), 0.004°C decade−1 (0.9%), and 0.005°C decade−1 (1.7%) for Ts-max, Ts-min, and Ts-mean, respectively. The results indicate that the derived Ts can be used to evaluate the long-term trend in land surface temperature.

By comparing the time series before and after correction (see Figs. S9a,b), the anomalous values in the observed Ts nearly disappear after the correction was applied to Ts during the automatic period.

b. Correcting the nonclimatic shifts in the total observed Ts

After correcting the biases caused by the change in observation methods, inhomogeneity issues may remain in some records of Ts, which were caused by instrument replacement and station relocation. Wang et al. (2010) improved an efficient method (RHtest) to homogenize the observed records, which has been implemented in many studies and performs fairly well (L. Cao et al. 2016; Li and Yan 2010). To obtain better results, the records need to match a reference series that is homogeneous and has variabilities similar to those of the target records (X. Wang et al. 2007). However, due to the transformation to automatic stations, there are not enough homogeneous records for use as reference series.

According to the analysis in section 4a, Ts derived from Eq. (3) with homogeneous inputs is a preferable selection as a reference series to homogenize the observed Ts. As shown in Figs. 5a–c, based on the monthly anomalies, the correlation coefficients between the observed Ts and derived Ts at the same station are high across China. By comparing the correlation coefficients between the observed Ts and derived Ts at the same station to those between the observed Ts and that of nearby stations (see Figs. 5d–f), the former is similar to the latter in areas with many stations and is much higher than the latter in areas with few stations such as the Tibetan Plateau. This result indicates that Ts derived by Eq. (3) with homogeneous inputs is a better reference series to detect and adjust the inhomogeneity issues in the observed records of Ts.

Fig. 5.

(a)–(c) Correlation coefficients between the observed Ts and Ts derived from Eq. (3) for 1960–2002. (d)–(f) Maximum correlation coefficients between the time series of observed Ts at a given station with nearby stations for 1960–2014; the station selection method is described in the supplemental material. (g) The annual variability in the significant nonclimatic shifts in daily maximum temperature Ts-max, daily maximum temperature Ts-min, and daily mean temperature Ts-mean. Also shown are the spatial distribution of total shifts in stations for (h) Ts-max, (i) Ts-min, and (j) Ts-mean.

Fig. 5.

(a)–(c) Correlation coefficients between the observed Ts and Ts derived from Eq. (3) for 1960–2002. (d)–(f) Maximum correlation coefficients between the time series of observed Ts at a given station with nearby stations for 1960–2014; the station selection method is described in the supplemental material. (g) The annual variability in the significant nonclimatic shifts in daily maximum temperature Ts-max, daily maximum temperature Ts-min, and daily mean temperature Ts-mean. Also shown are the spatial distribution of total shifts in stations for (h) Ts-max, (i) Ts-min, and (j) Ts-mean.

Based on the penalized maximal t test (X. Wang et al. 2007), 3855, 4716, and 4429 shifts are detected in the corrected records for Ts-max, Ts-min, and Ts-mean, respectively. As shown in Fig. 5g, the number of shifts from 2003 to 2005 is obviously larger than those in other periods, which indicates that inhomogeneities remain in the corrected Ts records around the time of the transformation to automated measurements, despite the correction performed in this study. The spatial distributions of the shifts are displayed in Figs. 5h–j, and most of observed Ts values showed nonclimatic shifts. There are not enough homogeneous records available as references to adjust the inhomogeneous records of Ts; thus, it is better to select the derived Ts from the homogeneous input variables at the same station. Then, the quantile-match (QM) method in the RHtest is used to adjust the significant shifts in the corrected records for Ts, which is the same as the method described in previous studies and introduced in the online supplemental material (Vincent et al. 2002; Xu et al. 2013).

Most inhomogeneities induced by relocations and instrumental changes are relatively small compared to the interannual variability. Therefore, the differences in the results with/without homogenization are not obvious in the time series of daily Ts values (see Fig. S9b). However, these inhomogeneity issues may have a large effect on the estimation of the long-term trend of Ts.

Comparing Fig. 6a with Fig. 2a, we can see that the jump in the time series for Ts-min and Ts-mean from 2003 to 2005 disappeared after corrections. Based on the homogeneous records of Ts, the national mean trends in Ts-max, Ts-min, and Ts-mean were 0.208°, 0.338°, and 0.246°C decade−1, respectively, during 1960–2014 (see Table 3). Similar to Ta, the warming trend of Ts-min is higher than that of Ts-max. Spatially, the warming in Ts follows an obvious latitudinal gradient from a faster rate in Northwest China to a slower rate in Southeast China. However, the warming rate of Ts-max is significantly lower in South China and the North China Plain than in other regions (see Fig. 6b), which is consistent with previous studies and is attributed to the significant decrease in Rs in these regions (Du et al. 2017; Zhou et al. 2017). In contrast, the warming trend of Ts-min in the North China Plain is higher than the trends in nearby regions (see Fig. 6c), which may be caused by increasing aerosol emissions that could reduce the efficiency of longwave cooling during nighttime (C. Cao et al. 2016). Consequently, for Ts-mean, the warming rate is higher in northwest China, northeast China, and the Tibetan Plateau, while the rate is relatively lower in South China, especially in the Sichuan Basin and Yunnan-Guizhou Plateau (see Fig. 6d).

Fig. 6.

(a) Yearly anomalies in Ts after homogenization from 1960–2014. (b)–(d) Warming patterns of Ts from 1960–2014 across China. (e)–(g) The long-term trends in observed Ts in the raw dataset [Ts (Raw)], Ta after homogenization [Ta (Homo)], and Ts after homogenization [Ts (Homo)] during 1960–77, 1978–97, 1998–2014, and 1960–2014. The whiskers denote the 95% confidence intervals for the long-term trends.

Fig. 6.

(a) Yearly anomalies in Ts after homogenization from 1960–2014. (b)–(d) Warming patterns of Ts from 1960–2014 across China. (e)–(g) The long-term trends in observed Ts in the raw dataset [Ts (Raw)], Ta after homogenization [Ta (Homo)], and Ts after homogenization [Ts (Homo)] during 1960–77, 1978–97, 1998–2014, and 1960–2014. The whiskers denote the 95% confidence intervals for the long-term trends.

Table 3.

The warming trends of land surface temperatures based on the original observations (Raw) and observations after corrections (Homo), and the warming trend of homogenized near-surface air temperature observations during 1960–2014 over China. Note that the ±95% confidence interval and all of the trend are significant at 95% confidence level.

The warming trends of land surface temperatures based on the original observations (Raw) and observations after corrections (Homo), and the warming trend of homogenized near-surface air temperature observations during 1960–2014 over China. Note that the ±95% confidence interval and all of the trend are significant at 95% confidence level.
The warming trends of land surface temperatures based on the original observations (Raw) and observations after corrections (Homo), and the warming trend of homogenized near-surface air temperature observations during 1960–2014 over China. Note that the ±95% confidence interval and all of the trend are significant at 95% confidence level.

5. Warming contrast between Ts and Ta

a. Warming contrast between Ts and Ta in observations

From 1960 to 2014, the warming of both Ts and Ta showed significant interdecadal changes, which can be roughly divided into three stages (see Figs. 2a and 6a). During the first stage (1960–77), both Ts and Ta were in a cooling trend, especially for maximum temperature, and the main reason for the cooling trend is attributed to the global dimming that occurred in this period (He et al. 2018). However, the cooling rate of Ts-max is faster than that of Ta-max in this period. Site migrations and instrument alterations were frequent during this period and could cause serious inhomogeneity issues in the observations (see Fig. 5g), which may be the reason for the trend difference between Ts-max and Ta-max; hence, the trend difference was greatly reduced after corrections (see Fig. 6e).

In the second stage (1978–98), both Ts and Ta showed a significant warming trend with an increasing rate. The main reason is that the global temperature entered a rapid warming stage with global diming gradually disappearing, and faster urbanization resulted in a warming effect on the station observations over China during this period (Sun et al. 2016; Wild 2009). However, there was no significant difference in the warming rate between Ts and Ta during 1978–98 (see Figs. 6e–g).

The third stage occurred after 1998, where the global warming slowed to a hiatus, and both Ts and Ta over China stopped rising, reaching stagnation or even decline. In this stage, the trend difference between Ts-min and Ta-min was mainly caused by the transformation from manual to automatic observations, and this trend difference almost disappeared after corrections (see Fig. 6f).

Spatially, by comparing Figs. 7b–d and 2h–j, the abnormal trend difference between Ts and Ta in northern China almost disappeared after corrections. The trend difference between Ts and Ta generally disappears and there is no obvious seasonal and diurnal contrast in the trend difference between Ts and Ta over China during 1960–2014 (see Fig. 7g).

Fig. 7.

(a) Yearly anomalies of the difference between Ts and Ta in the original observations (Raw) and the observations after homogenization (Homo) from 1960–2014. (b)–(d) Warming patterns of the difference between Ts and Ta after homogenization across China during 1960–2014. (e)–(g) The long-term trends of the difference between Ts and Ta, showing annual, warm season (warm), and cold season (cold) values. The whiskers denote the 95% confidence intervals for the long-term trends.

Fig. 7.

(a) Yearly anomalies of the difference between Ts and Ta in the original observations (Raw) and the observations after homogenization (Homo) from 1960–2014. (b)–(d) Warming patterns of the difference between Ts and Ta after homogenization across China during 1960–2014. (e)–(g) The long-term trends of the difference between Ts and Ta, showing annual, warm season (warm), and cold season (cold) values. The whiskers denote the 95% confidence intervals for the long-term trends.

Overall, based on the dataset after corrections, there is no significant trend difference between Ts and Ta during 1960–2014 (see Table S2). Existing studies have shown that the warming rate of Ta after 2000 is significantly higher than that of Ts, which is mainly caused by the inhomogeneity in the observed Ts that results from the automatic transformation (see Fig. 7a).

b. Warming contrast between Ts and Ta in the reanalysis

As shown in Fig. 8, after corrections, the observed Ts and Ts from the reanalysis products have good consistency in terms of interannual variability. Compared with the reanalysis products before 1979 (only NCEP-R2 and JRA-55), the consistency in the observation and reanalysis products after 1979 (in the satellite age) is better. This finding suggests that the assimilation of satellite products can notably improve the ability of the reanalysis system to reproduce Ts-mean. By comparing the raw observations of Ts-mean with the reanalysis, the raw Ts-mean was consistent with the reanalyses before 2003 but became significantly higher than the reanalyses after 2003, especially in the cold seasons. This phenomenon is consistent with the change in observation practices between the manual and automated Ts measurements explored in this paper. However, after corrections, the observed Ts-mean agreed well with the reanalyses, which indicated the rationality and effectiveness of our homogenization method (see Fig. 8). In terms of the reanalysis products, compared with the observed Ts-mean, the simulation capacity of Ts-mean in ERA5 is significantly higher than that of other reanalysis products, followed by JRA-55, ERA-Interim, and CFSR, while MERRA2 and NCEP-R2 are relatively poor.

Fig. 8.

(left) Yearly anomalies of Ts-mean in five reanalyses [CFSR, ERA-Interim (ERA-I), ERA5, JRA-55, MERRA2, and NCEP-R2] and in situ observations without homogenization (Obs-Raw) and those after homogenization (Obs-Homo) showing (a) annual, (c) warm season, and (e) cold season values in China. (right) Taylor diagrams for the annual anomalies of the raw observed Ts and five reanalyzed Ts from 1980 to 2014 in China: (b) annual, (d) warm seasons, and (f) cold seasons. The correlation coefficient, standard deviation, and the root-mean-square deviation (RMSD) are calculated against the yearly anomaly of the observed Ts after homogenization.

Fig. 8.

(left) Yearly anomalies of Ts-mean in five reanalyses [CFSR, ERA-Interim (ERA-I), ERA5, JRA-55, MERRA2, and NCEP-R2] and in situ observations without homogenization (Obs-Raw) and those after homogenization (Obs-Homo) showing (a) annual, (c) warm season, and (e) cold season values in China. (right) Taylor diagrams for the annual anomalies of the raw observed Ts and five reanalyzed Ts from 1980 to 2014 in China: (b) annual, (d) warm seasons, and (f) cold seasons. The correlation coefficient, standard deviation, and the root-mean-square deviation (RMSD) are calculated against the yearly anomaly of the observed Ts after homogenization.

For the long-term trend, the long-term trend of Ta-mean in the reanalysis is generally consistent with the observations (see Fig. 9). However, the raw Ts-mean is significantly higher than those in the reanalyses, especially in the cold seasons (see Fig. 9a). After corrections, this trend difference between observations and reanalyses essentially disappears. For all reanalyses, the differences between Ts-mean and Ta-mean have no significant trend at the 95% confidence level during 1980–2014, except ERA-Interim, which shows a slight downward trend (see Fig. 9c). The temperature difference between Ts-mean and Ta-mean in the raw observations increased significantly, which is consistent with previous studies but disagrees with those in the reanalyses. There is no significant warming contrast between Ts-mean and Ta-mean based on the observation after corrections, which agrees well with the reanalyses (see Fig. 9c).

Fig. 9.

The long-term trend of (a) Ts-mean and (b) Ta-mean, and (c) the difference between Ts-mean and Ta-mean in the five reanalyses [CFSR, ERA-Interim (ERA-I), ERA5, JRA-55, MERRA2, and NCEP-R2] and in situ observations before and after homogenization in China during 1980–2014. The whiskers denote the 95% confidence intervals for the long-term trends.

Fig. 9.

The long-term trend of (a) Ts-mean and (b) Ta-mean, and (c) the difference between Ts-mean and Ta-mean in the five reanalyses [CFSR, ERA-Interim (ERA-I), ERA5, JRA-55, MERRA2, and NCEP-R2] and in situ observations before and after homogenization in China during 1980–2014. The whiskers denote the 95% confidence intervals for the long-term trends.

Note that Ts in the reanalyses is defined as the temperature of the surface at the radiative energy balance, which refers to the interface between the surface of the soil, snow or ice, and the atmosphere. However, the influence of vegetation on the radiative equilibrium of the surface is considered in the Ts calculation in the reanalysis, which is different from the observed Ts over the vegetation area. Based on the above comparison, the influence of vegetation on Ts may be obvious in the multiyear mean or variabilities at small spatiotemporal scales, but not significant in the interannual variability and long-term trend (see Figs. 8 and 9).

c. Warming contrast between Ts and Ta in the model simulations

As shown in Fig. 10, from 1960 to 2014, CMIP5 models significantly underestimated the interannual variabilities in Ts-mean and Ta-mean, but captured their interdecadal variabilities well (see Figs. 10a,c). For the long-term trend, all models underestimated the warming trend of Ts-mean and Ta-mean in China during 1960–2014 (see Figs. 10b,d). The difference between the simulated Ts-mean and Ta-mean in all models shows a downward trend (see Fig. 10f), which is different from the results of this study. However, the trend difference between Ts-mean and Ta-mean is only from approximately −1.36% to −3.19% (from −1.34% to −3.09%; the percentages are negative because the Ts-meanTa-mean trend is negative; see Fig. 10f) of the trend of Ts-mean (Ta-mean) and generally negligible, which is consistent with the results derived from the homogeneous observations and reanalysis productions.

Fig. 10.

Yearly anomalies of (a) Ts-mean and (c) Ta-mean, and (e) the difference between Ts-mean and Ta-mean in CMIP5 models (see Table 2) and in situ observations after corrections during 1960–2100 in China. The time series from the CMIP5 models consists of the historical periods (1960–2005) and four future scenarios (RCP2.6, RCP4.5, RCP6.0, and RCP8.5 during 2006–2100). Also shown are long-term trends of (b) Ts-mean and (d) Ta-mean, and (f) the difference between Ts-mean and Ta-mean in future scenarios and in situ observations after homogenization in China during 1960–2014 and 2015–2100. The solid lines represent the multimodel mean and observations, and the shading indicates the spread of the ensemble mean.

Fig. 10.

Yearly anomalies of (a) Ts-mean and (c) Ta-mean, and (e) the difference between Ts-mean and Ta-mean in CMIP5 models (see Table 2) and in situ observations after corrections during 1960–2100 in China. The time series from the CMIP5 models consists of the historical periods (1960–2005) and four future scenarios (RCP2.6, RCP4.5, RCP6.0, and RCP8.5 during 2006–2100). Also shown are long-term trends of (b) Ts-mean and (d) Ta-mean, and (f) the difference between Ts-mean and Ta-mean in future scenarios and in situ observations after homogenization in China during 1960–2014 and 2015–2100. The solid lines represent the multimodel mean and observations, and the shading indicates the spread of the ensemble mean.

During 2015–2100, the interannual variabilities and long-term trends of Ts-mean and Ta-mean showed no significant difference in each scenario (Figs. 10b,d). Therefore, the difference between Ts and Ta shows a slight upward trend by 0.005°–0.006°C decade−1, accounting for 1.10%–10.14% (1.11%–11.28%) of the trend for Ts-mean (Ta-mean), and there is no significant difference in the warming contrast between Ts and Ta in different scenarios (see Fig. 10f). Therefore, under different global warming scenarios, the trend difference between Ts-mean and Ta-mean is very weak and much smaller than the trend of Ts-mean or Ta-mean. Therefore, model simulations from CMIP5 indicate that the global warming resulting from increasing greenhouse gas concentrations will have no significant effect on the difference between Ts-mean and Ta-mean over China.

The different responses between Ts and Ta to global change may result in a systemic bias in the evaluation of climate simulations. Previously, Cowtan et al. (2015) compared the simple GMST (Ta over both land and ocean) and the blended GMST (Ta over land and Ts over ocean) and found that the warming rate of the former is faster than that of the latter, indicating that the warming rate of Ts is substantially slower than that of Ta over the ocean. However, their study did not compare the surface temperature and near-surface temperature over the ocean. Our study provides a good complement to this research. We did compare Ts and Ta over mainland China in our study, but more work is needed to explore whether there is a different response between Ts and Ta to climate change over all global landmasses.

6. Conclusions

As a primary metric indicating changes in the land surface energy balance, the land surface temperature plays an important role in the response and feedback of land–atmosphere interactions with respect to climate change. Recently, a high-quality dataset of historical records for Ts across China has been developed by the CMA, which makes it possible to explore the variability in Ts under climate change. However, the raw records of Ts have significant inhomogeneity issues that were ignored in most previous studies, which resulted in some abnormalities in the long-term variability in Ts. The most significant abnormality is a jump in warming for Ts-min and Ts-mean during the cold seasons after 2003, while the near-surface air temperature exhibits a “warming hiatus.”

By analysis, we find that the abnormal jump in observed Ts is caused by a change in the measurement method: the land surface temperature was measured on the snow surface though manual observations until 2003 and then transformed gradually to be measured at the snow bottom by automatic observations. Additionally, we find a significant and stable linear relationship between Ts and Ta, especially on snow days, and this relationship is sensitive to changes in SD and surface solar radiation. Therefore, we develop a linear regression model between Ta, SD, Rs, and Ts based on observations during 1960–2003, and then derive Ts on snowy days during 2003–14. By comparing the observed and derived Ts, evaluating the long-term variability in Ts based on records with derived values is reliable. Additionally, by using the derived Ts from Eq. (3) with homogeneous inputs as a reference, we homogenized the observed Ts using RHtest.

After adjustment for the interannual variability, Ts rapidly warmed until 1998 and then slowed, which is consistent with the warming hiatus in Ta during this period. The warming trends for Ts-max, Ts-min, and Ts-mean are 0.208°, 0.338°, and 0.246°C decade−1, respectively, during 1960–2014 without an abnormal jump in Ts in 2003. Spatially, the warming trend in Ts is fastest in North China and decreases from northwest to southeast. Additionally, northeastern China and the North China Plain show slower warming for Ts-max but faster warming for Ts-min, which may be caused by the heavy air pollution in these regions. Note that the warming trends in both Ts and Ta slow after 1998, which may imply that land–atmosphere interaction plays an important role in the warming hiatus during this period.

Based on homogeneous observations, reanalyses, and model simulations, we analyzed the warming contrast between Ts and Ta in China under the background of global warming. The results show that there was no significant difference between Ts and Ta in the warming trend during 1960–2014 for observations and reanalyses. Although the CMIP5 model simulations showed a slight downward trend from 1960 to 2014 and a slight upward trend from 2015 to 2100 in different scenarios, the magnitude of all trends was much smaller than that of Ts or Ta. These results further proved that the large warming contrast between Ts and Ta over China found in previous studies is mainly caused by the inhomogeneity in observed Ts caused by automatic transformation.

Acknowledgments

This study was funded by the National Basic Research Program of China (2017YFA0603601), Postdoctoral Innovative Talent Program (212400247), Postdoctoral Science Foundation of China (212400249), and the National Natural Science Foundation of China (41525018). The latest observational data were obtained from the CMA (http://www.cma.gov.cn).

REFERENCES

REFERENCES
Anderson
,
M. C.
,
J. M.
Norman
,
W. P.
Kustas
,
R.
Houborg
,
P. J.
Starks
, and
N.
Agam
,
2008
:
A thermal-based remote sensing technique for routine mapping of land-surface carbon, water and energy fluxes from field to regional scales
.
Remote Sens. Environ.
,
112
,
4227
4241
, https://doi.org/10.1016/j.rse.2008.07.009.
Barry
,
R. G.
,
1996
:
The parameterization of surface albedo for sea ice and its snow cover
.
Prog. Phys. Geogr.
,
20
,
63
79
, https://doi.org/10.1177/030913339602000104.
Bosilovich
,
M. G.
,
2006
:
A comparison of MODIS land surface temperature with in situ observations
.
Geophys. Res. Lett.
,
33
,
L20112
, https://doi.org/10.1029/2006GL027519.
Cao
,
C.
,
X.
Lee
,
S.
Liu
,
N.
Schultz
,
W.
Xiao
,
M.
Zhang
, and
L.
Zhao
,
2016
:
Urban heat islands in China enhanced by haze pollution
.
Nat. Commun.
,
7
,
12509
, https://doi.org/10.1038/ncomms12509.
Cao
,
L.
,
Y.
Zhu
,
G.
Tang
,
F.
Yuan
, and
Z.
Yan
,
2016
:
Climatic warming in China according to a homogenized data set from 2419 stations
.
Int. J. Climatol.
,
36
,
4384
4392
, https://doi.org/10.1002/joc.4639.
Catherinot
,
J.
,
C.
Prigent
,
R.
Maurer
,
F.
Papa
,
C.
Jimenez
,
F.
Aires
, and
W.
Rossow
,
2011
:
Evaluation of “all weather” microwave-derived land surface temperatures with in situ CEOP measurements
.
J. Geophys. Res.
,
116
,
D23105
, https://doi.org/10.1029/2011JD016439.
Cowtan
,
K.
, and et al
,
2015
:
Robust comparison of climate models with observations using blended land air and ocean sea surface temperatures
.
Geophys. Res. Lett.
,
42
,
6526
6534
, https://doi.org/10.1002/2015GL064888.
Cramer
,
W.
, and et al
,
2013
: Detection and attribution of observed impacts. Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects, C. B. Field et al., Eds., Cambridge University Press, 979–1037, https://www.ipcc.ch/site/assets/uploads/2018/02/WGIIAR5-Chap18_FINAL.pdf.
Dewey
,
K. F.
,
1977
:
Daily maximum and minimum temperature forecasts and the influence of snow cover
.
Mon. Wea. Rev.
,
105
,
1594
1597
, https://doi.org/10.1175/1520-0493(1977)105<1594:DMAMTF>2.0.CO;2.
Du
,
J.
,
K.
Wang
,
J.
Wang
, and
Q.
Ma
,
2017
:
Contributions of surface solar radiation and precipitation to the spatiotemporal patterns of surface and air warming in China from 1960 to 2003
.
Atmos. Chem. Phys.
,
17
,
4931
4944
, https://doi.org/10.5194/acp-17-4931-2017.
Feng
,
F.
, and
K.
Wang
,
2018
:
Does the Modern-Era Retrospective Analysis for Research and Applications-2 aerosol reanalysis introduce an improvement in the simulation of surface solar radiation over China?
Int. J. Climatol.
,
39
,
1305
1318
, https://doi.org/10.1002/joc.5881.
Fuchs
,
M.
, and
C.
Tanner
,
1968
:
Surface temperature measurements of bare soils
.
J. Appl. Meteor.
,
7
,
303
305
, https://doi.org/10.1175/1520-0450(1968)007<0303:STMOBS>2.0.CO;2.
Good
,
E. J.
,
D. J.
Ghent
,
C. E.
Bulgin
, and
J. J.
Remedios
,
2017
:
A spatiotemporal analysis of the relationship between near-surface air temperature and satellite land surface temperatures using 17 years of data from the ATSR series
.
J. Geophys. Res. Atmos.
,
122
,
9185
9210
, https://doi.org/10.1002/2017JD026880.
Goodrich
,
L.
,
1982
:
The influence of snow cover on the ground thermal regime
.
Can. Geotech. J.
,
19
,
421
432
, https://doi.org/10.1139/t82-047.
Hansen
,
J.
,
R.
Ruedy
,
M.
Sato
, and
K.
Lo
,
2010
:
Global surface temperature change
.
Rev. Geophys.
,
48
,
RG4004
, https://doi.org/10.1029/2010RG000345.
He
,
Y.
,
K.
Wang
,
C.
Zhou
, and
M.
Wild
,
2018
:
A revisit of global dimming and brightening based on the sunshine duration
.
Geophys. Res. Lett.
,
45
,
4281
4289
, https://doi.org/10.1029/2018GL077424.
Jiménez
,
C.
,
C.
Prigent
,
J.
Catherinot
,
W.
Rossow
,
P.
Liang
, and
J. L.
Moncet
,
2012
:
A comparison of ISCCP land surface temperature with other satellite and in situ observations
.
J. Geophys. Res.
,
117
,
D08111
, https://doi.org/10.1029/2011JD017058.
Jin
,
M. S.
,
2012
:
Developing an index to measure urban heat island effect using satellite land skin temperature and land cover observations
.
J. Climate
,
25
,
6193
6201
, https://doi.org/10.1175/JCLI-D-11-00509.1.
Jin
,
M. S.
,
R.
Dickinson
, and
A.
Vogelmann
,
1997
:
A comparison of CCM2–BATS skin temperature and surface-air temperature with satellite and surface observations
.
J. Climate
,
10
,
1505
1524
, https://doi.org/10.1175/1520-0442(1997)010<1505:ACOCBS>2.0.CO;2.
Jones
,
P. D.
,
M.
New
,
D. E.
Parker
,
S.
Martin
, and
I. G.
Rigor
,
1999
:
Surface air temperature and its changes over the past 150 years
.
Rev. Geophys.
,
37
,
173
199
, https://doi.org/10.1029/1999RG900002.
Karnieli
,
A.
,
N.
Agam
,
R. T.
Pinker
,
M.
Anderson
,
M. L.
Imhoff
,
G. G.
Gutman
,
N.
Panov
, and
A.
Goldberg
,
2010
:
Use of NDVI and land surface temperature for drought assessment: Merits and limitations
.
J. Climate
,
23
,
618
633
, https://doi.org/10.1175/2009JCLI2900.1.
Knutti
,
R.
, and
J.
Sedláček
,
2013
:
Robustness and uncertainties in the new CMIP5 climate model projections
.
Nat. Climate Change
,
3
,
369
373
, https://doi.org/10.1038/nclimate1716.
Koch
,
J.
,
A.
Siemann
,
S.
Stisen
, and
J.
Sheffield
,
2016
:
Spatial validation of large-scale land surface models against monthly land surface temperature patterns using innovative performance metrics
.
J. Geophys. Res. Atmos.
,
121
,
5430
5452
, https://doi.org/10.1002/2015JD024482.
Koike
,
T.
,
2004
:
The coordinated enhanced observing period-an initial step for integrated global water cycle observation
.
WMO Bull.
,
53
,
115
121
.
Li
,
Z.
, and
Z.
Yan
,
2010
:
Application of multiple analysis of series for homogenization to Beijing daily temperature series (1960–2006)
.
Adv. Atmos. Sci.
,
27
,
777
787
, https://doi.org/10.1007/s00376-009-9052-0.
Liu
,
B.
,
M.
Henderson
,
L.
Wang
,
X.
Shen
,
D.
Zhou
, and
X.
Chen
,
2017
:
Climatology and trends of air and soil surface temperatures in the temperate steppe region of North China
.
Int. J. Climatol.
,
37
,
1199
1209
, https://doi.org/10.1002/joc.5076.
Liu
,
Y.
,
T.
Hiyama
, and
Y.
Yamaguchi
,
2006
:
Scaling of land surface temperature using satellite data: A case examination on ASTER and MODIS products over a heterogeneous terrain area
.
Remote Sens. Environ.
,
105
,
115
128
, https://doi.org/10.1016/j.rse.2006.06.012.
Mann
,
M. E.
,
S. K.
Miller
,
S.
Rahmstorf
,
B. A.
Steinman
, and
M.
Tingley
,
2017
:
Record temperature streak bears anthropogenic fingerprint
.
Geophys. Res. Lett.
,
44
,
7936
7944
, https://doi.org/10.1002/2017GL074056.
Mühlbauer
,
S.
,
A. C.
Costa
, and
M.
Caetano
,
2016
:
A spatiotemporal analysis of droughts and the influence of North Atlantic Oscillation in the Iberian Peninsula based on MODIS imagery
.
Theor. Appl. Climatol.
,
124
,
703
721
, https://doi.org/10.1007/s00704-015-1451-9.
Rogelj
,
J.
,
M.
Meinshausen
, and
R.
Knutti
,
2012
:
Global warming under old and new scenarios using IPCC climate sensitivity range estimates
.
Nat. Climate Change
,
2
,
248
253
, https://doi.org/10.1038/nclimate1385.
Schwarz
,
N.
,
S.
Lautenbach
, and
R.
Seppelt
,
2011
:
Exploring indicators for quantifying surface urban heat islands of European cities with MODIS land surface temperatures
.
Remote Sens. Environ.
,
115
,
3175
3186
, https://doi.org/10.1016/j.rse.2011.07.003.
Screen
,
J.
, and
D.
Williamson
,
2017
:
Ice-free Arctic at 1.5°C?
.
Nat. Climate Change
,
7
,
230
231
, https://doi.org/10.1038/nclimate3248.
Serreze
,
M. C.
,
R. C.
Schnell
, and
J. D.
Kahl
,
1992
:
Low-level temperature inversions of the Eurasian Arctic and comparisons with Soviet drifting station data
.
J. Climate
,
5
,
615
629
, https://doi.org/10.1175/1520-0442(1992)005<0615:LLTIOT>2.0.CO;2.
Singh
,
R.
,
C.
Singh
,
S. P.
Ojha
,
A. S.
Kumar
,
C.
Kishtawal
, and
A.
Kumar
,
2016
:
Land surface temperature from INSAT-3D imager data: Retrieval and assimilation in NWP model
.
J. Geophys. Res. Atmos.
,
121
,
6909
6926
, https://doi.org/10.1002/2016JD024752.
Sokratov
,
S. A.
, and
R. G.
Barry
,
2002
:
Intraseasonal variation in the thermoinsulation effect of snow cover on soil temperatures and energy balance
.
J. Geophys. Res.
,
107
,
4093
, https://doi.org/10.1029/2001JD000489.
Sun
,
Y.
,
X. B.
Zhang
,
G. Y.
Ren
,
F. W.
Zwiers
, and
T.
Hu
,
2016
:
Contribution of urbanization to warming in China
.
Nat. Climate Change
,
6
,
706
709
, https://doi.org/10.1038/nclimate2956.
Taylor
,
K. E.
,
R. J.
Stouffer
, and
G. A.
Meehl
,
2012
:
An overview of CMIP5 and the experiment design
.
Bull. Amer. Meteor. Soc.
,
93
,
485
498
, https://doi.org/10.1175/BAMS-D-11-00094.1.
Trigo
,
I.
,
S.
Boussetta
,
P.
Viterbo
,
G.
Balsamo
,
A.
Beljaars
, and
I.
Sandu
,
2015
:
Comparison of model land skin temperature with remotely sensed estimates and assessment of surface–atmosphere coupling
.
J. Geophys. Res. Atmos.
,
120
,
12 096
12 111
, https://doi.org/10.1002/2015JD023812.
Tsuang
,
B.-J.
,
M.-D.
Chou
,
Y.
Zhang
,
A.
Roesch
, and
K.
Yang
,
2008
:
Evaluations of land–ocean skin temperatures of the ISCCP satellite retrievals and the NCEP and ERA reanalyses
.
J. Climate
,
21
,
308
330
, https://doi.org/10.1175/2007JCLI1502.1.
Vincent
,
L. A.
,
X.
Zhang
,
B. R.
Bonsal
, and
W. D.
Hogg
,
2002
:
Homogenization of daily temperatures over Canada
.
J. Climate
,
15
,
1322
1334
, https://doi.org/10.1175/1520-0442(2002)015<1322:HODTOC>2.0.CO;2.
Wang
,
A.
,
M.
Barlage
,
X.
Zeng
, and
C. S.
Draper
,
2014
:
Comparison of land skin temperature from a land model, remote sensing, and in situ measurement
.
J. Geophys. Res. Atmos.
,
119
,
3093
3106
, https://doi.org/10.1002/2013JD021026.
Wang
,
K.
, and
R. E.
Dickinson
,
2013
:
Contribution of solar radiation to decadal temperature variability over land
.
Proc. Natl. Acad. Sci. USA
,
110
,
14 877
14 882
, https://doi.org/10.1073/pnas.1311433110.
Wang
,
K.
,
Z.
Wan
,
P.
Wang
,
M.
Sparrow
,
J.
Liu
,
X.
Zhou
, and
S.
Haginoya
,
2005
:
Estimation of surface long wave radiation and broadband emissivity using Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature/emissivity products
.
J. Geophys. Res.
,
110
,
D11109
, https://doi.org/10.1029/2004JD005566.
Wang
,
K.
,
J.
Wang
,
P.
Wang
,
M.
Sparrow
,
J.
Yang
, and
H.
Chen
,
2007
:
Influences of urbanization on surface characteristics as derived from the moderate-resolution imaging spectroradiometer: A case study for the Beijing metropolitan area
.
J. Geophys. Res.
,
112
,
D22S06
, https://doi.org/10.1029/2006JD007997.
Wang
,
K.
,
Q.
Ma
,
Z. J.
Li
, and
J. K.
Wang
,
2015
:
Decadal variability of surface incident solar radiation over China: Observations, satellite retrievals, and reanalyses
.
J. Geophys. Res. Atmos.
,
120
,
6500
6514
, https://doi.org/10.1002/2015JD023420.
Wang
,
L.
,
M.
Henderson
,
B.
Liu
,
X.
Shen
,
X.
Chen
,
L.
Lian
, and
D.
Zhou
,
2018
:
Maximum and minimum soil surface temperature trends over China, 1965–2014
.
J. Geophys. Res. Atmos.
,
123
,
2004
2016
, https://doi.org/10.1002/2017JD027283.
Wang
,
X.
, and
Y.
Feng
,
2013
: RHtests V4 User Manual. Climate Research Division, Atmospheric Science and Technology Directorate Science and Technology Branch, Environment Canada, 29 pp., http://etccdi.pacificclimate.org/RHtest/RHtestsV4_UserManual_20July2013.pdf.
Wang
,
X.
,
Q.
Wen
, and
Y.
Wu
,
2007
:
Penalized maximal t test for detecting undocumented mean change in climate data series
.
J. Appl. Meteor. Climatol.
,
46
,
916
931
, https://doi.org/10.1175/JAM2504.1.
Wang
,
X.
,
H.
Chen
,
Y.
Wu
,
Y.
Feng
, and
Q.
Pu
,
2010
:
New techniques for the detection and adjustment of shifts in daily precipitation data series
.
J. Appl. Meteor. Climatol.
,
49
,
2416
2436
, https://doi.org/10.1175/2010JAMC2376.1.
Wang
,
Y.
,
Z. Z.
Hu
, and
F.
Yan
,
2017
:
Spatiotemporal variations of differences between surface air and ground temperatures in China
.
J. Geophys. Res. Atmos.
,
122
,
7990
7999
, https://doi.org/10.1002/2016JD026110.
Wild
,
M.
,
2009
:
Global dimming and brightening: A review
.
J. Geophys. Res.
,
114
,
D00D16
, https://doi.org/10.1029/2008JD011470.
Wuebbles
,
D.
, and et al
,
2014
:
CMIP5 climate model analyses: Climate extremes in the United States
.
Bull. Amer. Meteor. Soc.
,
95
,
571
583
, https://doi.org/10.1175/BAMS-D-12-00172.1.
Xu
,
W.
,
Q.
Li
,
X. L.
Wang
,
S.
Yang
,
L.
Cao
, and
Y.
Feng
,
2013
:
Homogenization of Chinese daily surface air temperatures and analysis of trends in the extreme temperature indices
.
J. Geophys. Res. Atmos.
,
118
,
9708
9720
, https://doi.org/10.1002/JGRD.50791.
Xu
,
W.
,
C.
Sun
,
J.
Zuo
,
Z.
Ma
,
W.
Li
, and
S.
Yang
,
2019
:
Homogenization of monthly ground surface temperature in China during 1961–2016 and performances of GLDAS reanalysis products
.
J. Climate
,
32
,
1121
1135
, https://doi.org/10.1175/JCLI-D-18-0275.1.
Yang
,
K.
,
T.
Koike
, and
B.
Ye
,
2006
:
Improving estimation of hourly, daily, and monthly solar radiation by importing global data sets
.
Agric. For. Meteor.
,
137
,
43
55
, https://doi.org/10.1016/j.agrformet.2006.02.001.
Ying
,
W.
,
L.
Xiaoning
, and
J.
Xiaohui
,
2006
: Differences between automatic and manual meteorological observation. WMO Tech. Conf. on Instruments and Methods of Observations, Geneva, Switzerland, WMO, 4–6.
Zhang
,
H.
,
E.
Wang
,
D.
Zhou
,
Z.
Luo
, and
Z.
Zhang
,
2016
:
Rising soil temperature in China and its potential ecological impact
.
Sci. Rep.
,
6
,
35530
, https://doi.org/10.1038/srep35530.
Zhang
,
T.
,
2005
:
Influence of the seasonal snow cover on the ground thermal regime: An overview
.
Rev. Geophys.
,
43
,
RG4002
, https://doi.org/10.1029/2004RG000157.
Zhou
,
C.
,
K.
Wang
, and
Q.
Ma
,
2017
:
Evaluation of eight current reanalyses in simulating land surface temperature from 1979 to 2003 in China
.
J. Climate
,
30
,
7379
7398
, https://doi.org/10.1175/JCLI-D-16-0903.1.
Zhou
,
C.
,
Y.
He
, and
K.
Wang
,
2018
:
On the suitability of current atmospheric reanalyses for regional warming studies over China
.
Atmos. Chem. Phys.
,
18
,
8113
8136
, https://doi.org/10.5194/acp-18-8113-2018.
For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).