Contributions of surface solar radiation and precipitation to the spatiotemporal patterns of surface and air warming in China from 1960 to 2003

Although global warming has been attributed to increases in atmospheric greenhouses gases, the mechanisms underlying spatiotemporal patterns of warming trends remain under debate. Herein, we analyzed surface and air warming observations recorded at 1977 stations in China from 1960 to 2003. Our results showed a significant spatial pattern for the warming of the daily maximum surface (Ts-max) and air (Ta-max) temperatures, and the pattern was stronger in northwest and northeast China and weaker or negative in South China and the North China Plain. These warming spatial patterns were attributed to surface shortwave solar radiation (Rs) and precipitation (P), which play a key role in the surface energy budget. During the study period, Rs decreased by −1.50± 0.42 W m−2 10 yr−1 in China, which reduced the trends of Ts-max and Ta-max by about 0.139 and 0.053 C 10 yr−1, respectively. More importantly, the decreasing rates in South China and the North China Plain were stronger than those in other parts of China. The spatial contrasts in the trends of Ts-max and Ta-max in China were significantly reduced after adjusting for the effect of Rsand P . For example, after adjusting for the effect of Rs and P , the difference in the Ts-max and Ta-max values between the North China Plain and the Loess Plateau was reduced by 97.8 and 68.3 %, respectively; the seasonal contrast in Ts-max and Ta-max decreased by 45.0 and 17.2 %, respectively; and the daily contrast in the warming rates of the surface and air temperature decreased by 33.0 and 29.1 %, respectively. This study shows that the land energy budget plays an essential role in the identification of regional warming patterns.


Introduction
Increases in observational data and rapid developments in simulation capacity of climate models have provided evidence for the phenomenon of global warming (Hartmann et al., 2013), and the increases in anthropogenic greenhouse gases and other anthropogenic effects are considered as the primary causes.However, significant spatial and temporal heterogeneities in climate warming have been observed.For example, faster warming rates occur in semiarid regions and a "warming hole" has been identified in the central United States (Boyles and Raman, 2003;Huang et al., 2012).These spatiotemporal heterogeneities represent a major barrier to the reliable detection and attribution of global warming (Tebaldi et al., 2005;Mahlstein and Knutti, 2010).Furthermore, uncertainties in model simulations generally increase from global to regional scales because of uncertainty in regional climatic responses to global change (Hingray et al., 2007;Mariotti et al., 2011).Therefore, investigations of the spatial and temporal patterns of regional climate changes and regional climatic response mechanisms to global change are crucial for increasing the accuracy of models designed to detect and explain the causes of global climate change and predictions of future regional climate change.
The spatial heterogeneity of climate warming can be attributed to local climate factors and anthropogenic factors (Karl et al., 1991).For the local climate factors, determining factors such as cloud cover and precipitation (P ) can significantly influence the speed of regional warming (Hegerl and Zwiers, 2007;Lauritsen and Rogers, 2012).Spatial het-erogeneities in climate-factor trends have an important influence on various changes in the land surface energy balance.Studies have demonstrated that an increase in cloud cover can diminishes the surface solar radiation (R s ) and therefore reduces the daytime temperature (Dai et al., 1997;Zhou et al., 2010;Taylor et al., 2011), although it has the potential to increase nighttime temperature by intercepting outgoing longwave radiation (Campbell and VonderHaar, 1997;Shen et al., 2014).
Precipitation (P ) can alter the proportion of surface absorbed energy partitioned into sensible heat flux and latent heat flux; therefore, it has an inevitable effect on both land surface and near-surface air temperatures (Wang and Dickinson, 2012;Wang and Zhou, 2015).Additionally, P has a significant effect on soil thermal inertia and the response of surface vegetation, which results in an important feedback for regional and global warming (Seneviratne et al., 2010;Wang and Dickinson, 2012;Ait-Mesbah et al., 2015;Shen et al., 2015).
In addition to local climate factors, regional climate systems are significantly affected by the anthropogenic emissions of aerosols.Studies have indicated that improvements in air quality in recent decades over North America and Europe have led to a brightening effect (Vautard et al., 2009;Wild, 2012), whereas East Asia and India have experienced declines in R s (Xia, 2010;Menon et al., 2002;Wang et al., 2012Wang et al., , 2015)).Consequently, variations in R s may have an effect on both local and global climate change (Wild et al., 2007;Wang and Dickinson, 2013a).
Changes in land cover can also alter the energy exchange between the land surface and the atmosphere, and such changes have the potential to affect regional climates (Bounoua et al., 1999;Zhou et al., 2004;Falge et al., 2005).Previous studies have suggested that urbanization and other land-use changes contribute to promoting the warming effect caused by greenhouse gases (Kalnay and Cai, 2003;Lim et al., 2005;Chen et al., 2015).Overall, the effects of these factors on climate change may be very important at the regional scale and could lead to marked spatial differences in regional climate change; however, they are usually omitted from the detection and attribution of climate change at the global scale (Karoly and Stott, 2006).
China is a vast territory that has an abundance of climactic zones stretching from tropical to cold temperate, and a special alpine climate is observed over the Tibetan Plateau.Additionally, the dramatic economic development and explosive population growth in China in recent decades have caused significant changes in land cover and severe air pollution, including frequent haze events (Yin et al., 2017;Cheng et al., 2014;Wang et al., 2016).The climatic diversity and intensive human activity in this region will likely lead to a unique response to global warming with obvious spatial differences in climate change.Karl et al. (1991) analyzed the observational records for the period 1951-1989 and found that warming trends in China were faster than those of the United States but slower than those of the former Soviet Union.Several studies have revealed that the warming rate in Northwest China was approximately 0.33-0.39• C 10 yr −1 during the second half of the last century (Zhang et al., 2010;Li et al., 2012), which was significantly higher than the average warming rate over China of 0.25 • C 10 yr −1 (Ren et al., 2005) or the average global rate of 0.13 • C 10 yr −1 (Hegerl and Zwiers, 2007).The air temperature (T a ) over the Tibetan Plateau has increased by 0.44 • C 10 yr −1 over the last 30 years (Duan and Xiao, 2015), and this rate is considerably faster than the overall warming rate in the Northern Hemisphere (0.23 • C 10 yr −1 ) and worldwide (0.16 • C 10 yr −1 ) (Hartmann et al., 2013).To provide insights into global warming and improve the accuracy of future climate change predictions, understanding the characteristics and mechanisms of regional climate change is critical.
T a is a common metric for determining climate change on the global or regional scales.The land surface temperature (T s ) is also important in climate change research because of its direct relationship with the land surface energy budget.Previously, T s values used in regional climate research are primarily derived from satellite retrievals or reanalysis datasets (Weng et al., 2004;Peng et al., 2014), which both have satisfactory global coverage but questionable accuracy and integrity.Furthermore, satellite-derived T s values are only available under clear-sky conditions, thus limiting their applicability in climate change studies.
In China, both T s and T a are measured as conventional meteorological observation parameters by nearly all weather stations.An analysis of the spatiotemporal patterns of these parameters identified a close relationship between T s and T a , which indicates that T s and T a present equivalent accuracy when used to determine the characteristics of climate change.More importantly, T s is more sensitive than T a to the local land surface energy budget.
It is well known that the diurnal cycles in T a and T s are primarily determined by the surface energy budget.After sunrise, the surface absorbs solar radiation, and the surface net radiation becomes positive and heats the surface first.As a result, the air above the surface becomes unstable.Surface net radiation can be partitioned into three parts: ground heat flux, sensible heat flux, and latent heat flux.Ground heat flux heats the surface and stores energy during the daytime, and this energy may be re-emitted at night.Sensible heat flux directly heats the air above the surface.Latent heat flux is the energy employed to vaporize water during the surface water evaporation and vegetation transpiration processes.How surface net radiation partitions into ground heat flux, sensible heat flux, and latent heat flux is determined by both surface and atmospheric conditions (Wang et al., 2010a, b;Wang and Dickinson, 2012), i.e., surface wetness.Daytime surface net radiation is primarily determined by R s (Wang and Liang, 2008) and precipitation or surface wetness control partition of surface net radiation into latent and sensible fluxes (Wang and Liang, 2008).Therefore, it is expected that changes in R s and P play a key role in the variability of T s and T a (Hartmann et al., 1986;Wild, 2012;Manara et al., 2015).
However, quantitative assessments of the impact of R s on T s and T a are still lack due the paucity of high quality of long-term estimates of R s .In this study, we used sunshineduration-derived R s (Wang, 2014;Wang et al., 2015) to quantify the impact of R s on the spatial pattern of T a and T s .To our knowledge, this study presents the first analysis of the relationship between R s (and P ) and T a (and T s ) based on their spatiotemporal patterns and we further quantified the effect of variations in R s (and P ) on T a (and T s ) in China for the period 1960-2003.This article is organized as follows.Section 2 introduces the data and methods used in the study.Section 3 describes the spatial and temporal patterns of climate warming over China, analyses the effect of the variation in R s and P on T a and T s , and examines the spatial and temporal patterns of the warming trend of T a and T s after adjusting for the effects of R s and P , which eliminated the effects of R s and P on warming and highlighted the effects of large-scale warming caused by elevated concentrations of atmospheric greenhouse gases.Moreover, the spatial contrast in the warming trends of T a and T s in China was substantially reduced after adjusting for the effect of R s and P , and this result is consistent with the expectations under global warming.Finally, Sect. 4 presents a summary and discussion.
2 Data and methods

Data
The meteorological observational data used in this study are included recently released daily meteorological datasets, such as the China National Stations' Fundamental Elements Datasets V3.0 (CNSFED V3.0), and they were downloaded from China's National Meteorological Information Center (http://data.cma.cn/data/cdcdetail/dataCode/)(Cao et al., 2016).These datasets included observations of T s , T a , barometric pressure, relative humidity, and sunshine duration.All of the data acquisition and compilation of the climate variables were subjected to quality control measures.
As shown in Fig. 1, the number of stations used in this study (1977 stations selected from a total of 2479 stations) was significantly higher than that of previous studies (i.e., 57-852 stations) (Kukla and Karl, 1993;Shen and Varis, 2001;Liu et al., 2004;Li et al., 2015).Therefore, the observational data provided better spatial coverage and higher confidence in the detection of regional climate change than in previous studies (Fig. 1).
Observations of T s from weather stations are different from T s data retrieved via other approaches, such as satellite images and reanalysis.The T s observations were performed in 4 × 2 m bare land plots proximal to the weather stations.The surface of the observational fields was loose, grassless and flat, and at the same level as the ground surface of the weather station.Three thermometers, including a surface thermometer, a surface maximum thermometer, and a surface minimum thermometer were placed horizontal to the surface of the observational field, with half of each thermometer embedded in the soil and the other half exposed to the air.When the observational field was covered by snow, the thermometers were placed on the snow surface.Additionally, the exposed parts of the thermometers were cleaned to remove dust and dew.
We verified the reliability of the T s observational records by analyzing the relationship between T a and T s during 1960-2003.As shown in Fig. S1 in the Supplement, the mean Pearson correlation coefficients between daily maximum land surface temperature (T s-max ) and daily maximum air temperature (T a-max ) calculated from the monthly anomalies were 0.775, 0.843, and 0.806 for the annual, warm, and cold seasonal scales, respectively, and these values were statistically significant (99 % confidence level) for all stations.The mean correlation coefficients between the daily minimum land surface temperature (T s-min ) and daily minimum air temperature (T a-min ) were 0.861, 0.842, and 0.865 for the annual, warm, and cold seasonal scales, respectively, and these values were statistically significant (99 % confidence level) for all stations.The high correlations indicated that observations of either T s or T a could be used for climate change detection.
The most fundamental energy resource for T s and T a is R s .In most previous studies, the observed R s have been used to analyze the relation between the variation in R s and T a over China.However, fewer sites were used for R s observations than for other climatic variables; for example, only 85 sites were used for R s observations in Liu et al. (2004) and only 90 sites were used in Li et al. (2015).
Importantly, sensitivity drift the instruments used for the R s observations led to a faster dimming rate before 1990, and instrument replacements from 1990 to 1993 resulted in a false sharp increase in R s (Wang, 2014;Wang et al., 2015).The limited distribution and low quality of R s observations have impeded the wide scientific application of this parameter.
Therefore, we used sunshine-duration-derived R s , which is based on an effective hybrid model developed by Yang et al. (2006).This model has subsequently been improved (Wang, 2014;Vose et al., 2005) and it has performed well in regional and global applications (Tang et al., 2011;Wang et al., 2012).Sunshine-duration-derived R s not only accurately reflects the effects of clouds and aerosols on the R s but also more exactly reveals long-term trends (Wang, 2014;Wang et al., 2015).Additionally, sunshine-duration-derived R s values are better correlated with the satellite retrievals, reanalysis, and climate model simulations than R s values observed from observation (Wang et al., 2015).
The data are collected by a total of 2474 meteorological stations; however, the lengths of the effective observation records for the stations are different.Additionally, only a small number of stations were installed before 1960, and the observational records of T s at many stations were anomalous after 2003 because of automation.Therefore, in our analysis, we selected 1977 meteorological stations (see Fig. 1) for which the observation records with valid data were longer than 30 years during the 43 years between 1960 and 2003.
The monthly anomalies relative to the 1961-1990 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 (Li et al., 2015;Sun et al., 2016).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 the next April.

Methods
As shown in Fig. 1, the spatial distribution of the weather stations throughout China is extraordinarily asymmetric and the density of weather stations in east China is far greater than that in west China.We used the area-weight average method to reduce these biases when calculating the national mean.First, we divided the study region into 1 • × 1 • grids (see Fig. S2) for a total 953 grids covering China.Second, we assigned all selected stations to the grids, and this resulted in 627 grids containing stations, which accounted for 65.79 % of the total.Finally, the grid box value was the average of all stations in the grid, and the national mean was the area-weight average of all effective grids (Jones and Moberg, 2003).
The linear trends reported in this study were calculated via linear regression based on the monthly anomalies of T , R s , and P .Two national mean trends were calculated from the anomalies of the grids.In the first method (Method I), the national mean monthly anomalies were calculated using the area weight of each grid first, and then the national mean trend based on the time series of the national average anomalies was calculated.In the second method (Method II), the trend at each grid was calculated first, and then the national mean trend was calculated from the grid trends.
In our study, we calculated the national mean trends of the temperatures using Method I and II because both methods have been used in previous studies (Gettelman and Fu, 2008).The results for the two methods are expected to be the same when the time series of all grids is integrated and data are not missing (Zhou et al., 2009); however, when data are missing, small differences may occur (See Table 1).As shown in Table 1, the absolute value of the difference between Method I and Method II ranged from 0.011 to 0.033 • C 10 yr −1 , which represents 3.4 to 14.3 % of the trends (using the results of Method I as the reference).For purposes of clarification, the trends derived from Method I are discussed in the main text, whereas the results from both methods are shown in Table 1.
The effect of R s and P on T s-max and T a-max was determined via a multiple linear regression (Roy and Haigh, 2011) of the monthly anomalies using the following equation: where T represents the monthly anomalies of T s-max , T s-min , T a-max , and T a-min ; R s and P represent the monthly anomalies of surface solar radiation and precipitation, respectively; S R s and S P are the sensitivities of temperatures to R s and P , respectively; c is the constant term; and ε indicates the residuals of the equation, which represents the contribution from other factors such as longwave radiation flux and internal variability.The coefficients of determination (R 2 ) for the multilinear regression equation (Eq. 1) are shown in Fig. S3, and they indicate the portion of the variance of T that could be attributed to that of R s and P .High coefficients of determination were obtained, which showed that the linear regression performed well, particularly for South China and the North China Plain.To separate the contributions of R s and P , we further calculated the partial correlation coefficients between R s and T (or P and T ), which are shown in Figs.S4 and S5.
To determine the effect of R s and P on the analyzed temperatures, we removed their effects from their original time series of T s-max and T a-max based on the multilinear relationship calculated in Eq. ( 1).Then, we calculated the trends from both the original and adjusted time series.By comparing the derived trends of the original and adjusted time series, we quantitatively assessed the effect of R s and P on T s-max Table 1.Warming rates (unit: • C 10 yr −1 ) of the temperatures (T s-max , T s-min , T a-max , T a-min ) for the annual, warm and cold seasonal scales."Raw" and "Adjusted" represent the warming rates calculated for the data before and after adjusting for the effect of surface solar radiation (R s ) and precipitation (P ), respectively.In Method I, the national mean anomalies were calculated first and then the national mean trend based on this time series was calculated.In Method II, the trend of each grid was calculated first and then the national mean value of the trends of all grids was calculated using the area-weight average method.We calculated the national mean trends of the temperatures using both methods.and T a-max , particularly for the spatiotemporal pattern of their trends.

Results
3.1 Trends of surface temperature and air temperature

Temporal patterns in temperature variability
The long-term changes in T s-max and T a-max and T s-min and T a-min from 1960 to 2003 are shown in Figs. 2 and 3, respectively.In addition to the annual variability (Figs.2a and  3a), the temperature variability in both warm seasons (May-October; Figs.2b and 3b) and cold seasons (November to the following April; Figs.2c and 3c) were analyzed.In the annual records, all temperatures exhibited an obvious warming trend throughout China (Figs. 2a and 3a).
As shown in Table 1, the national mean warming rate from 1960 to 2003 for T s-max was 0.227 ± 0.091 • C 10 yr −1 (95 % confidence level) and T a-max was 0.167 ± 0.068 • C 10 yr −1 (95 % confidence level).The warming rate of T a-max based on the 1977 stations examined in the current study was slightly higher than the global average (0.141 • C 10 yr −1 ) from 1950 to 2004 (Vose et al., 2005) and the rate obtained from a previous analysis (0.127 • C 10 yr −1 ) of temperatures from 1955 to 2000 based on 305 stations in China (Liu et al., 2004).Additionally, the increases in T s-max and T a-max in the cold seasons were much larger than those in the warm seasons, which is consistent with previous studies of China and other regions (Vose et al., 2005;Ren et al., 2005;Shen et al., 2014).
Similarly, the warming rates of T s-min and T a-min in the warm seasons were also clearly lower than those in the cold seasons.As shown in Fig. 3, T s-min increased by 0.315 ± 0.058 • C 10 yr −1 (95 % confidence level) and T a-min increased by 0.356 ± 0.057 • C 10 yr −1 (95 % confidence level) (see Fig. 3a) from 1960 to 2003.The warming trend of T a is generally consistent with earlier studies (Liu et al., 2004;Shen et al., 2014;Li et al., 2015); however, these trends are considerably larger than the rates reported for the global average (0.204 • C 10 yr −1 ) (Vose et al., 2005).For the seasonal scales, the warming rate of T s-min and T a-min in the cold seasons was almost double that of the warm seasons from 1960 to 2003 (see Table 1).
The warming rate of T s-min (T a-min ) was significantly faster than that of T s-max (T a-max ) and the warming rates of all temperatures in the cold seasons were substantially greater than those in the warm seasons (Easterling et al., 1997;Liu et al., 2004;Li et al., 2015).Although previous studies have indicated that the microclimate (e.g., urban heat island) has a larger effect on minimum temperatures because of the lower and more stable boundary layer at night (Christy et al., 2009;Zhou and Ren, 2011), many investigators argue that variability in R s is the primary reason for the daily contrast in warming rates (Makowski et al., 2009;Sanchez-Lorenzo and Wild, 2012).

Spatial patterns in temperature variability
As shown in Fig. 4, clear spatial heterogeneity was demonstrated in the warming rates for T s-max and T a-max in China from 1960 to 2003.The trends of T s-max and T a-max were sta-  tistically higher for the Tibetan Plateau as well as Northwest and Northeast China (see Fig. S6) compared with the North China Plain and South China.Cooling trends in T s-max even detected for the Sichuan Basin, the Yangtze River Delta, and the Pearl River Delta.Lower rates of warming of T a-max in South China and the North China Plain have also been previously reported (Liu et al., 2004;Li et al., 2015).
The warming rates of T s-max and T a-max in South China and the North China Plain in the warm seasons were considerably lower than those in the cold seasons, which resulted in stronger spatial heterogeneity in the warm seasons (Fig. 4b  and h).The spatial and seasonal patterns of T a-max were similar, although they were not as similar as the patterns of T s-max .The spatial contrast in the trends between T s-min and T a-min was much less than that between T s-max and T a-max , although a strong dependence on latitude was observed (Fig. 4d and j).Related studies suggested that this dependence was strongly associated with the mode variability in large-scale circulation, such as a negative trend in the North Atlantic Oscillation during this period (Wallace et al., 2012;Ding et al., 2014).
The correlation between T s and T a was highly significant.Based on the time series of the national mean yearly anomalies (see Figs. 2 and 3), the correlation coefficient between T s-max and T a-max was 0.877 and between T s-min and T a-min was 0.976 on the annual scale.In the spatial pattern of the trends (Fig. 4), the correlation coefficient between T s-max and T a-max was 0.488 and between T s-min and T a-min was 0.638 on the annual scale.All of these correlations between T s and T a were significant at the 95 % significance level, which indicated a close relation between T s and T a for both interannual fluctuations and secular trends.
The correlation between T s-min and T a-min was significantly higher than that between T s-max and T a-max .T s-min is closely related to the land-atmosphere longwave wave radiation balance at night, which is closely associated with the atmospheric greenhouse effect (Dai et al., 1999).During the day, T s is directly determined by the land surface energy balance, i.e., the incoming energy (including R s ) and atmospheric longwave radiation (Wang and Dickinson, 2013b), and it is partitioned into latent and sensible heat fluxes (Zhou and Wang, 2016).Although T a is dependent on the landatmosphere sensible heat flux, it is also affected by local and/or large-scale circulation.Thus, the changes in the land surface energy balance caused by R s and P have different levels of effect on T s and T a during the day, which most likely caused the lower correlation between T s-max and T a-max than that between T s-min and T a-min .
3.2 Effect of R s and P on temperatures

Effect of R s
As shown in Fig. S4, R s is closely linked with T s-max and T a-max but not with T s-min and T a-min , and the correlation between T s-max and R s was higher than that between T a-max and R s .For the seasonal scales, the partial correlation between T s-max and T a-max and R s in the warm seasons was stronger than that in the cold seasons, and this correlation was stronger in South China and the North China Plain.South China has high soil moisture; therefore, the relationship between the energy used for evapotranspiration and R s is approximately linear (Zhou et al., 2007;Wang and Dickinson, 2013a).However, northwest China presents dry soil over most of the year; thus, the energy used for evapotranspiration is more dependent on P .As a result, the energy available for heating the surface and air temperatures is not as closely correlated with R s .Therefore, the correlation co-  (d-f)), daily maximum air temperature (T a-max , (g-i)), and daily minimum air temperature (T a-min , (j-l)) for the annual, warm (May-October), and cold (November-next April) seasonal scales.All trends reported here were calculated using a linear regression based on the least-squares method.
efficients between R s and T s-max and T a-max were lower in the northwest China.
To quantify the effect of R s on temperature, the sensitivity of the studied temperatures to changes in R s was calculated (Eq.1).As shown in Fig. S7, T s-max was the most sensitive to R s , followed by T a-max , and the national means for T s-max was 0.092 ± 0.018 • C (W m −2 ) −1 (95 % confidence level) and T a-max was 0.035 ± 0.010 • C (W m −2 ) −1 (95 % confidence level).T s-min and T a-min were not sensitive to R s because these temperatures are primarily affected by atmospheric longwave radiation at night.
Based on the above analysis, we calculated the effect of changes in R s on the studied temperatures.From 1960 to 2003, the calculations of the monthly anomalies at 1977 stations indicated that the national mean rate of decrease in R s was −1.502 ± 0.42 W m −2 10 yr −1 (95 % confidence level), and the trend was significant in most regions of China (see Fig. S8).Our rate of decrease was considerably less than the global average diminishing rate (form approximately −2.3 to −5.1 W m −2 10 yr −1 ) between the 1960s and the 1990s (Gilgen et al., 1998;Stanhill and Cohen, 2001;Liepert, 2002;Ohmura, 2006) and the national mean dimming rate across China (from approximately −2.9 to −5.2 W m −2 10 yr −1 ) between the 1960s and the 2000s based on radiation station observations (Che et al., 2005;Liang and Xia, 2005;Shi et al., 2008;Wang et al., 2015).
As noted in the data section, the sensitivity drift and replacement of instruments used for the R s observations resulted in a significant homogenization of the station observation records (Wang, 2014;Wang et al., 2015), which introduced considerable uncertainty into the trend estimations.Tang et al. (2011) used quality-controlled observational data from 72 stations and two radiation models based on 479 stations to determine that the rate in China decreased from approximately −2.1 to −2.3 W m −2 10 yr −1 during 1961-2000, and they also showed that R s values have remained essentially unchanged since 2000.These findings are generally consistent with our results.
Because of the decreasing trend in R s , the national mean warming trends of T s-max and T a-max decreased by 0.139 • C 10 yr −1 and 0.053 • C 10 yr −1 , respectively.Spatially, the decreasing rate of R s in South China and the North China Plain was significantly higher than that in other regions, particularly in the warm seasons (Fig. 5b).Therefore, the cooling effect of decreasing R s on T s-max and T a-max was more significant in South China and the China North Plain, and it resulted in significantly lower warming rates of T s-max and T a-max in those regions than in the other regions (see Fig. 4).The spatial consistency between the decreasing R s trend and the slowdown of T s-max and T a-max warming implies that variations in R s were the primary reason for the spatial heterogeneity of the warming rate in T s-max and T a-max .

Effect of P
As shown in Fig. S5, a significant negative correlation was detected between T s-max and P , and the correlation was more significant in the warm seasons than in the cold seasons.P negatively correlated with temperature because P reduces temperatures by increasing the surface evaporative cooling (Dai et al., 1997;Wang et al., 2006).The national mean sensitivities of T s-max and T a-max to P were −0.321 ± 0.098 • C 10 mm −1 and −0.064 ± 0.054 • C 10 mm −1 (95 % confidence level), respectively.As shown in Fig. S9, seasonal and spatial changes in the sensitivity of T s-max and T a-max to P were apparent (Fig. S9a-c and g-i).The sensitivities of T s-max and T a-max were significantly higher in arid regions (dry seasons) than humid regions (rainy seasons) (Wang and Dickinson, 2013a).As expected, T s-min and T a-min were both less sensitive to variations in P .
The trend in P from 1960 to 2003 over the 1977 stations showed obvious spatial heterogeneities.A slight increasing trend in P was observed in China during this period at rate of 0.112 ± 0.718 mm 10 yr −1 (95 % confidence level).An increasing P trend was observed in northwestern China and southeastern China, whereas a decreasing trend was observed in the North China Plain, the Sichuan Basin, and parts of northeastern China.However, the P trends were not significant in most regions (see Fig. S8).Variations in P significantly differed by season (see Fig. 6b and 6c).The seasonal and spatial variations in P are consistent with those of previous studies (Zhai et al., 2005;Wang et al., 2015).
For T a-max and T s-max , the warming trend in the North China Plain, the Sichuan Basin, and parts of northeastern China was aggravated by the reduction in P , whereas the warming trends in northwestern China and in the Mongolian Plateau were slowed by increases in P (Fig. 6d).For the national average, the effect of increasing P resulted in decreases in the warming trends of T s-max and T a-max by −0.007 and −0.002 • C 10 yr −1 , respectively.However, the effect of P on T s-max was approximately an order of magnitude less than that of R s .

Trends of surface and air temperature after
adjusting for the effect of R s and P Based on the above analysis of the effect of R s and P on temperatures, we found that variations in R s and P had little effect on T s-min and T a-min .However, R s and P had an important effect on the trends of T s-max and T a-max (see Fig. S3), particularly in central and South China, where T s-max and T a-max were more closely related to R s (see Fig. S4).Therefore, only the effects of R s and P on T s-max and T a-max were analyzed.After adjusting for the effect of R s and P (Fig. 7), the warming rates of T s-max and T a-max increased by 0.146 • C 10 yr −1 (64.3 %) and 0.055 • C 10 yr −1 (33.0 %), respectively.Additionally, the increasing amplitude of warming rates in the warm seasons was significantly higher than that in the cold seasons, which resulted in a seasonal contrast in warming rates, with T s-max and T a-max decreasing by 45.0 and 17.2 %, respectively (see Table 1).More importantly, after adjusting for the effect of R s and P , the spatial coherence of the warming rates of T s-max and T a-max in South China and the North China Plain clearly improved (Fig. 8).The regional differences among the North China Plain, South China, and other regions in China significantly decreased because of the increase in warming rates in South China and the North China Plain.Additionally, the warming trends of T s-max and T a-max became more statistically significant in the North China Plain and South China (see Fig. S10).
To clearly illustrate these changes, we selected two regions in China for further investigation: R1 primarily included the North China Plain and R2 primarily included the Loess Plateau (see Fig. 9a).Although these regions share the same latitudes, the trend for R s were substantially different (see Fig. 9b).After adjusting for the effect of R s and P , the annual trends for T s-max and T a-max in R1 increased by 0.304 and 0.118 • C 10 yr −1 , respectively, whereas those in R2 increased by only 0.025 and 0.016 • C 10 yr −1 , respectively.Therefore, after the adjustment, the contrasts in the warming rates of T s-max and T a-max between R1 and R2 were significantly reduced (see Fig. 9d).
After the adjustment in R1, the seasonal and diurnal contrasts in the warming rates of T s-max and T a-max significantly decreased.The contrasts in warming rates between the warm seasons and cold seasons decreased by 68.7 % for T s-max and by 50.8 % for T a-max after the adjustment.Additionally, the contrasts in the warming rates between T s-max and T s-min decreased by 93.4 % and between T a-max andT a-min decreased by 59.6 % in R1.In R2, the adjustment did not significantly change the seasonal and diurnal contrasts in temperatures.Overall, the trends for R1 and R2 became more consistent after adjusting for difference in R s and P (see Fig. 9d).(d-f) and the third line (g-i) are the trend changes caused by secular variations in R s on T s-max and T a-max .Equation ( 1) was used to strip away the effect of R s on temperatures, and we calculated the trend difference ( Trend, (d-i)) between the time series of temperatures before and after adjusting for the effect of R s .Finally, the effect of R s on the trends of T s-max and T a-max was quantified and analyzed (Sect.3.2.1).(d-f) and the third line (g-i) are the trend changes caused by secular variations in P on T s-max and T a-max .We used Eq. ( 1) to remove the effects of P on the temperatures, then calculated the trend difference ( Trend, (d-i) between the time series of temperatures before and after adjusting for the effect of P .Finally, the effect of P on the trends of T s-max and T a-max was quantified and analyzed (Sect.3.2.2).  1) to simultaneously adjust for the effects of surface solar radiation (R s ) and precipitation (P ) on T s-max and T a-max and then analyzed the changes in the interannual variation of T s-max and T a-max (Sect.3.3).

Conclusions and discussion
Although a general warming trend has been observed throughout China, the regional warming trends show significant spatial and temporal heterogeneity.In this study, we analyzed the spatial and temporal patterns of T s and T a from 1960 to 2003 and further analyzed and quantified the effects of R s and P on these temperatures.The primary results of the study are as follows.
The national mean warming rates from 1960 to 2003 of T s-max , T s-min , T a-max , and T a-min were 0.227 ± 0.091, 0.315 ± 0.058, 0.167 ± 0.068, and 0.356 ± 0.057 • C 10 yr −1 , respectively.The warming rates of T s-max and T a-max in South China and the North China Plain were significantly lower than those in the other regions, and the spatial heterogeneity in the warm seasons was greater than that in the cold seasons.
During the study period, the R s value decreased by −1.502 ± 0.042 W m −2 10 yr −1 (95 % confidence level), and higher dimming rates were observed in South China and the North China Plain.Using a partial regression analysis, we found that R s plays a distinctly important role in the spatial warming patterns of T s-max and T a-max .
After adjusting for the effect of R s and P , the warming rates of T s-max and T a-max in South China and the North China Plain significantly increased and the regional differences in warming rates in China clearly decreased (see Fig. 8).After the adjustments, the warming rates of T s-max and T a-max in the North China Plain increased by 0.304 and 0.118 • C 10 yr −1 , respectively, whereas those on the Loess Plateau increased only by 0.025 and 0.016 • C 10 yr −1 , respectively.Therefore,  (b, d, f)) for the annual, warm, and cold seasonal scales after adjusting for the effects of surface solar radiation (R s ) and precipitation (P ).We used Eq. ( 1) to simultaneously adjust the effects of R s and P on T s-max and T a-max and then analyzed the changes in the secular trends of T s-max and T a-max (Sect.3.3).
the differences in warming rates of T s-max and T a-max between the North China Plain and the Loess Plateau were almost eliminated (see Fig. 9d).
After adjusting for the effect of R s and P , the warming trend of T s-max increased by 0.146 • C 10 yr −1 and that of T a-max increased by 0.055 • C 10 yr −1 .In addition, the trends of T s-max and T a-max became 0.373 ± 0.068 and 0.222 ± 0.062 • C 10 yr −1 , respectively.Reduction in R s resulted in decreases in the warming rates of T s-max and T a-max by 0.139 • C 10 yr −1 and 0.053 • C 10 yr −1 , respectively, which accounted for 95.0 and 95.8 % of the total effect of R s and P , respectively.For the seasonal contrast, the warming rates of T s-max and T a-max decreased by 45.0 and 17.2 %, respectively.For the daily contrast, the warming rates of T s and T a decreased by 33.0 and 29.1 %, respectively.
In addition to R s and P , temperature warming rates may be affected by many other factors, such as land cover and land use changes; however, those factors have not been discussed in this study because of lack of data (Liu et al., 2005;Zhang et al., 2016).After adjusting for the effect of changes in R s and P changes, the spatial differences in the warming trends clearly decreased; however, certain regional differences remained.The warming rate of T s-max in the Sichuan Basin remained significantly lower than that in other regions after adjusting for these effects.Additionally, the differences in the warming rates of T s-min and T a-min between the northern and southern areas were not explained by the effects of R s and P ; further study is required.>. -"'"

Figure 1 .
Figure 1.Elevation map of mainland China and spatial distribution of the 1977 meteorological stations used in this study.The datasets were provided by China's National Meteorological Information Center (You et al., 2016) (http://data.cma.cn/data/cdcdetail/dataCode/).

Figure 2 .
Figure 2. National mean yearly anomalies of daily maximum land surface temperature (T s-max , blue line) and daily maximum air temperature (T a-max , red line) for the annual (a), warm (b), and cold (c) seasonal scales for the reference period from 1961 to 1990.

Figure 3 .
Figure3.National mean yearly anomalies of daily minimum land surface temperature (T s-min , blue line) and daily minimum air temperature (T a-min , red line) for the annual (a), warm (b), and cold (c) seasonal scales for the reference period1961-1990.

Figure 4 .
Figure 4. Maps of the trends of the monthly anomalies for daily maximum land surface temperature (T s-max , (a-c)), daily minimum land surface temperature (T s-min ,(d-f)), daily maximum air temperature (T a-max , (g-i)), and daily minimum air temperature (T a-min , (j-l)) for the annual, warm (May-October), and cold (November-next April) seasonal scales.All trends reported here were calculated using a linear regression based on the least-squares method.

Figure 5 .
Figure 5. Maps of the trends in surface solar radiation (R s , (a-c)) and their effect on the warming rates of daily maximum land surface temperature (T s-max , (d-f)) and daily maximum air temperature (T a-max , (g-i)).The first line (a-c) is the trends of R s from 1960 to 2003; the second line(d-f) and the third line (g-i) are the trend changes caused by secular variations in R s on T s-max and T a-max .Equation (1) was used to strip away the effect of R s on temperatures, and we calculated the trend difference( Trend, (d-i)) between the time series of temperatures before and after adjusting for the effect of R s .Finally, the effect of R s on the trends of T s-max and T a-max was quantified and analyzed (Sect.3.2.1).

Figure 6 .
Figure 6.Maps of the trends in precipitation (P ) (a-c) and their effect on the warming rates for daily maximum land surface temperature (T s-max , (d-f)) and daily maximum air temperature (T a-max , (g-i)).The first line (a-c) is the trends of P during 1960-2003; the second line(d-f) and the third line (g-i) are the trend changes caused by secular variations in P on T s-max and T a-max .We used Eq.(1) to remove the effects of P on the temperatures, then calculated the trend difference( Trend, (d-i)  between the time series of temperatures before and after adjusting for the effect of P .Finally, the effect of P on the trends of T s-max and T a-max was quantified and analyzed (Sect.3.2.2).

Figure 7 .
Figure 7. Regional average anomalies of daily maximum land surface temperature (T s-max , blue line) and daily maximum air temperature (T a-max , red line) for the annual (a), warm (b), and cold (c) seasonal scales for the reference period from 1961 to 1990.We used Eq.(1) to simultaneously adjust for the effects of surface solar radiation (R s ) and precipitation (P ) on T s-max and T a-max and then analyzed the changes in the interannual variation of T s-max and T a-max (Sect.3.3).

Figure 8 .
Figure 8. Maps of the trends of the monthly anomalies for the daily maximum land surface temperature (T s-max , (a, c, e)) and daily maximum air temperature (T a-max ,(b, d, f)) for the annual, warm, and cold seasonal scales after adjusting for the effects of surface solar radiation (R s ) and precipitation (P ).We used Eq.(1) to simultaneously adjust the effects of R s and P on T s-max and T a-max and then analyzed the changes in the secular trends of T s-max and T a-max (Sect.3.3).

-
Figure 9. (a) Maps of the trends of surface solar radiation (R s ) and the location of the regions selected for further analysis: R1 (latitude: 30-40 • N; longitude: 110-120 • W) and R2 (latitude: 30-40 • N; longitude: 100-110 • W).(b) National mean trends for R1 and R2.(c) Annual, warm, and cold seasonal-scale trends calculated based on the data before adjusting the effects of R s and P .(d) Annual, warm, and cold seasonal-scale trends calculated based on the data after adjusting the effect of R s and P .All error bars indicate the 95 % confidence interval.