The Prediction of the Tibetan Plateau Thermal Condition with Machine Learning and Shapley Additive Explanation

Abstract

The thermal condition over the Tibetan Plateau (TP) plays a vital role in the South Asian high (SAH) and the Asian summer monsoon (ASM); however, its prediction skill is still low. Here, two machine learning models are employed to address this problem. Expert knowledge and distance correlation are used to select the predictors from observational datasets. Both linear and nonlinear relationships are considered between the predictors and predictands. The predictors are utilized for training the machine learning models. The prediction skills of the machine learning models are higher than those of two state-of-the-art dynamic operational models and can explain 67% of the variance in the observations. Moreover, the SHapley Additive exPlanation method results indicate that the important predictors are mainly from the Southern Hemisphere, Eurasia, and western Pacific, and most show nonlinear relationships with the predictands. Our results can be applied to find potential climate teleconnections and improve the prediction of other climate signals.

1. Introduction

The Tibetan Plateau (TP) plays an important role in the climate system [1]. Existing studies have indicated the great impacts of the TP on the regional climate, especially the Asian summer monsoon (ASM) [2,3,4,5,6]. In boreal winter, the mechanical effects of the TP contribute to the Asian winter monsoon, while in boreal summer, the thermal effects of the TP influence the ASM [1,3,7]. The sensible heat air pump over the TP influences the formation, onset, and evolution of the ASM [5,6,8]. In boreal summer, the formation of the South Asian high (SAH), also known as the Tibetan Plateau high, is attributed to the thermal effects of the TP [1]. The tropospheric temperature over the TP largely reflects the variability in the SAH and the thermal conditions over the TP [9]. The boreal summer TP tropospheric temperature (TPTT) is utilized to represent the thermal condition of the TP [10,11]. The TPTT is closely related to the variabilities in the SAH and East Asian summer monsoon (EASM; Figure S1). The variance explained by the TPTT for the SAH and EASM is greater than 30%, indicating that the TPTT can well characterize the movement of the SAH and the EASM variability. Therefore, the TPTT index is the predictand of our study.

Despite the importance of the TP thermal condition, the precursors and prediction of TP heating still need further research. There are few studies on the prediction of TP thermal conditions, and the prediction results come from only dynamic models. Due to the high altitude and complex terrain, the simulated bias over the TP makes the prediction of its thermal condition more difficult [12]. Existing studies on the precursors of TP heating mainly focus on the tropical oceans, the North Atlantic Oscillation, and Arctic Sea ice [13,14,15,16,17,18]. Few studies have focused on the Southern Hemisphere precursors of TP heating [19,20]. Novel prediction methods need to be employed to improve the prediction of the TP thermal condition and analyze the vital predictors for the TP thermal condition.

Machine learning models have been employed in meteorology in recent years. Previous studies have confirmed that machine learning models are superior to traditional statistical methods and dynamic models in climate problems [21,22,23,24,25,26]. The application of machine learning in meteorology has been extended into climate mode/signal and precipitation prediction [21,25,26,27], the subgrid of climate models [24], the revision of model outputs [28], and climate analysis [22]. In contrast to deep learning methods, some machine learning models, such as gradient boosting tree methods, do not rely on large amounts of training data and computing resources. Medium or small datasets are sufficient to achieve good performance [23,25,28]. Therefore, it is possible to use limited, high-quality observational data in machine learning. When dealing with regression problems, the input of the decision tree of boosting tree models depends on the previous decision tree, which is updated cyclically through gradient descent to reduce the errors. Existing studies have confirmed that boosting tree models often outperform other machine learning models, such as random forest, in Geoscience [29,30,31,32,33,34]. To address the inexplicability of machine learning models, also known as the ‘black boxes,’ the SHapley Additive exPlanation method (SHAP) is employed to improve the interpretability of the machine learning models [35]. Depending on the additivity and comparability between different machine learning models, SHAP can analyze the contributions and the interactions among different predictors, which are also called features in machine learning [36,37].

We aim to solve the following problems: (1) Can machine learning models improve the prediction skills of the TP thermal condition (here referred to as TPTT)? If so, what is the improvement relative to the dynamic models? (2) Can interpretive methods be used to quantitatively characterize the importance of the features? If so, are these important features consistent with previous theoretical studies? Here, we concentrate on possible boreal spring precursors of TPTT because seasonal predictions can provide a reference for the operational forecasts. In addition, boreal spring is the season that most existing studies related to TP heating predictors have focused on [15,19,38,39].

2. Materials and Methods

2.1. Data

The Met Office Hadley Centre’s monthly sea ice concentration (SIC) dataset (v2.2.0.0) from 1850 to 2020 is utilized, which has a horizontal resolution of 0.25$°\times$0.25$°$[40]. The weekly snow cover extent (SC) data in the Northern Hemisphere, with a Cartesian grid from 1966 to 2021, are used [41]. The monthly Berkeley Earth land/ocean temperature (TS) record from 1850 to the present is employed with a horizontal resolution of 1$°\times$1$°$ [42]. We use the monthly air temperature from the European Centre for Medium-Range Weather Forecasts (ECWMF) fifth-generation reanalysis product (ERA5; 1950–present) with a resolution of 0.25$°\times$ 0.25$°$. The TPTT index is defined by the boreal summer (June–August; JJA) eddy air temperature over the TP averaged over 200–500 hPa (Figure 1a blue box; $25°–45°$N, $65°–105°$E) [10,11]. Here, the eddy air temperature is the air temperature minus its zonal mean.

The above datasets have different time coverages; their common dataset coverage is from 1967 to 2020. We use the dataset from 1967 to 2010 to train the machine learning models and the dataset from 2011 to 2020 to validate the performance of the different models (Figure 1b). We also analyze the hindcast results of the operational dynamic models ECMWF-SEAS5 and DWD-GCFS2.1, which are two state-of-the-art models from the Copernicus Climate Change Service seasonal multisystem. These hindcast results are compared with the machine learning prediction results over the 2011–2020 period. All variables are detrended to remove the global warming effect.

2.2. Prediction Models

The eXtreme Gradient Boosting model (XGBoost) is a classification and regression model based on the gradient boosting decision tree method [43]. The Light Gradient Boosting Machine (LightGBM) is also a tree-based gradient boosting method that can solve high-dimensional input variable problems [44]. Existing studies have confirmed the good prediction performance of these two models [25,33,45,46]. The two machine learning models consist of many simple weak learners (also known as the small regression models), and the final predictions are the weighted sum of the predictions of all weak learners. Moreover, as boosting tree models, XGBoost and LightGBM are not sensitive to multicollinearity because the response of boosting tree models to their features relies nonlinearly on the upper-level tree outputs, which reduces the simultaneous interaction of features [47].

A hyperparameter optimization framework for machine learning named optuna is employed [48]. Optuna is based on ‘Bayesian Optimization’ and is more efficient than other optimization methods, such as ‘grid search’ and ‘randomized search’. The hyperparameters of XGBoost and LightGBM are optimized by optuna within 500 iterations to obtain the optimal hyperparameters. A 5-fold cross-validation is used to prevent overfitting.

2.3. SHAP Method

Due to the ‘black box’ nature of machine learning models, the SHAP method is utilized to evaluate the concise causality relationships in the machine learning models (also known as improving the interpretability of the XGBoost and LightGBM) [35]. This method has been widely used in many fields [37,49,50]. The SHAP method has three desirable properties: local accuracy, missingness, and consistency. The consistency indicates that when the feature importance of a model increases or is maintained, the attribution value of this feature will increase or remain unchanged. Here, ‘model’ means the machine learning model. Missingness means that the feature attribution equals 0 when the feature is missing ($z_i^{\prime }=0$). The local accuracy means that the sum of all the feature attributions equals the local specific outputs of the model, shown as follows:

$$\mathcal{F}\left(x_i\right)\approx ℊ\left(z^{\prime }\right)={\varphi }_0+{\displaystyle \sum }_{i=1}^M{\varphi }_i{z}_i^{\prime }$$

(1)

where $ℊ\left(z^{\prime }\right)$ is the explanation model, also known as the output of the model $\mathcal{F}\left(x_i\right)$; ${\varphi }_0$ is the output of the dummy model with no features; ${\varphi }_i$ is the attribution of feature $i$, also known as the SHAP value; $z_i^{\prime }$ equals 1 when feature $i$ is used, and is otherwise 0; and M is the number of all features. The SHAP value is defined as follows:

$${\varphi }_i={\displaystyle \sum }_{S\subseteq N\backslash \left\{i\right\}}\frac{\left|S\right|!\left(M-\left|S\right|-1\right)!}{M!}\left[\mathcal{F}\left(S\cup \left\{i\right\}\right)-\mathcal{F}\left(S\right)\right]$$

(2)

where N is the set of all input features and S is the set of nonzero indices in $z^{\prime }$. SHAP consistently interprets the outputs of different machine learning models. The SHAP value can provide comparisons between different models. As with the SHAP method, the feature-related SHAP values are additive (local accuracy property). In addition, the Tree SHAP explainer we used to explain the tree models in this study can provide the SHAP value of the feature importance of a single feature and the SHAP interaction value to explain the interaction of different local features [36].

Different models provide different indicators to measure feature importance, and the indicators from different models cannot be compared and superimposed. Under different indicators, the feature importance ranking obtained by the same model is not the same, so a comparable feature importance indicator is needed. Due to the local accuracy of SHAP, the SHAP value has the same physical meaning as the predictand, and the SHAP values of different models are also comparable and additive. The feature importance values directly provided by XGBoost and LightGBM have no physical meaning. To eliminate the difference in feature importance caused by the selection of different indicators and obtain comparable and additive feature importance, we select the SHAP value as the measurement indicator.

2.4. Evaluation Metrics

The distance correlation $DCORR\left(x,y\right)$ and root-mean-square error (RMSE) are used to evaluate the prediction skills, which are defined as follows:

$$\begin{array}{c}DCORR\left(x,y\right)=\frac{dcov\left(x,y\right)}{\sqrt{dcov\left(x,x\right)dcov\left(y,y\right)}} \end{array}$$

(3)

$$dco{v}^2\left(x,y\right)=\widehat{S_1}+\widehat{S_2}-2\widehat{S_3}$$

(4)

$$\begin{array}{c}\widehat{S_1}=\frac{1}{n^2}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{\parallel x_i-x_j\parallel }_{dx}{\parallel y_i-y_j\parallel }_{dy} \end{array}$$

(5)

$$\begin{array}{c}\widehat{S_2}=\frac{1}{n^2}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{\parallel x_i-x_j\parallel }_{dx}\frac{1}{n^2}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{\parallel y_i-y_j\parallel }_{dy} \end{array}$$

(6)

$$\begin{array}{c}\widehat{S_3}=\frac{1}{n^2}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{\displaystyle {\displaystyle \sum }_{l=1}^n}{\parallel x_i-x_l\parallel }_{dx}{\parallel y_i-y_l\parallel }_{dy} \end{array}$$

(7)

$$\begin{array}{c}\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-y_i\right)}^2}{n}} \end{array}$$

(8)

where $n$ is the length of $x$ and $y$ [51]. Compared with the Pearson correlation, the distance correlation can evaluate linear and nonlinear relationships. If the distance correlation between $x$ and $y$ equals 0, the two variables are independent of each other. Previous studies have verified the good performance of distance correlation in feature screening [52]. In addition, the Monte Carlo significance test is employed to test the significance of the distance correlation.

3. Results

Figure 1a shows a boreal summer climatological tropospheric eddy temperature anomaly at 200–500 hPa over the Asian continent, with the high-temperature center extending from North Africa to the western North Pacific. The peak temperature can reach 5 $°K$ and is located over the southern TP, with a large warm center found by Li and Yanai [53]. Existing studies indicate that the warming rate of the TP is higher than the average rate in the Northern Hemisphere and almost twice the global mean [54]. Here, the blue box contains the whole TP, which is used to represent the thermal condition over the TP [10,11]. The standardized detrended TPTT sequence is shown in Figure 1b. The value range of the test dataset is from −1.4 to 2.0, which is within the range values of the training dataset (from −2.4 to 2.0). This also partly shows that the training dataset and the test dataset have the same distribution, so it is reasonable for us to divide the training data and test data in this way.

Because XGBoost and LightGBM cannot extract features automatically, we need to select features to train the machine learning models using our expert knowledge. We focus on the seasonal timescale prediction, so we mainly consider the variables with the storage capability of the signals. Considering previous studies on TP heating, the SIC, SC, and TS are selected as the key variables [13,16,17,18,39]. Figure 2 shows the DCORR between the TPTT and the previous boreal spring variables. The significantly correlated regions of the SIC are mainly located around the Antarctic, and the highly correlated regions of SC are over the eastern TP and the Eurasian continent (Figure 2a,b). Compared with the SIC and SC, the TS is closely related to the TPTT at low, middle, and high latitudes (Figure 2c). To ensure the close relationships between these features and the TPTT, we define the regional averaged value in these correlated regions that pass the 99% Monte Carlo test as the features (also known as the predictors). The 15 selected features all show significant connections with the TPTT index (Figure S2). For convenience, the 15 features are abbreviated with variable names and locations in Figure 2. Qian, Jia, Lin, and Zhang [25] also used a similar method for screening features and found it feasible. The 15 features are utilized for training the XGBoost and LightGBM.

Note that significant boreal spring signals corresponding to TPTT exist in two polar areas and middle-to-low latitude areas, especially in the Southern Hemisphere (Figure 2a,c). However, existing studies pay little attention to the Southern Hemisphere precursors of TP heating. In the following, to clarify the importance of each precursor for the predictands, the outputs of the machine learning models are analyzed to evaluate the models’ prediction skills and quantify the importance of the features for the TPTT.

The observed and predicted results of TPTT are shown in Figure 3a. The DCORR between the observation and the prediction of both XGBoost and LightGBM can pass the 95% significance test, and the prediction skills of LightGBM (DCORR = 0.77; RMSE = 0.79) are better than these of the XGBoost (DCORR = 0.7; RMSE = 0.79). The ensemble of both machine learning results can be closer to the observed results (DCORR = 0.82; RMSE = 0.71), which can pass the 99% Monte Carlo test. The hindcasts of the two state-of-the-art models, ECMWF-SEAS5 and DWD-GCFS2.1, are chosen for comparison with the machine learning results, and these two dynamic models show excellent prediction skills on seasonal timescales [55,56]. The DCORRs (RMSEs) between the dynamic model results and the observations are 0.48 (1.13) in DWD-GCFS2.1 and 0.54 (0.98) in ECMWF-SEAS5, which are lower than the machine learning results and cannot pass the 90% significance test. The sign test (sign consistency between the predictions and the observations) of the first-order difference results also indicates that the results of XGBoost and LightGBM are perfectly consistent with the observations. However, only 80% (ECMWF-SEAS5) and 60% (DWD-GCFS2.1) of the predictions are consistent with the observations in these two dynamic models. Therefore, this further illustrates the superiority of the machine learning models over dynamic models for TPTT prediction.

Figure 3a shows that the predictions of the machine learning models are highly consistent with the observations, even in years with an extremely high TPTT index, such as 2018, and the predicted value is very close to the observed value in XGBoost. The extreme year of 2018 is analyzed to understand the prediction model. The decision routine (also known as the decision plot) of the XGBoost model for 2018 is shown in Figure 3b. The trend of the curve represents the accumulation of SHAP values of all features. For example, the value of ts_wd is −0.039 in 2018, and the largest positive SHAP value is in the feature ts_wd, which means that ts_wd contributes greatly to increasing the predicted TPTT value in 2018. In contrast, ts_ac attributes to the decreasing prediction in 2018. There are other features that greatly contribute to the 2018 TPTT value, such as ts_ea, ts_aa, ts_nz, and sc_ea, and these features are mainly from the Southern Hemisphere.

The SHAP analysis of a single sample is not enough to fully understand the logic of the overall prediction of the machine learning models. Therefore, we calculate the global and local feature importance of all the samples. The XGBoost and LightGBM ensemble results of feature importance are shown in Figure 4. The global feature importance in Figure 4a shows the absolute SHAP value average in all the samples, which represents the feature attributions in all the samples, and the ranking importance results are like those of the single sample analysis (Figure 3b). The local feature importance in Figure 4b indicates the distributions of the feature attribution in all the samples. As space is limited in this paper, the most important six features for which the SHAP values are larger than 0.1 are analyzed (Figure 4a). The dependence plots of the six features are shown in Figure 5a–f and are detailed and displayed in Figure 4b. The dependence plot shows the relationship between the feature values and the related SHAP values and can help clarify the influence of the features on the predictands. The large and small feature values of ts_wd are distributed on both sides in Figure 4b. The increase in ts_wd increases the SHAP value of ts_wd, which increases the predicted TPTT index. However, when ts_wd is larger than 0.5, its influence on TPTT weakens, but it still attributes to a larger prediction, showing a nonlinear relationship (Figure 5a). The effect of ts_nz is like that of ts_wd, except that it is a turning point of 0.2 (Figure 5e). The positive feature values of ts_ea and ts_aa correspond to a decrease in the TPTT index and vice versa. In addition, the distribution of these two features exhibits the characteristics of a sign function (Figure 5b,c). As sc_ea increases, the TPTT prediction decreases, but its attribution tends to remain stable after the feature values tend to be positive (Figure 5d). The influence of ts_wp has a long right-tailed distribution, and ts_wp feature values less than −0.5 have great attributions to the increasing TPTT prediction and, there is an obvious linear relationship between the SHAP values of ts_wp and its feature values (Figure 4b and Figure 5f). The above lists the influence of a single feature on the TPTT index. The interaction between two features also affects the TPTT prediction; however, the interaction effects are relatively weak compared to single features (Figure S3). According to the interaction dependence plot, the interaction between ts_nz and ts_wp is chosen. When ts_nz is less than −0.3, and ts_wp is greater than 0, the attribution of this interaction leads to an increase in the TPTT index (Figure 5g).

Most of the important features ranked by the machine learning models and SHAP are in the Southern Hemisphere. The TS anomaly in the eastern Antarctic Peninsula is closely related to the Antarctic Oscillation (AAO), and the in-phase relationship between ts_wd and TPTT further confirms the results of existing studies [19,57]. Ding et al. [58] points out that there is an internal sea surface temperature (SST) quadruple mode in the South Pacific that is driven by the Pacific South American mode, which is similar to the results shown in Figure 2c. Liu et al. [59] and Hsu and Chen [60] indicate that SST anomaly in the South Pacific can influence the western North Pacific and East Asia, which explains the attributions of ts_ea and ts_nz to the subsequent TPTT. Vernekar, Zhou, and Shukla [39] utilize composite analysis to suggest that Eurasian snow cover can decrease the air temperature over the TP, and their results are consistent with our analysis. The Philippine Sea anticyclone teleconnection of the El Niño–South Oscillation leads to a positive SST anomaly in the western North Pacific (ts_wp), which can also induce decreased TPTT [14,61]. In addition, ts_wp may also be linked to the Victoria mode and Pacific Decadal Oscillation and has a vital impact on the SAH and the TPTT [38,62].

4. Discussion

4.1. Uncertainty

SIC data is a product of the fusion of observation data and satellite data, TS data is station-based data, and ERA5 data is a product of the assimilation of observation data, satellite data, and other data. The SC data is pure satellite-based data with uncertainty, so we mainly analyze the uncertainty of SC data. Previous uncertainty assessments of SC suggest that the bias due to illumination, cloud cover, and confounding of SC and other surfaces can lead to overestimation of boreal spring and summer SC in the Arctic and northern Canada [63,64,65]. Although SC data overestimates Arctic SC in June, the SC in boreal spring is consistent with other independent data sets, so our analysis using Eurasian SC data in boreal spring is reliable [66].

4.2. Comparison, Limitations and Contributions

In this paper, the limited observational data are used for seasonal-scale climate prediction and high prediction skills are achieved through two boosting tree models. The random forest model is also utilized. Considering the prediction skill of random forest is poor compared to the boosting tree models, it is not shown in this paper. In meteorology, deep learning and relatively simple machine learning models have been applied to predictions. Due to the powerful fitting ability of deep learning, deep learning models can usually obtain high prediction skills [21,67]. Owing to the large amount of data demanded by deep learning, predictions using numerical simulation data and deep learning need to be further developed, which may further improve prediction skills.

In terms of feature selection, similar feature selection methods (Pearson correlation coefficient screening) and boosting tree models are employed by other researchers for seasonal prediction [25]. Since the Pearson correlation cannot measure the nonlinear relationship, it may miss the features with high nonlinear correlation, so the distance correlation is selected to overcome this problem [52]. The other feature selection methods include the empirical orthogonal function and maximum variance unfolding, etc. [68,69]. These methods are advantageous for the screening of local features. For the selection of global predictors, these methods have obvious limitations. In addition, we also screened the features according to the global climatic annual precipitation, but the correlations between the selected features and the predictand are extremely poor, and the prediction skills were also very low (figures not shown).

In this research, due to the temporal limitation of the data, we use the data from 1967 to 2010 to train the machine learning models, and thus the models only learn the features from the 44 years of data. For signals that do not appear in these 44 years, the prediction skills of the machine learning models are consequentially poor. Therefore, longer time scale datasets should be considered in future studies to improve the generalization ability of the machine learning models. In this paper, we focus on the surface climate variables (SC, SIC, and TS) and do not consider factors of stratospheric signals and other non-meteorological elements, which can also be further studied in the future [70,71].

Previous studies indicate that the Indian Ocean SST has important effects on the SAH and TPTT [13,72,73,74,75]. However, the impact of the Indian Ocean (ts_io) is underestimated in the machine learning models, and similar results also appear in the deep learning method related to the EASM [26]. This underestimation problem is still unresolved in machine learning, so we need to be careful to avoid the underestimation of the Indian Ocean in future studies. Nonlinear relationships in the climate system are often ignored in quantitative studies. However, the vital impacts of these nonlinear interactions should be given sufficient attention (Figure 5). By comparing our results with previous studies, we find that some features, such as ts_aa, show a very important predictive significance in machine learning models, but there is no corresponding research on this subject. Therefore, further research is needed to verify whether such a physical connection exists.

In this paper, the data used are heterogeneous, such as SC, SIC, and TS, with different physical meanings. SHAP is used to transform heterogeneous features with different physical meanings into uniform ones under the nature of additivity. The observed datasets from 1967 to 2010 are used to train the statistical machine learning models; however, the results still need to be examined in terms of climate dynamics theories and numerical simulations. Therefore, the analysis results obtained by SHAP are a quasi-explanation, not the final explanation, but it is nevertheless meaningful for revealing objective phenomena (connections) that are not yet discovered, which helps to fill the gap in the interpretability of machine learning models of the TPTT. In addition, there are few studies on the prediction of TP heat source and its interpretability, our prediction based on machine learning and the SHAP method may provide ideas for the operational prediction and precursor studies.

5. Conclusions

This paper employs two machine learning models (XGBoost and LightGBM) and the SHAP method to predict the TPTT index and analyze its key precursors. Expert knowledge and distance correlation are applied to extract highly correlated features, which are supported by theoretical knowledge. Both linear and nonlinear relationships are considered in choosing the features. The features are fed into the machine learning models to obtain a final prediction. The prediction skills of the two machine learning models can pass the 99% significance test (DCORR = 0.82), with an RMSE of 0.71. The prediction skills of both state-of-the-art dynamic models (ECMWF-SEAS5 and DWD-GCFS2.1) are lower than the prediction skills of the machine learning models, which indicates that machine learning models can significantly improve the prediction of TPTT. By using the SHAP method, we find that the importance of the features in all the samples is ranked from high to low as follows: ts_wd, ts_ea, ts_aa, sc_ea, ts_nz, ts_wp, ts_am, sc_tp, ts_ac, ts_np, ts_af, ts_sp, ts_ep, sic_s, and ts_io. The six most important features that attribute the most to the TPTT prediction are mainly located in the Southern Hemisphere, the western North Pacific, and Eurasia, and most of these features have nonlinear relationships with the TPTT. The influences of these most important features on the TPTT are generally consistent with the existing studies.