Developing a novel hybrid model for the estimation of surface 8 h ozone (O₃) across the remote Tibetan Plateau during 2005–2018

We developed a two-stage model called the random-forest–generalised additive model (RF–GAM), based on satellite data, meteorological factors, and other geographical covariates, to predict the surface 8 h O₃ concentrations across the remote Tibetan Plateau. The 10-fold cross-validation result suggested that RF–GAM showed excellent performance, with the highest R² value (0.76) and lowest root-mean-square error (RMSE) (14.41 µg m⁻³), compared with other seven machine-learning models. The predictive performance of RF–GAM showed significant seasonal discrepancy, with the highest R² value observed in summer (0.74), followed by winter (0.69) and autumn (0.67), and the lowest one in spring (0.64). Additionally, the unlearning ground-observed O₃ data collected from open-access websites were applied to test the transferring ability of the novel model and confirmed that the model was robust in predicting the surface 8 h O₃ concentration during other periods (R²=0.67, RMSE = 25.68 µg m⁻³). RF–GAM was then used to predict the daily 8 h O₃ level over the Tibetan Plateau during 2005–2018 for the first time. It was found that the estimated O₃ concentration displayed a slow increase, from 64.74±8.30 µg m⁻³ to 66.45±8.67 µg m⁻³ from 2005 to 2015, whereas it decreased from the peak to 65.87±8.52 µg m⁻³ during 2015–2018. Besides this, the estimated 8 h O₃ concentrations exhibited notable spatial variation, with the highest values in some cities of the northern Tibetan Plateau, such as Huangnan (73.48±4.53 µg m⁻³) and Hainan (72.24±5.34 µg m⁻³), followed by the cities in the central region, including Lhasa (65.99±7.24 µg m⁻³) and Shigatse (65.15±6.14 µg m⁻³), and the lowest O₃ concentration occurred in a city of the southeastern Tibetan Plateau called Aba (55.17±12.77 µg m⁻³). Based on the 8 h O₃ critical value (100 µg m⁻³) provided by the World Health Organization (WHO), we further estimated the annual mean nonattainment days over the Tibetan Plateau. It should be noted that most of the cities on the Tibetan Plateau had excellent air quality, while several cities (e.g. Huangnan, Haidong, and Guoluo) still suffered from more than 40 nonattainment days each year, which should be given more attention in order to alleviate local O₃ pollution. The results shown herein confirm that the novel hybrid model improves the prediction accuracy and can be applied to assess the potential health risk, particularly in remote regions with few monitoring sites.

1 Introduction

Along with the rapid economic development and urbanisation, the anthropogenic emissions of nitrogen oxides (NO_x) and volatile organic compounds (VOCs) displayed high-speed growth. The chemical reactions between NO_x and VOCs in the presence of sunlight triggered ambient ozone (O₃) formation (Wang et al., 2017, 2019). As a strong oxidant, ambient O₃ could play a negative role in human health through aggravating the cardiovascular and respiratory function (Ghude et al., 2016; Marco, 2017; Yin et al., 2017a). Apart from the effect on human health, O₃ also posed a great threat to vegetation growth (Emberson et al., 2018; Feng et al., 2015, 2019; Qian et al., 2018). Moreover, the tropospheric O₃ can perturb the radiative energy budget of the earth–atmosphere system, as it is the third most important greenhouse gas next to carbon dioxide (CO₂) and methane (CH₄), thereby changing the global climate (Bornman et al., 2019; Fu et al., 2019; Wang et al., 2019). Recently, the particulate matter with a concentration of less than 2.5 µm (PM_2.5) showed a persistent decrease, while the O₃ issue has been increasingly prominent in China (Li et al., 2017b, 2019b). Therefore, it was critical to accurately reveal the spatiotemporal variation in O₃ pollution and assess its heath risk in China.

A growing body of studies began to investigate the spatiotemporal variation in the O₃ level worldwide. Wang et al. (2014b) demonstrated that the 8 h O₃ concentrations in nearly all of the provincial cities experienced remarkable increases during 2013–2014. Following this work, Li et al. (2017b) reported that the annual mean O₃ concentration over China increased by 9.18 % during 2014–2016. In other Asian countries except China, Vellingiri et al. (2015) performed long-term observation and found that the O₃ concentration in Seoul, South Korea, has displayed a gradual increase in recent decades. In the southeastern United States, Li et al. (2018) observed that the surface O₃ concentration has displayed a gradual decrease in the past 10 years. Although the number of ground-level monitoring sites has been increasing globally, the limited monitoring sites still cannot accurately reflect the fine-scale O₃ pollution status because each site shows little spatial representativeness (0.25–16.25 km²) (Shi et al., 2018). Furthermore, the number of monitoring sites in many countries (e.g. China and the United States) displays an uneven distribution characteristic at the spatial scale. In China, most of these sites are concentrated in the North China Plain (NCP) and Yangtze River Delta (YRD), while western China has an extreme lack of ground-level O₃ data, which often increases the uncertainty of health assessment. Therefore, many studies used various models to estimate the O₃ concentrations without monitoring sites. Chemical transport models (CTMs) were often considered to be the typical methods to predict the surface O₃ level. Zhang et al. (2011) employed the GEOS-Chem model to simulate the surface O₃ concentration over the United States, suggesting that the model could capture the spatiotemporal variation in surface O₃ concentration at a large spatial scale. Later on, Wang et al. (2016) developed a hybrid model called land use regression (LUR) coupled with CTMs to predict the surface O₃ concentration in the Los Angeles Basin, California. In recent years, these methods were also applied to estimate the surface O₃ level over China. Liu et al. (2018) used the Community Multiscale Air Quality (CMAQ) model to simulate the nationwide O₃ concentration over China in 2015. Nonetheless, the high-resolution O₃ prediction using CTMs might have widely deviated from the measured value, owing to the imperfect knowledge about the chemical mechanism and the higher uncertainty of the emission inventory. Moreover, the continuous emission data of NO_x and VOCs were not always open access, which restricted the long-term estimation of the surface O₃ concentration using CTMs.

Fortunately, the daily satellite data enable the fine-scale estimations of the O₃ level at a regional scale due to broad spatial coverage and high temporal resolution (McPeters et al., 2015). Shen et al. (2019) confirmed that the satellite-retrieved O₃ column amount could accurately reflect the spatiotemporal distribution of the surface O₃ level. Therefore, some studies tried to use traditional statistical models coupled with high-resolution satellite data to estimate the ambient O₃ level. Fioletov et al. (2002) used the satellite measurement to investigate the global distribution of O₃ concentrations based on a simple linear model. Recently, Kim et al. (2018) employed the integrated empirical geographic regression method to predict the long-term (1979–2015) variation in ambient O₃ concentration over the United States based on O₃ column amount data. Although the statistical modelling of ambient O₃ concentration is widespread around the world, most of the traditional statistical modelling only utilised the linear model to predict the ambient O₃ concentration, which generally decreased the prediction performance because the nonlinearity and high-order interactions between O₃ and predictors cannot be managed by a simple linear model.

As an extension of traditional statistical model, machine-learning methods have been widely applied to estimate the pollutant levels in recent years because of their excellent predictive performances. Among these machine-learning algorithms, decision tree models such as random forest (RF) and extreme gradient boosting (XGBoost) generally showed fast training speed and excellent prediction accuracy (Li et al., 2020; Zhan et al., 2018). Furthermore, decision tree models can obtain the contribution of each predictor to air pollutants, which was beneficial to the parameter adaption and model optimisation. Chen et al. (2018b) has firstly employed the RF model to simulate the PM_2.5 level in China since 2005. Following this work, we recently used the XGBoost model to estimate the 8 h O₃ concentration on the island of Hainan for the first time and captured the moderate predictive performance (R²=0.59) (Li et al., 2020). While the decision tree model showed many advantages in predicting the pollutant level, the spatiotemporal autocorrelation of pollutant concentration was not a concern in these studies. Li et al. (2019a) confirmed that the prediction error by the decision tree model varied greatly with space and time. Thus, it is imperative to incorporate the spatiotemporal variables into the original model to further improve the performance. To resolve the defects of decision tree models, Zhan et al. (2018) developed a hybrid model called RF-spatiotemporal Kriging (STK) to predict the O₃ concentration over China and achieved better performance (overall – R²=0.69; southwestern China – R²=0.66). Unfortunately, the RF–STK model still showed some weaknesses in predicting O₃ concentration. First of all, the predictive performance of the STK model was strongly dependent on the number of monitoring sites and their spatial densities. The model often showed worse predictive performance in regions with few monitoring sites (Gao et al., 2016). Moreover, the ensemble model cannot simulate the O₃ level during the periods without ground-level-measured data. In contrast, the generalised additive model (GAM) not only considers the time autocorrelation of O₃ concentration but also shows better extrapolation ability (Chen et al., 2018a; Ma et al., 2015). Thus, the ensemble model of RF and GAM was proposed to predict the spatiotemporal variation in the surface 8 h O₃ concentration.

The Tibetan Plateau, the highest plateau in the world, shows higher surface solar radiation compared with the regions outside the plateau. It was well documented that high solar radiation tended to generate a large amount of the OH radical, resulting in the O₃ formation via the reaction of VOC and the OH radical (Ou et al., 2015). While the total O₃ column amount on the Tibetan Plateau has displayed a slight decrease since the 1990s, the convergent airflow formed by subtropical anticyclones could bring ozone-rich air surrounding the plateau to the low atmosphere (Lin et al., 2008), thereby leading to a higher surface O₃ concentration over the plateau. Most studies focused on the stratosphere–troposphere transport of O₃ on the Tibetan Plateau, whereas limited effort was given to investigating the ground-level O₃ level over this region. To date, only several studies were concerned with the spatiotemporal variation in the surface O₃ concentration in this region based on field-observation data (Chen et al., 2019; Shen et al., 2014; Yin et al., 2017b). Unfortunately, the few monitoring sites on the Tibetan Plateau cannot capture the real O₃ pollution status, especially in remote areas (e.g. the northern part of the Tibetan Plateau), because each site only possessed limited spatial representativeness. Apart from these field measurements, Liu et al. (2018) (R=0.60) and Zhan et al. (2018) (R²=0.66) used CTMs and the machine-learning model to simulate the surface O₃ concentration over China in 2015, respectively. Both of these studies included the predicted O₃ level on the Tibetan Plateau. Although they finished the pioneering work, the predictive performances of both studies were not excellent. Therefore, it was imperative to develop a higher-quality model to enhance the modelling accuracy.

Here, we developed a new hybrid-method (RF–GAM) model integrating satellite data, meteorological factors, and geographical variables to simulate the gridded 8 h O₃ concentrations over the Tibetan Plateau for the first time. Based on the estimated surface O₃ concentration, we clarified the long-term variation (2005–2018) of the surface O₃ concentration and quantified the key factors for the annual trend. Filling the gap of statistical estimation of the 8 h O₃ level in a remote region, this study provides useful datasets for epidemiological studies and air quality management.

2 Materials and methods

2.1 Study area

The Tibetan Plateau is located in southwestern China, which ranges from 26.00 to 39.58^∘ N and from 73.33 to 104.78^∘ E. The Tibetan Plateau is surrounded by the Taklamakan Desert to the north and Sichuan Basin to the southeast. The land area of the Tibetan Plateau reaches 2.50×10⁶ km² (Chan et al., 2006). Based on the air circulation pattern, the Tibetan Plateau can be roughly classified into the monsoon-influenced region and the westerly-wind-influenced region (Wang et al., 2014a). The annual mean air temperature in most regions is below 0 ^∘C. The annual mean rainfall amount on the Tibetan Plateau ranges from 50 to 2000 mm. The terrain conditions are complex, and higher altitudes are concentrated in the central region. The Tibetan Plateau is generally treated as a remote region lacking in anthropogenic activity, and most of the residents are concentrated in the southeastern and southern parts of the Tibetan Plateau. The Tibetan Plateau consists of 19 prefecture-level cities, and their names and corresponding geographical locations are shown in Figs. 1 and S1.

Figure 1 The geographical locations and annual mean 8 h O₃ concentrations in the ground-observed sites (red dots) over the Tibetan Plateau during 2014–2018. The elevation data are collected from geographical and spatial data cloud at a 30 m spatial resolution.

2.2 Data preparation

2.2.1 Ground-level 8 h O₃ concentration

The daily 8 h O₃ data in 37 monitoring sites over the Tibetan Plateau from 13 May 2014 to 31 December 2018 were collected from the national air quality monitoring network. The O₃ levels in all of these sites were determined using an ultraviolet-spectrophotometry method. The highest 8 h moving average O₃ concentration each day was calculated as the daily 8 h O₃ level after data quality assurance. The data quality of all the monitoring sites was assured on the basis of the HJ 630-2011 specifications. The data with no more than two consecutive hourly measurements missing in all the days were treated as the valid data.

2.2.2 Satellite-retrieved O₃ column amount

The O₃ column amounts (DU: total molecules cm⁻²) during 2005–2018 were downloaded from the Ozone Monitoring Instrument O₃ (OMI O₃) level-3 data with a 0.25^∘ spatial resolution from the website of the National Aeronautics and Space Administration (NASA) (https://acdisc.gsfc.nasa.gov/data/Aura_OMI_Level3/OMDOAO3e.003/, last access: 19 May 2020). The OMI O₃ product shows global coverage and traverses the earth once a day. The O₃ column amount with a cloud radiance fraction >0.5, terrain reflectivity >30 %, and solar zenith angles >85^∘ should be removed. In addition, the cross-track pixels significantly influenced by the row anomaly should be deleted.

2.2.3 Meteorological data and geographical covariates

The daily meteorological data were obtained from ERA-Interim datasets with 0.125^∘ resolution. These meteorological data consisted of the 2 m dew-point temperature (d2m), 2 m temperature (t2m), 10 m U wind component (u10), 10 m V wind component (v10), boundary layer height (blh), sunshine duration (sund), surface pressure (sp), and total precipitation (tp). The 30 m resolution elevation data (DEM) were downloaded from the China Resource and Environmental Science Data Center (CRESDC). The data of the gross domestic product (GDP) and population density with 1 km resolution were also extracted from CRESDC. Population density and GDP in 2005, 2010, and 2015 were integrated into the model to predict the surface 8 h O₃ concentration over the Tibetan Plateau because these data were available every 5 years. Additionally, the land use data of 30 m resolution (e.g. water, grassland, urban, forest) were also extracted from CRESDC. Lastly, the latitude, longitude, and time were also incorporated into the model.

All of the explanatory variables collected were resampled to 0.25^∘ × 0.25^∘ grids to predict the O₃ level. The original meteorological data with 0.125^∘ resolution were resampled to the 0.25^∘ grid. The land use area, elevation, GDP, and population density in each grid were calculated using spatial clipping. Lastly, all of the predictors were integrated into an intact table to train the model.

2.3 Model development and assessment

The RF–GAM model was regarded as the hybrid model of RF and GAM. The RF–GAM model was a two-stage model in which the prediction error estimated by the RF model was then simulated by GAM. The prediction results of RF and GAM were summed as the final result of the RF–GAM model (Fig. 2). The detailed equation is as follows:

$$\begin{array}{}\text{(1)}& Z(s,t)=P(s,t)+E(s,t),\end{array}$$

where Z(s, t) is the estimated 8 h O₃ level at the location s and time t, P(s, t) represents the 8 h O₃ concentration predicted by the RF model, and E(s, t) denotes the prediction error by GAM.

Figure 2 The workflow for predicting the spatiotemporal distributions of 8 h O₃ levels.

Figure 3 Density scatterplots of model fitting and cross-validation results at a daily level. Panels (a), (b), and (c) represent RF–GAM, RF–STK, and RF models, respectively. The red dotted line denotes the fitting linear-regression line. MPE, RMSE, and RPE are mean prediction error (µg m⁻³), root-mean-square error (µg m⁻³), and relative percentage error (%), respectively.

In the RF model, a large number of decision trees were planted based on the bootstrap sampling method. At each node of the decision tree, the random samples of all predictors were applied to determine the best split among them. Following the procedure, a simple majority vote was employed to predict the 8 h O₃ level. The RF model avoided a priori linear assumption of O₃ concentration and predictors, which was often not in good agreement with the actual state. The RF model has two key parameters, including n_tree (the number of trees grown) and m_try (the number of explanatory variables sampled for splitting at each node). The prediction performance of the RF model was strongly dependent on the two parameters. The optimal n_tree and m_try were determined based on the least out-of-bag (OOB) errors. Based on the iteration result, the optimal n_tree and m_try reached 500 and 5, respectively. Besides this, the backward variable selection method was performed on the RF submodel to achieve better performance. At each step of the predictor selection, the variable with the least important value was excluded from the next step. This one-variable-at-a-time exclusion method was repeated until only two explanatory variables remained in the submodel. Finally, all of the selected variables except the area of water were integrated into the model to achieve the best prediction performance. The detailed RF model is as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathrm{O}}_{\mathrm{3}}& ={\mathrm{O}}_{\mathrm{3}}\mathrm{column}+\mathrm{Elevation}+\mathrm{Agr}+\mathrm{Urban}+\mathrm{Forest}\\ & +\mathrm{GDP}+\mathrm{Grassland}+\mathrm{Population}+\mathrm{Prec}+T+\mathrm{WS}\\ & +P+t_{\mathrm{sun}}+\mathrm{RH},\end{array}\end{array}$$

where O₃ denotes the observed 8 h O₃ level in the monitoring site; the O₃ column represents the O₃ column amount in the corresponding grid; Elevation denotes the corresponding elevation of the site; and Agr, Urban, Forest, and Grassland are the agricultural land, urban land, forest land, and the grassland, respectively. Population represents the population density in the corresponding site. Prec, T, WS, P, t_sun, and RH are precipitation, air temperature, wind speed, air pressure, sunshine duration, and relative humidity, respectively. Additionally, another five models, including the RF, generalised regression neutral network (GRNN), backward-propagation neural network (BPNN), Elman neural network (ElmanNN), and extreme learning machine (ELM), also used the backward variable selection method. The R² value was treated as an important parameter for adding or reducing the variable. The variable should be removed when the R² value of the submodel showed a remarkable decrease with the integration of this variable. Lastly, the optimal variable group was applied to establish the submodel.

Following the RF submodel, the prediction error estimated by the RF submodel was further modelled by the GAM. GAM could reflect the time autocorrelation of the predictive error of RF model, and thus the ensemble model of RF and GAM might decrease the modelling error of the one-stage model. All of the variables were incorporated into the models to establish the second-stage model, and the backward variable selection was also used to determine the optimal variable group.

The 10-fold cross-validation (CV) technique was employed to evaluate the predictive performances for all of the machine-learning models. All of the training datasets were randomly classified into 10 subsets uniformly. In each round of validation, nine subsets were used to train, and the remaining subset was applied to test the model performance. The process was repeated 10 times until every subset has been tested. Some statistical indicators, including the R², root-mean-square error (RMSE), mean prediction error (MPE), relative percentage error (RPE), and the slope, were calculated to assess the model performance. The optimal model with the best performance was used to estimate the 8 h O₃ concentration in recent decades.

Table 1 The R² values, RMSE, MPE, and RPE of RF–GAM in different years during 2014–2018 over the Tibetan Plateau.

Figure 4 The transferring-ability validation of RF–GAM method based on the measured daily 8 h O₃ concentration during December 2013–May 2014.

Figure 5 The variable importance of predictors in the final RF–GAM model.

3 Results and discussion

3.1 The validation of model performance

Figures 3 and S2 show the density scatterplots of the fitting and 10-fold cross-validation results for eight machine-learning models for China. The 10-fold cross-validation R² values followed the order of RF–GAM (R²=0.76)>RF–STK (R²=0.63)>RF (R²=0.55)>GRNN (R²=0.53)>BPNN (R²=0.50)>XGBoost (R²=0.48)>ElmanNN (R²=0.47)>ELM (R²=0.32). The RMSE values of RF–GAM, RF–STK, RF, GRNN, XGBoost, BPNN, ElmanNN, and ELM were 14.41, 17.79, 19.13, 19.41, 20.73, 20.06, 20.61, and 23.36 µg m⁻³, respectively. Both the MPE and RPE showed similar characteristics to RMSE of the order of RF–GAM (10.97 µg m⁻³ and 26.50 %) < RF–STK (13.48 µg m⁻³ and 35.15 %) < RF (14.71 µg m⁻³ and 35.51 %) < GRNN (14.89 µg m⁻³ and 35.82 %) < BPNN (15.43 µg m⁻³ and 36.19 %) < ElmanNN (15.75 µg m⁻³ and 37.05 %) < XGBoost (15.80 µg m⁻³ and 38.13 %) < ELM (18.23 µg m⁻³ and 44.05 %) (Figs. 3 and S2). Besides this, the slope of the RF–GAM model was closer to 1 compared with other models. It was well documented that the RF model generally showed better performance than other models because this method did not need to define complex relationships between the explanatory variables and the O₃ concentration (e.g. linear or nonlinear). Furthermore, the variable importance indicators calculated by the RF model can help the user to distinguish the key variables from the noise ones and make full use of the strength of each predictor to assure the model robustness. Although BPNN, GRNN, XGBoost, ElmanNN, and ELM have been widely applied to estimate the air pollutant concentrations (Chen et al., 2018c; Zang et al., 2018; Zhu et al., 2019), these methods suffered from some weaknesses in predicting the pollutant level. For instance, both the BPNN and ElmanNN models could capture the locally optimal solution when the training subsets were integrated into the final model, which decreased the predictive performance of the model (Wang et al., 2015). Moreover, BPNN generally showed slow training speed, especially with the huge training subsets (Li and Park, 2009; Wang et al., 2015). ELM often consumed more computing resources and experienced the overfitting issue due to the increase in sampling size (Huang et al., 2015; Shao et al., 2015). The GRNN method advanced the training speed compared with the BPNN model and avoided the locally optimal solution during the modelling process (Zang et al., 2019), whereas the predictive performance is still worse than that of the RF model. XGBoost was often considered to be robust in predicting the air pollutant level (Li et al., 2020), while the model did not display excellent performance in the present study. This might be attributable to the sampling size in the present study not being big enough because the model generally showed better performance with big samples. Moreover, we found that the two-stage model was superior to the one-way model in the predictive performance. This encouraging result suggested that the relationship between the predictors and the 8 h O₃ concentration varied with space and time. The two-stage model used the GAM method to further adjust the prediction error of the RF model and considered the spatiotemporal correlation of the predictor error on the Tibetan Plateau. Although the STK model incorporated space and time into the model simultaneously, the RF–GAM model outperformed the RF–STK model. It was assumed that the STK model showed higher uncertainty in predicting the O₃ concentration in regions with few sampling sites (Gao et al., 2016; Li et al., 2017a). Overall, the ensemble RF–GAM model showed significant improvement in predictive performance.

Table 2 The R² values, RMSE, MPE, and RPE of RF–GAM in four seasons over the Tibetan Plateau.

Table 3 The predictive performances of RF–GAM in different provinces over the Tibetan Plateau.

Figure 6 The mean value of estimated 8 h O₃ concentration during 2005–2018 over the Tibetan Plateau.

The performances of RF–GAM displayed slight differences for each year during 2014–2018. As shown in Table 1, the R² value showed the highest value (0.76) in 2016, followed by that in 2018 (0.75), 2017 (0.73), and 2015 (0.72), and showed the lowest one in 2014 (0.69). Both the RMSE and MPE exhibited the lowest values in 2014, while these parameters did not show significant variation during 2015–2018. The lowest R² value and the highest RPE were found in 2014 due to having the lowest sample size, while the highest R² value and lowest RPE in 2016 were due to the maximum sample size. Geng et al. (2018) found that the predictive performance of the machine-learning model was strongly dependent on the number of training samples and sampling frequency. The lower RMSE and MPE in 2014 might be attributable to the lack of measured O₃ data in spring, which decreased the higher value of O₃ concentration. The performances of the RF–GAM model in four seasons were also assessed by 10-fold cross validation (Table 2). The predictive performance of the RF–GAM model showed significant seasonal differences, with the highest R² value observed in summer (0.74), followed by winter (0.69) and autumn (0.67), and the lowest one in spring (0.64). However, both the RMSE and MPE displayed different seasonal characteristics from the R² value. Both the RMSE and MPE for RF–GAM followed the order of spring (15.32 and 11.94 µg m⁻³) > summer (15.13 and 11.75 µg m⁻³) > winter (14.58 and 11.44 µg m⁻³) > autumn (13.23 and 10.52 µg m⁻³). The lowest R² value in spring might be caused by multiple O₃ sources and complicate O₃ formation mechanisms. On the one hand, the O₃ in spring might be generated from the local anthropogenic emission or long-range transport (Li et al., 2017b, 2019b). On the other hand, a strong stratosphere–troposphere exchange process due to the lower height of the troposphere on the Tibetan Plateau might lead to the higher O₃ concentration in spring (Skerlak et al., 2014). Unfortunately, both the long-range transport and stratosphere–troposphere exchange process were missing in the RF–GAM model, which restricted the accuracy of O₃ estimation in spring. The large estimation errors (e.g. RMSE, MPE, and RPE) in spring and summer were attributable to the high 8 h O₃ concentration in these seasons, while the low prediction error observed in autumn was due to the low O₃ level.

Figure 7 The inter-annual variation in predicted 8 h O₃ level (µg m⁻³) from 2005 to 2018 across the Tibetan Plateau.

Apart from the seasonal variation, we also investigated the spatial variabilities in the predictive accuracy for the RF–GAM model. The Tibetan Plateau was classified into five provinces, and then the predictive performance of RF–GAM model in each province was calculated. Among the five provinces, Gansu displayed the highest R² value (0.74), followed by Sichuan Province (0.71), Qinghai Province (0.70), the autonomous region of Tibet (0.69), and Yunnan Province (0.54) (Table 3). The results shown herein were not in agreement with the previous studies by Geng et al. (2018), who confirmed that the predictive performance of the machine-learning model was positively associated with the sampling size. It was assumed that the spatial distribution of the sampling sites in Tibet was uneven and the sampling density was low, though Tibet possessed the highest number of monitoring sites of the provinces. The prediction errors (RMSE and MPE) did not exhibit the same characteristics as the R² value. The higher RMSE and MPE were found in the autonomous region of Tibet (14.81 and 11.24 µg m⁻³) and Qinghai Province (14.83 and 11.33 µg m⁻³) due to the higher values of blh and sund. The lowest values of the RMSE and MPE could be observed in Yunnan Province, which was due to the higher rainfall amount. The highest RPE was found in Yunnan Province (25.85 %), followed by Tibet (22.90 %), Qinghai (22.65 %), and Sichuan (22.62 %), and the lowest one was found in Gansu Province (22.51 %), which might be linked with the sample size.

Figure 8 The trend analysis of predicted 8 h O₃ concentration. Panels (a) and (b) represent the result of Mann–Kendall method and discrepancy of estimated O₃ level during 2005–2018 across the Tibetan Plateau.

Although 10-fold cross-validation verified that the RF–GAM model showed better predictive performance in estimating the surface 8 h O₃ concentration, the test method cannot validate the transferring ability of the final model. The monitoring site on the Tibetan Plateau before May 2014 is very limited, and only the daily 8 h O₃ data in Lhasa from the open-access website (https://www.aqistudy.cn/historydata/daydata.php?city=%E6%8B%89%E8%90%A8&month=2013-12, last access: 19 May 2020) were available to compare with the simulated data. As depicted in Fig. 4, the R² value of the unlearning 8 h O₃ level against the predicted 8 h O₃ concentration reached 0.67, which was slightly lower than that of the 10-fold cross-validation R² value. Overall, the extrapolation ability of the RF–GAM model is satisfactory, and thus it was assumed that the model could be applied to estimate the O₃ concentration in other years. Both the RMSE and MPE for the unlearning 8 h O₃ level against the predicted 8 h O₃ concentration were significantly higher than those of the 10-fold cross validation. It was assumed that Lhasa showed a higher surface 8 h O₃ concentration over the Tibetan Plateau.

To date, some previous studies also simulated the surface O₃ concentration on the Tibetan Plateau using statistical models (Zhan et al., 2018). For instance, Zhan et al. (2018) employed the RF–STK model to estimate the surface O₃ concentration over China and explained the 66 % spatial variability in the O₃ level on the Tibetan Plateau. Apart from these statistical models, some classical CTMs were also applied to estimate the O₃ concentration in remote areas. Both Liu et al. (2018) and Lin et al. (2018) used CMAQ to estimate the O₃ level across China, while the R² values in most of cities were lower than 0.50. In terms of the predictive performance, the RF–GAM model in our study showed significant advantages compared with previous studies. It should be noted that our RF–GAM model could outperform most of current models, chiefly because of (1) accounting for the temporal autocorrelation of the surface O₃ concentration and (2) the use of high-quality satellite data.

Figure 9 The seasonal variability in estimated 8 h O₃ level across the Tibetan Plateau. Panels (a–d) represent the predicted 8 h O₃ concentrations in spring, summer, autumn, and winter, respectively.

Figure 10 The spatial distributions of nonattainment days on the Tibetan Plateau during 2005–2018.

3.2 Variable importance

The results of variable importance for key variables are depicted in Fig. 5. In the final RF–GAM model, it was found that time was the dominant factor for the 8 h O₃ concentration on the Tibetan Plateau, indicating that the ambient O₃ concentration displayed significant temporal correlation. Following the time, meteorological factors served as the main factors for the O₃ pollution in remote regions. The sum of sund, sp, d2m, t2m, and tp constituted 34.43 % of the overall variable importance. Among others, sund was considered to be the most important meteorological factor for the O₃ pollution. It was assumed that strong solar radiation and long duration of sunshine favoured the photochemical generation of ambient O₃ (Malik and Tauler, 2015; Stähle et al., 2018). Tan et al. (2018) demonstrated that the chemical reaction between NO_x and VOCs was strongly dependent on the sunlight. Besides this, the atmospheric pressure (sp) was also treated as a major driver for the O₃ pollution over the Tibetan Plateau. Santurtún et al. (2015) have demonstrated that sp was closely linked to the atmospheric circulation and synoptic-scale meteorological pattern, which could influence the long-range transport of ambient O₃. Apart from sund and sp, d2m and t2m played significant role in the O₃ pollution, which was consistent with many previous studies (Zhan et al., 2018). Zhan et al. (2018) observed that cold temperatures were not favourable to the O₃ formation. d2m can affect the surface O₃ pollution through two aspects. On the one hand, RH affected heterogeneous reactions of O₃ and particles (e.g. soot, mineral) (He et al., 2017; He and Zhang, 2019; Yu, 2019). On the other hand, high RH could increase the soil moisture and evaporation, and thus the water-stressed plants tended to emit more biogenic isoprene, thereby promoting the elevation of O₃ concentration (Zhang and Wang, 2016). It should be noted that the effect of precipitation on O₃ pollution was weaker than those of other meteorological factors. Zhan et al. (2018) also found a similar result and believed that rain scavenging served as the key pathway for the O₃ removal only when O₃ pollution was very serious. The effect of the O₃ column amount on surface O₃ concentration seemed to be lower than those of most meteorological factors, suggesting that vertical transport of ambient O₃ was complex. Although socioeconomic factors and land use types were not dominant factors for the O₃ pollution on the Tibetan Plateau, they still cannot be ignored in the present study because the predictive performance would worsen if these variables were excluded from the model. It was widely acknowledged that the anthropogenic emissions focused on developed urban areas with high population density, especially on the remote plateau (Zhang et al., 2007; Zheng et al., 2017). Compared with the urban land, the grassland played a more important role in the O₃ pollution on the Tibetan Plateau. It was thus assumed that the grassland was widely distributed on the Tibetan Plateau, which could release a large amount of biogenic volatile organic compounds (BVOCs) (Fang et al., 2015). It was well known that photochemical reactions of BVOCs and NO_x in the presence of sunlight caused the O₃ formation (Calfapietra et al., 2013; Yu et al., 2006). Furthermore, Fang et al. (2015) confirmed that BVOC emission on the Tibetan Plateau displayed a remarkable increase in the wet seasons.

3.3 The spatial distribution of estimated 8 h O₃ concentration over the Tibetan Plateau

Figure 6 depicts the spatial distribution of the 8 h O₃ level estimated by the novel RF–GAM model. The spatial distribution pattern modelled by the RF–GAM model showed a similar characteristic to the results found by previous studies except on the northern Tibetan Plateau (Liu et al., 2018). The estimated 8 h O₃ concentration displayed the highest value in some cities of the northern Tibetan Plateau, such as Huangnan (73.48±4.53 µg m⁻³) and Hainan (72.24±5.34 µg m⁻³), followed by the cities in the central region, including Lhasa (65.99±7.24 µg m⁻³) and Shigatse (65.15±6.14 µg m⁻³), and the lowest one was found in a city of the southeastern Tibetan Plateau (Aba) (55.17±12.77 µg m⁻³). The spatial pattern of the 8 h O₃ concentration is highly consistent with the result predicted by Liu et al. (2018) using the CMAQ model, while it is not in agreement with the result estimated by Zhan et al. (2018) using the RF–STK model. The difference between the present study and Zhan et al. (2018) is seen on the northern Tibetan Plateau, which lacks a monitoring site and still has the higher uncertainty. Firstly, this might be due to the weakness of RF–STK mentioned above. Moreover, Zhan et al. (2018) only used the ground-level-measured data in 2015 to establish the model, and the data from new sites since 2015 were not incorporated into the model, which could increase the model uncertainty (Zhan et al., 2018). As shown in Fig. 6, most of the cities in Qinghai Province (e.g. Huangnan, Hainan, and Guoluo) generally showed a higher 8 h O₃ concentration over the Tibetan Plateau, which was in a good agreement with the spatial distribution of the O₃ column amount (Fig. S3). Besides this, some cities in Tibet, such as Shigatse and Lhasa, also showed higher 8 h O₃ levels. It was assumed that the precursor emissions in these regions were significantly higher than those in other cities of the Tibetan Plateau (Fig. S4). Zhang et al. (2007) used the satellite data to observe that the higher VOCs and NO_x emission was concentrated in the residential areas with high population density on the remote Tibetan Plateau. Apart from the effect of anthropogenic emission, the meteorological conditions could also be important factors for the 8 h O₃ concentration. As shown in Figs. S5–S10, a higher blh and sp on the northeastern Tibetan Plateau might promote the O₃ formation through the reaction of the VOC and OH radical, leading to a higher 8 h O₃ concentration in these cities (Ou et al., 2015). In addition, a lower tp occurred on the northern Tibetan Plateau and the northeastern Tibetan Plateau, both of which were unfavourable to the ambient O₃ removal (Yoo et al., 2014). In contrast, the higher tp observed on the southeastern Tibetan Plateau resulted in slight O₃ pollution.

Table 4 The estimated 8 h O₃ concentration in 19 prefecture-level cities over the Tibetan Plateau during four seasons, including spring, summer, autumn, and winter.

3.4 The temporal variation in the simulated 8 h O₃ concentration over the Tibetan Plateau

The annual mean estimated 8 h O₃ concentration on the Tibetan Plateau displayed a slow increase, from 64.74±8.30 to 66.45±8.67 µg m⁻³ in 2005–2015 (Table S1), whereas it decreased from the peak to 65.87±8.52 µg m⁻³ during 2015–2018 (Fig. 7). Based on the Mann–Kendall method (Fig. 8a), it was found that the surface O₃ concentration exhibited a slight increase on the whole, while the degree of increase was not significant (p>0.05). Besides this, it should be noted that the O₃ concentrations in various regions showed different rates of increase. As depicted in Fig. 8b, we found that the 8 h O₃ concentrations on the northern, western, and eastern Tibetan Plateau displayed significant an increasing trend at the rate of 1–3 µg m⁻³ during 2005–2018. The middle region of the Tibetan Plateau showed a moderate increase trend at the rate of 0–1 µg m⁻³. However, the 8 h O₃ concentration in Shigatse and Shannan even displayed a decreasing trend in 2005–2018.

Besides this, the 8 h O₃ concentrations on the Tibetan Plateau displayed significantly seasonal discrepancy. The estimated 8 h O₃ level on the Tibetan Plateau followed the order of spring (75.00±8.56 µg m⁻³) > summer (71.05±11.13 µg m⁻³) > winter (56.39±7.42 µg m⁻³) > autumn (56.13±8.27 µg m⁻³) (Fig. 9 and Table 4). The 8 h O₃ concentrations in most of prefecture-level cities showed similar seasonal characteristics, with overall seasonal variation on the Tibetan Plateau. Based on the result summarised in Table S2, it was found that the key precursors of ambient O₃ generally displayed higher emissions in winter compared with other seasons. However, the seasonal distribution of ambient O₃ concentration was not in accordance with the precursor emissions, suggesting that the meteorological factors might play more important roles in ambient O₃ concentration. It was well known that the higher air temperatures in spring and summer were closely related to the low sp and high sund, both of which promoted O₃ formation (Sitnov et al., 2017). Although summer showed the highest air temperature and the longest sunshine duration, the higher rainfall amount in summer decreased the ambient O₃ concentration via wet deposition (Li et al., 2017a, 2019b). Moreover, the highest blh occurred in spring, which was favourable to the strong stratosphere–troposphere exchange process on the Tibetan Plateau (Skerlak et al., 2014). Therefore, the 8 h O₃ concentrations in summer and winter were lower than that in spring. Nonetheless, the 8 h O₃ levels in Diqing, Shannan, and Nyingchi displayed the highest values in spring (56.38±7.87, 73.90±5.97, and 73.22±2.77 µg m⁻³), followed by winter (45.88±7.05, 61.71±4.32, and 62.24±3.63 µg m⁻³) and summer (44.35±5.90, 61.00±5.86, and 59.60±2.33 µg m⁻³), and the lowest ones in autumn (37.45±5.76, 54.70±3.13, and 53.84±2.06 µg m⁻³). The lower O₃ level in summer than winter was mainly attributable to the higher precipitation observed in the summer of these cities (Fig. S11). In addition, it should be noted that the NO_x and VOC emissions of the southern Tibetan Plateau (e.g. Shannan) exhibited higher values in winter compared with other seasons.

3.5 The nonattainment days over the Tibetan Plateau during 2005–2018

The annual mean nonattainment days in the 19 prefecture-level cities over the Tibetan Plateau are summarised in Table 2. The value of 100 µg m⁻³ was regarded as the critical value for the 8 h O₃ level by the World Health Organization (WHO). Nonattainment days refer to total days with the 8 h O₃ concentration higher than 100 µg m⁻³. Although the annual mean 8 h O₃ concentrations in all of the cities over the Tibetan Plateau did not exceed the critical value, not all of the regions experienced excellent air quality in the long term (2005–2018). Some cities of Qinghai Province, including Huangnan, Haidong, and Guoluo, suffered from 45, 40, and 40 nonattainment days each year (Fig. 10 and Table 5). Besides this, some cities on the southern Tibetan Plateau, such as Shigatse and Shannan, also experienced more than 40 nonattainment days each year, suggesting that the Tibetan Plateau still faced the risk of O₃ pollution. Fortunately, some remote cities, such as Ali, Ngari, and Qamdo, did not experience excessive O₃ pollution all the time, which was ascribed to low precursor emissions and appropriate meteorological conditions. It should be noted that the nonattainment days in regions with high O₃ concentration showed significant seasonal difference, whereas the seasonal difference was not remarkable in cities with low O₃ pollution. As shown in Table 2, it should be noted that nearly all of the nonattainment days could be detected in spring and summer, which was in good agreement with the O₃ levels in different seasons, indicating that the O₃ pollution issue should be given more attention in spring and summer.

Table 5 The mean nonattainment days (8 h O₃ level >100 µg m⁻³) in 19 prefecture-level cities over the Tibetan Plateau each year.

The determination of nonattainment days showed some uncertainties, owing to the predictive error of modelled O₃ concentration. First of all, meteorological data used in RF–GAM model were collected from reanalysis data, and these gridded data often showed some uncertainties, which could increase the uncertainty of O₃ estimation. Second, the O₃ column amount used in the present study reflected the vertical O₃ column amount rather than the surface O₃ concentration. Thus, it could decrease the predictive performance of the surface O₃ level.

4 Summary and implications

In the present study, we developed a novel hybrid model (RF–GAM) based on multiple explanatory variables to estimate the surface 8 h O₃ concentration across the remote Tibetan Plateau. The 10-fold cross-validation method demonstrated that RF–GAM achieved excellent performance, with the highest R² value (0.76) and lowest root-mean-square error (RMSE) (14.41 µg m⁻³), compared with other models, including the RF–STK, RF, BPNN, XGBoost, GRNN, ElmanNN, and ELM models. Moreover, the unlearning ground-level-measured O₃ data validated the fact that the RF–GAM model showed better extrapolation performance (R²=0.67, RMSE =25.68 µg m⁻³). The result of variable importance suggested that time, sund, and sp were key factors for the surface 8 h O₃ concentration over the Tibetan Plateau. Based on the RF–GAM model, we found that the estimated 8 h O₃ concentration exhibited notable spatial variation, with higher values in some cities of the northern Tibetan Plateau, such as Huangnan (73.48±4.53 µg m⁻³) and Hainan (72.24±5.34 µg m⁻³), and lower values in some cities of the southeastern Tibetan Plateau, such as Aba (55.17±12.77 µg m⁻³). Besides this, we also found that the O₃ level displayed a slow increase, from 64.74±8.30 to 66.45±8.67 µg m⁻³ from 2005 to 2015, while the O₃ concentration decreased to 65.87±8.52 µg m⁻³ in 2018. The estimated 8 h O₃ level on the Tibetan Plateau showed significant seasonal discrepancy in the order of spring (75.00±8.56 µg m⁻³) > summer (71.05±11.13 µg m⁻³) > winter (56.39±7.42 µg m⁻³) > autumn (56.13±8.27 µg m⁻³). Based on the critical value set by the WHO, most of the cities on the Tibetan Plateau had excellent air quality, while several cities (e.g. Huangnan, Haidong, and Guoluo) still suffered from more than 40 nonattainment days each year.

The RF–GAM model for O₃ estimation has several limitations. First of all, the O₃ estimation of the northern Tibetan Plateau might show some uncertainties because there are few ground-level monitoring sites, and thus we cannot validate the reliability of predicted values in regions without a monitoring site. Secondly, our approach did not include data on the emission inventory or traffic count because the continuous emissions of NO_x and VOCs were not open access. Lastly, we only focused on the temporal variation in the surface O₃ concentration in the past 10 years, and the short-term O₃ data cannot reflect the response of O₃ pollution to climate change. In the future work, we should combine more explanatory variables such as long-term NO_x and VOC emissions to retrieve the surface O₃ level over the Tibetan Plateau in recent decades.