The study area (83.87°E–112.07°E and 21.14°N–36.49°N) is located in Southwest China, including Guangxi, Guizhou, Chongqing, Yunnan, Sichuan, southern Qinghai and some parts of the Tibet Autonomous Region, with a total area of 2.33 × 10⁶ km² (Xu et al., 2020) (Figure 1). The elevation in Southwest China presents a considerable variation with a difference of ∼8,000 m, exhibiting complicated topographic structures (e.g., mountains, plateaus, hills, basins, and plains). Due to the wide range of latitudes and complex topography, Southwest China is characterized by a variety of climates and environments, ranging from monsoon region in the southeast zone to semi-arid region in the northwest zone (Jin and Wang, 2016). There are various ecosystems, including tropical rain forest, tropical seasonal rain forest, subtropical evergreen broad-leaved forest, and alpine vegetation (Gao et al., 2018).

The ANUSPLIN is created using thin-plate smoothing splines, which make it suitable for interpolating climate data with large noises (Hutchinson and Gessler, 1994; Price et al., 2000; Guo et al., 2020). The noisy multivariable climatic data are treated as a function with one or more independent variables during the fitting process of a climatic surface, and thus can produce mean error lower than other interpolation methods (Islam and Déry, 2017; Zhao et al., 2019; Cheng et al., 2020). The theoretical statistical model is expressed as: Z i = f ( x i ) + b T y i + e i ( i = 1 , … , N ) (1) where Z _i represents the predicted value at location i; x _i is the spline independent variable as a multidimensional vector, and f represents a smoothing function of x _i which needs to be estimated; y _i is the independent covariable as a multidimensional vector, and b is the unknown coefficients for the y _i ; n is the number of observational data. Each e _i is an independent, zero mean error term with variance w _i σ², where W _i is the known relative error variance and σ² is the error variance which is constant across all data points.

The traditional ANUSPLIN treats longitude and latitude as independent variables, with elevation as a covariate. To optimize the ANUSPLIN model, we improved it by incorporating more terrain-related factors, with slope, and aspect also as covariates (hereafter called M-ANUSPLIN).

The NSAT parameters in Southwest China are estimated using co-kriging, ANUSPLIN and M-ANUSPLIN models, respectively. Due to the complex topography of the domain over Southwest China, we tested the performance of these models using different combinations of covariates (Table 4) to identify the best model in this region, and to find which variable contributes more to NSAT variability.

To evaluate these models, a 10-fold cross validation test was conducted to assess the overall error of the interpolated NSAT grid. The advantage of 10-fold cross validation is that all observations are used for both training and validation, and each observation is used for validation exactly once. Thus, it is widely used to validate gridded observations (Appelhans et al., 2015; Yoo et al., 2018). In 10-fold cross validation for this study, the original observation data of 494 meteorological stations were randomly partitioned into ten subsamples. Of the ten subsamples, a single subsample was retained as the validation data for testing the model, and the remaining nine subsamples were used as training data. The cross validation process was then repeated ten times, with each of the ten subsamples used exactly once as the validation data. Hence, 10 different combinations of training and test sets were formed, and each of training and test pair was applied and evaluated. Final evaluation of 10-fold cross validation test was determined by the average mis-classification probability over the ten test sets to produce a single estimation.

Based on the results of 10-fold cross validation test, the statistical indices of Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and R-squared (R ²) between predicted and observed values were selected as interpolation performance evaluation criteria. Briefly, MAE provides a measure of how far the estimate can be in error; RMSE provides a measure that is, sensitive to outliers; and R ² provides the proportion of variation that is, explained by the predictor variables. The performance and bias were then compared by using the three indices, and the interpolation method with better performance was further selected. The calculation formulas of them are shown below: M A E = 1 n ∑ i = 1 n | P i − O i | (2) R M S E = 1 n ∑ i = 1 n ( P i − O i ) 2 (3) R 2 = 1 − ∑ i = 1 n ( P i − O i ) 2 ∑ i = 1 n ( O i − O ¯ ) 2 (4) where P _i and O _i are the estimated NSAT and original observational NSAT at each station, respectively;   O ¯ is the mean of observational NSAT; and n is the sample number.

To further examine the accuracy of obtained data set, two published and widely used air temperature products were compared. One is the 30 arc-second resolution grid product obtained from WorldClim Data Portal (hereafter called WorldClim2), the other is the gridded data sets cover China at 1 km × 1 km resolution (hereafter called HMTC). Mean temperature values in time series of 1970–2000 derived from all datasets were compared. The M-ANUSPLIN and WorldClim2 values were resampled to 1 km × 1 km to render them consistent with the HMTC reference products.

A trend slope ratio analysis was examined to investigate the changing trends of NSAT for each pixel in 1969–2018 at both annual and seasonal scales. The formula is as follows (Vogelsang and Nawaz, 2017): s l o p e = n × ∑ i = 1 n i × T e m i − ∑ i = 1 n i ∑ i = 1 n T e m i n × ∑ i = 1 n i 2 − ( ∑ i = 1 n i ) 2 (5) where slope is the degree of change in Tem; n is the number of studied years; i is the order of year from 1 to 50 in the study period; and Tem _i is the average Tem in the ith year. Slope >0 means that the air temperature over n years increased (warming trend); while slope<0 signifies a decreasing trend (cooling trend). To test the significance of these trends, a significance test (F-test) was applied. According to the F-test, the trends were divided into categories of extremely significant (p < 0.01), significant (0.01 < p< 0.05), and non-significant level (p > 0.05).

A coefficient of variation (CV) index was also considered to evaluate the spatio-temporal variation of NSAT. The CV was calculated as follows (Yang and Jiang, 2017): C V = 1 T e m ¯ × 1 n − 1 ∑ i = 1 n ( T e m i − T e m ¯ ) 2 × 100 % (6) where CV is the coefficient of variation for NSAT; T e m ¯ is the mean temperature; Tem _i is the temperature for year i; n is the number of studied years; i is the order of year from 1 to 50 in the study period. The significance test is carried out based on the p value. The CV is classified into three categories of weak variation (0 < CV ≤ 10%), medium variation (10% < CV ≤ 100%) and strong variation (CV > 100%).

For a long-term climatic dataset, it is expected to experience multiple changes rather than a single break (Khapalova et al., 2018). To detect and identify the time when significant changes for NSAT happened in the time series of 1969–2018, pruned exact linear time (PELT) was also conducted in this study. The superiority of PELT exists in the ability of accurate and fast detection and identification of multiple change-points (Killick et al., 2012).

Table 2 compares the MAE, RMSE and R ² of the predicted NSAT parameters using different models over Southwest China in 1969–2018. Apparently, the prediction of M-ANUSPLIN model gave lower MAE and RMSE compared to the ANUSPLIN and co-kriging models for both annual and seasonal NSATs. The R ² values of the M-ANUSPLIN model were 0.077–0.328°C above those of the other two models. Specifically, for annual parameters, relative to the ANUSPLIN and co-kriging models, the MAEs and RMSEs for the M-ANUSPLIN model were 0.02–0.04°C lower for mean temperature (T _mean); 0–0.02°C lower for maximum temperature (T _max); and 0.03–0.07°C lower for minimum temperature (T _min). For seasonal parameters, the MAEs and RMSEs for the M-ANUSPLIN model were also 0.01–0.04°C and 0–0.05°C lower than the ANUSPLIN and co-kriging models. These results indicate the improvement of the M-ANUSPLIN model by incorporating slope and aspect as interpolators. Moreover, Table 2 also reflects the seasonal dependence of the model performance. In general, the MAEs and RMSEs decline in the order from winter to autumn to spring then to summer, while the R ² basically has the opposite trend.

Figure 2 and Table 3 show the performances of co-kriging, ANUSPLIN and M-ANUSPLIN models along altitudinal gradients. It can be seen that in general, the M-ANUSPLIN model has obvious advantages in regions with altitude <4,000 m, followed by the ANUSPLIN and co-kriging models. However, with the exception of elevation >4000 m, the MAE and RMSE of the M-ANUSPLIN model were slightly higher than those of the ANUSPLIN model. Moreover, in areas with altitude >4,000 m, all the three models generally gave overestimated values.

A summary of the performances of the ANUSPLIN and M-ANUSPLIN models were further compared using different combinations of covariates (Table 4). It can be seen that the incorporation of additional terrain-related factors resulted in more accurate results, and slope angle contribute more to NSAT than slope orientation.

To further probe into the spatial variation patterns of annual and seasonal NSATs over Southwest China, the variations with significance test were analyzed at the pixel scale. The area proportions occupied by extremely significant, significant and non-significant NSAT related indices are shown in Figures 6, 7. It can be seen that the M-ANUSPLIN interpolation datasets can capture the detailed NSAT very well, and they can accurately represent the climate characteristics in Southwest China, such as the extremely significant cooling regions with high elevations (e.g., the northwest of Qinghai-Tibet Plateau), the non-significant warming regions with low elevations (e.g., the northeast of Sichuan Basin), and extremely significant warming regions in most part of the study area. Moreover, despite slight differences between the spatial variations of annual NSAT and those of the seasonal indices, they exhibited highly consistent characteristics. Specifically, for the mean annual NSAT, during the period 1969–2018, 85.66% of the study area showed extremely significant warming with the most notable increases in the southeast region. Areas occupied by non-significant warming accounted for 2.48%, followed by significant warming with 1.32%, leading to an overall warming tendency across Southwest China. Nevertheless, the northwest region experienced a contrary characteristic, with extremely significant, and significant cooling accounted for 8.35%. This could be attributed to the high altitude in this region, which belongs to the Qinghai-Tibetan Plateau. Regarding the mean seasonal NSAT, areas occupied by extremely significant warming in spring, summer, autumn and winter accounted for 63.48, 71.47, 72.36 and 54.77%, respectively. Otherwise, it can be seen that winter has the highest proportions for non-significant warming, comparing to other seasons (Figure 7).

Figure 8 shows the coefficient maps of variation (CV) index for NSAT over Southwest China. It indicates the high consistency of the CV index of annual NSAT with that of seasonal index over the past 50 years. Both annual and seasonal CV indices were generally stable and mainly dominated by weak or medium variations. Strong variations were mainly observed in high altitude regions, or in the combined section for encompassing plain and mountainous regions. For seasonal NSAT indices, it also mainly exhibited weak or medium variations, with different patterns in different seasons. Specifically, strong variations were identified in autumn and spring, followed by winter, while summer was the most stable, denoting that autumn maybe more susceptible to the warming.

Previous studies indicated that topographic information might be the key factor in climate indices prediction, especially in mountainous or high altitude areas with low density of meteorological stations (Gao et al., 2018; Peng et al., 2019). Therefore, it is possible that the precision of interpolation could be further improved by incorporating more detail terrain-related factors (Cheng et al., 2020). Nevertheless, traditional interpolation method is usually processed under the assumption that the NSAT is only dependent on the altitude. However, Šafanda (1999) demonstrated the strong dependence of NSAT on slope angle and orientation. Therefore, it is possible that the precision of interpolation for NSAT could be improved by considering slope angle and orientation in the interpolation process.

In this study, we included longitude and latitude as independent variables, and incorporated elevation, slope and aspect as covariates. We found that compared to the ANUSPLIN model, the MAE and RMSE values decreased by 0–5.77 and 1.96–8.86% for annual parameters, and 1.43–4.65 and 0–7.35% for seasonal parameters. As a result, the M-ANUSPLIN exhibited smaller deviation and performed better against other interpolation methods, especially in mountainous regions. However, our results also reflected the relatively poor performance of the M-ANUSPLIN model in regions with altitude >4,000 m, suggesting that it is not optimal for all regions.

Our study confirmed that slope is an important terrain-related factor in the interpolation process for NSAT (Table 4). As slope angle affects the temperature of surface objects by influencing the incidence angle and reflectivity of solar radiation, and then alters NSAT (Li et al., 2015; Peng et al., 2020). However, our study revealed the much smaller contribution of slope orientation to NSAT estimation than slope angle. This is not consistent with existing knowledge, which believes that slope orientation plays an important role in NSAT prediction. Earlier studies reported that in the middle latitudes of the Northern Hemisphere, the north slopes are generally colder at the same elevation than the south slopes because sunny aspects receive more direct solar radiation than northern aspects (Šafanda, 1999; Li et al., 2015). One reason for this inconsistent might be due to the different vegetation types in Southwest China. For example, Šafanda (1999) explained that the surface temperature in the meadow is higher than that in the forest because much of the Sun radiation is absorbed by the trees. Thus, the NSAT for meadows located at north slopes might be higher than the NSAT for forests located at south slopes. Another reason might be due to the existence and duration of the snow cover in high altitude regions, which can offset the effects of slope orientation on NSAT. Consequently, the effect of slope orientation on NSAT in Southwest China was smaller compared to those in other regions.

To examine the accuracy of M-ANUSPLIN interpolated dataset, the predicted results were compared to both the WorldClim 2.0 dataset (Fick and Hijmans, 2017) and the HMTC dataset (Peng et al., 2019). The 10-fold cross validation test was used to evaluate the overall error of the interpolated NSAT grid obtained from each dataset and the original observation data of 494 meteorological stations (Table 5).

Table 5 lists the mean, maximum, minimum values for annual and seasonal NSAT variables obtained from different datasets. It can be seen that the climatology anomaly of M-ANUSPLIN dataset is the lowest comparing to the WorldClim 2.0 and the HMTC datasets. Specifically, the anomalies are relatively high for HMTC (0.72–1.51°C for MAE and 1.19–1.97°C for RMSE); intermediate for WorldClim2 (0.66–0.87°C for MAE and 1.03–1.36°C for RMSE); and the lowest for M-ANUSPLIN (0.50–0.69°C for MAE and 0.82–1.14°C for RMSE). The seasonal NSAT variables shows the similar trend, where HMTC gives the highest estimation error, followed by WorldClim2 and M-ANUSPLIN. The good performance of M-ANUSPLIN can be mainly attributed to two reasons. Firstly, in case of WorldClim2 and HMTC datasets, only ∼300 sites and ∼500 sites of in-situ observation stations were respectively used for interpolation, and therefore significant biases occurred due to the complex terrain (Peng et al., 2019). Secondly, the M-ANUSPLIN model incorporated more detailed topographic information than traditional models, and thus can capture NSAT features with better precision.

Overall, the results indicate that the improved M-ANUSPLIN model can produce more accurate NSAT values than traditional interpolation models, and has apparent advantages over other interpolation methods in complex terrain areas like Southwest China.

Whilst this study demonstrated that the improved M-ANUSPLIN method can estimate more accurate NSAT than traditional models, especially in complex terrain areas like Southwest China, uncertainties still remain regarding the density of meteorological stations, input datasets and urbanization effect. Firstly, as gauge stations are relatively sparse in the western region of the study area, acquiring a correct distribution of NSAT through interpolation is difficult. The interpolated NSAT grid surface might contain some biases, which could introduce some uncertainty, especially in regions where the NSAT varies significantly in space and time. In addition, as the satellite SRTM product represents the average value in a 500 m × 500 m pixel, it does not provide the fine details and thus could not fully represent the real situation, which might also lead to additional biases.

Secondly, another source of uncertainty could be due to the input datasets. The NSAT is affected not only by topographical factors, but also closely related to other factors, e.g., vegetation and soil (Cho and Choi, 2014; Lensky et al., 2018). Moreover, as mentioned above, the effects of snow cover could be stronger than general expectation. It implies the need of careful interpretation of NSAT. Therefore, to further improve the accuracy of NSAT prediction, more information should be taken into account in the future.

Thirdly, the urbanization heat effect, which has a strong warming effect on NSAT (Kalnay and Cai, 2003; Ren et al., 2008; Luo and Lau, 2021), has not been considered in this study. This could cause an underestimation of NSAT, especially in the eastern region of Southwest China with a lot of big cities. Nonetheless, the results of this work should still, at least qualitatively, reveal the trends and spatial-temporal variations of NSAT in Southwest China over the past 50 years.

Lastly, elevation is an important factor affecting spatial variability of climate (Pepin et al., 2015). It is well known that for a stationary atmosphere, an increase in elevation leads to a subsequent decrease in air pressure and NSAT (You et al., 2008). According to EI Kenawy et al., 2009, different regions can experience different variations of NSAT, as each region has a unique terrain (Limsakul & Goes, 2008). Moreover, the variation of NSAT was not consistent between high and low land areas. Many studies suggested that NSAT increased more rapidly at higher than at lower elevations (Beniston, 2003; You et al., 2008; Pepin et al., 2015), which has been defined as elevation-dependent warming (EDW). However, this faster warming is not ubiquitous across the globe (Thakuri et al., 2019). In general, this study agrees with the previous studies of the EDW phenomenon with most high-altitude areas show significant warming. Nevertheless, we also observed a significant cooling trend in some high altitude regions like the northwestern part of Southwest China in the past 50 years. The reason for this phenomenon is still unclear, thus further investigation would be required.