A grid model for vertical correction of precipitable water vapor over the Chinese mainland and surrounding areas using random forest

Various ground-based observing techniques provide precipitable water vapor (PWV) products with different spatial resolutions. To effectively integrate these products, especially in terms of vertical orientation, spatial interpolation is essential. In this context, we have developed a model to characterize PWV variation with altitude over our study area. Our model, known as RF-PWV (a PWV vertical correction grid model with a 1° × 1° resolution), is constructed using random forest based on the relationship between the differences in different pressure level PWV data from the fifth-generation European Centre for Medium-Range Weather Forecasts reanalysis (ERA5) monthly average hourly data and corresponding differences in their height differences over time. When validated against 1 h ERA5 PWV profiles, RF-PWV exhibits a 99.84 % reduction in bias and a 63.41 % decrease in the RMSE compared with the most recent model, C-PWVC1. Furthermore, when validated against radiosonde data, RF-PWV shows a 96.36 % reduction in bias and a 5 % decrease in the RMSE compared with C-PWVC1. Additionally, RF-PWV outperforms C-PWVC1 in terms of resistance to seasonal and height difference interference. The model eliminates the need for meteorological parameters, allowing for high-precision PWV vertical correction by inputting only time and height differences. Consequently, RF-PWV can significantly reduce errors in vertical correction, enhance PWV fusion product accuracy, and provide insights into PWV vertical distribution, thereby contributing to climate research.

1 Introduction

Precipitable water vapor (PWV), the most abundant greenhouse gas, primarily resides in the troposphere and plays a pivotal role in the global energy budget, hydrological cycle, and climate change (Zhang et al., 2018; Li et al., 2022b; Dessler and Sherwood, 2009; Raval and Ramanathan, 1989; Rocken et al., 1997). Various observation platforms, including radiosondes (RSs), microwave water vapor radiometers (WVRs), satellite remote sensing, ground-based global navigation satellite system (GNSS) data, and reanalysis data, have amassed extensive PWV data through long-term data accumulation (Huang et al., 2022). Combining multisource data enables more accurate and comprehensive water vapor monitoring and meteorological research (B. Zhang et al., 2019; Li et al., 2022a; Alshawaf et al., 2015). However, inconsistent pressure levels (heights) for storing PWV data from different sources hinder the fusion and reliability analysis of PWV multisource data (Chen et al., 2023b; Yang et al., 2023). Therefore, precise PWV vertical corrections are indispensable with respect to the utilization of PWV fusion products. The vertical distribution of PWV is closely related to the formation and distribution of rainfall and clouds, which is the pivot of weather forecasting and is also one of the factors affecting convection and monsoon climates (Bevis et al., 1992; Keil et al., 2008; Rose and Rencurrel, 2016). The vertical distribution and temporal variability of PWV are essential for understanding regional weather and global climate, improving the climate models, and predicting future climate change (Jacob, 2001; Renju et al., 2015). Hence, proposing a more accurate and applicable PWV vertical correction model is of paramount importance.

Common methods for PWV vertical correction involve establishing empirical vertical correction models to enhance the applicability of PWV vertical correction (Emardson and Johansson, 1998; Dousa and Elias, 2014; Huang et al., 2023). Reitan (1963) introduced an empirical formula describing water vapor density's exponential decrease in the vertical direction, based on the relationship between PWV near the surface and at high altitudes. The PWV lapse rate (−0.5 mm km⁻¹), estimated by Kouba (2008) using the International GNSS Service (IGS) and the Vienna Mapping Function 1 (VMF1), has been widely adopted. However, considering the seasonal variations in the PWV lapse rate as constant introduces significant errors in PWV vertical correction (Tomasi, 1977; Leckner, 1978; H. X. Zhang et al., 2019; Zhang et al., 2022). Huang et al. (2021) developed a PWV vertical correction model that accounts for seasonal variations in the PWV lapse rate, offering greater accuracy and stability than the classic PWV vertical correction model (PWV lapse rate = −0.5 mm km⁻¹) in China. Wang et al. (2022) incorporated spherical harmonic functions to develop a PWV vertical correction model, achieving high accuracy in the Tibetan Plateau. Nevertheless, many existing models assume PWV's exponential decrease and represent PWV lapse rate variations using periodic functions, failing to address complex nonlinear variations beyond daily/sub-daily and seasonal variations in the PWV lapse rate.

Neural network techniques are well suited for handling nonlinear problems and have found applications in various industries (Zheng et al., 2022). Machine learning has demonstrated promising potential in modeling tropospheric parameters (Ravuri et al., 2021; Lam et al., 2022). Senkal (2015) developed a model for predicting PWV in Turkey (Türkiye) using a resilient propagation (RP) neural network, which provides PWV estimates for a given location. Validation with RS PWV data in the study area revealed good agreement between the new model and RS PWV data. Zhu et al. (2022) created a weighted-mean temperature (Tm) vertical correction grid model (CTm-FNN) employing a feedforward neural network in China. This model outperformed the Chinese Tropospheric Model (CTrop) and Global Pressure and Temperature 3 (GPT3), reducing the root-mean-square error (RMSE) by 86 % and 83 %, respectively.

Therefore, this paper presents a random-forest-based precipitable water vapor vertical correction grid model, termed RF-PWV, for China and surrounding areas, harnessing random forest's powerful nonlinear fitting capability and the high temporal resolution of monthly average hourly PWV data. With RF-PWV, PWV differences can be obtained by simply inputting time and height differences, allowing for high-precision PWV vertical correction. The model offers PWV vertical correction techniques for multisource PWV fusion, weather forecasting, and climate studies.

We begin by providing an overview of the study area and the experimental dataset. Subsequently, we describe the data processing strategy and modeling methodology. Next, we evaluate the performance of the RF-PWV model. Finally, we conclude our study and outline future directions.

2 Data and methods

2.1 Study area

The study area includes the region between 15 and 55° N and between 70 and 135° E, covering mainland China and its surrounding areas, characterized by extensive land and ocean. China's topography exhibits significant variation, with higher elevations in the west gradually sloping to lower elevations in the east. Influenced by the monsoon climate, the summer monsoon brings substantial moisture from the ocean into the region, while winter introduces cold, dry air inland (Sun et al., 2019; W. X. Zhang et al., 2019). These geographical and climatic factors contribute to a complex spatiotemporal variation in PWV. As a result, the vertical distribution of PWV in this area presents a challenging problem to characterize, making it a suitable choice for our experimental area.

2.2 Datasets

2.2.1 ERA5 PWV

ERA5, the fifth-generation atmospheric reanalysis product developed by the European Centre for Medium-Range Weather Forecasts (ECMWF), offers access to 1 h meteorological data across 37 pressure levels, with a horizontal resolution as fine as 0.25° × 0.25°. This dataset can be downloaded from https://cds.climate.copernicus.eu/ (last access: 29 March 2024​​​​​​​) (Albergel et al., 2018). ERA5 is renowned for its superior accuracy compared with its predecessor, ERA-Interim, and has gained widespread usage in meteorological research (Hersbach et al., 2020; Lu et al., 2023; Chen et al., 2023a). Moreover, the monthly averaged dataset, in terms of accuracy, rivals the daily dataset while demonstrating greater stability (Dogan and Erdogan, 2022). Additionally, the monthly average hourly dataset offers the advantage of capturing both seasonal variations in meteorological data and finer-grained sub-daily variations. In this study, we utilize the monthly average hourly dataset, which provides 1 h data at 37 pressure levels with a spatial resolution of 1° × 1°. The PWV for each pressure level is determined through integration, as described by Y. L. Zhang et al. (2019) and X. M. Wang et al. (2016):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{PWV}}_i={\sum }_i^{n-\mathrm{1}}{\displaystyle \frac{(q_i+q_{i+\mathrm{1}})\cdot (p_{i+\mathrm{1}}-p_i)}{\mathrm{2}\cdot {\mathit{\rho }}_{\mathrm{w}}\cdot }},\text{(2)}& {\displaystyle }g=\mathrm{9.780325}\cdot {\left[{\displaystyle \frac{\mathrm{1}+\mathrm{0.00193185}\cdot \sin (\mathit{\phi }{)}^{\mathrm{2}}}{\mathrm{1}-\mathrm{0.00669435}\cdot \sin (\mathit{\phi }{)}^{\mathrm{2}}}}\right]}^{\mathrm{0.5}}.\end{array}$$

Here, n represents the total number of layers; PWV_i, q_i, and p_i represent the PWV (mm), specific humidity (kg kg⁻¹), and pressure (Pa) at the i layer, respectively; ρ_w is the density of liquid water, which is standardized to 1000 kg m⁻³; g is the gravitational acceleration (m s⁻²); and φ denotes the latitude (rad).

It is crucial to emphasize that the upper boundary of the troposphere lies at approximately 10 km altitude (Ding, 2020). Consequently, PWV effectively approaches 0 mm when situated at elevations exceeding 12 km vertically. As a result, we restrict our PWV calculations to cover pressure levels within the range of 0 to 12 km above the grid point for all subsequent analyses and investigations.

2.2.2 RS PWV

RSs are widely recognized for their high-precision PWV measurements and are commonly considered a reference standard for evaluating other measurement techniques (Adeyemi and Joerg, 2012; Wang et al., 2021; Zhao et al., 2022). We obtained RS PWV data from the Integrated Global Radiosonde Archive (IGRA), accessible at https://www1.ncdc.noaa.gov/pub/data/IGRA (last access: 29 March 2024​​​​​​​), with a temporal resolution of 12 h. We made use of meteorological data from 148 stations, focusing on pressure levels within the 0–12 km altitude range (as illustrated in Fig. 1). The specific humidity at each pressure level was determined by employing Eqs. (3) and (4), which are as follows (Zhai and Eskridge, 1997; Ross and Elliott, 1996):

$$\begin{array}{}\text{(3)}& {\displaystyle }e={\displaystyle \frac{\mathrm{RH}\cdot e_{\mathrm{s}}}{\mathrm{100}}},\text{(4)}& {\displaystyle }q={\displaystyle \frac{\mathrm{0.622}\cdot e}{p-\mathrm{0.378}\cdot e}}.\end{array}$$

Here, RH represents relative humidity (%), e_s signifies saturated vapor pressure (Pa), e denotes water vapor pressure (Pa), and q represents specific humidity (kg kg⁻¹). Subsequently, RS PWV values for various pressure levels were calculated using Eq. (1).

Figure 1 Distribution of the selected radiosonde sites.

2.3 Establishment of the RF-PWV model

Random forest (RF), an ensemble learning method that combines multiple weak learners to form a single strong learner, typically improves generalization performance and model robustness (Breiman, 2001; Sagi and Rokach, 2018). Compared with the backpropagation neural network (BPNN), random forests are less prone to overfitting, especially with noisier datasets like PWV. Random forests handle noisy data and outliers more efficiently, making new models more robust and often easier to tune (L. A. Wang et al., 2016; Tyralis et al., 2019). In addition, our previous study has shown that RF outperforms BPNN with respect to modeling spatiotemporal variability in tropospheric parameters (Li et al., 2023a). Thus, RF is employed to model the height dependency of PWV. The equation governing random forest's prediction is expressed as follows:

$$\begin{array}{}\text{(5)}& Y\left(X\right)={\displaystyle \frac{\mathrm{1}}{B}}{\sum }_{b=\mathrm{1}}^B{T}_b\left(X\right).\end{array}$$

Here, Y(X) is the final prediction result, T_b(X) represents the predicted value of each decision tree, and B denotes the number of decision trees. The selection of an appropriate number of decision trees is pivotal in modeling; too few trees may lead to overfitting, whereas too many trees can result in excessively long modeling times (Sun et al., 2021; Probst and Boulesteix, 2017).

2.3.1 Defining the primary parameter

To assess the performance of machine learning models, 10-fold cross-validation is a commonly employed technique (Rodriguez et al., 2010; Zhang and Yao, 2021). In this context, 10-fold cross-validation was employed to ascertain the optimal number of decision trees based on the RMSE. The fundamental principle of 10-fold cross-validation entails dividing the input data into 10 groups. Subsequently, nine randomly selected groups are utilized as the training set, and the remaining group serves as the validation set. This process is reiterated 10 times to ensure that all data are included in both training and testing. This approach provides results that closely approximate the accuracy of the final model while guarding against overfitting (Santos et al., 2018). Based on our experience, we experimented with decision tree numbers ranging from 5 to 95, with a step size of 10, to train the model and evaluate its performance under varying decision tree quantities (Li et al., 2023a). The results, depicted in Fig. 2, exhibit a significant decline in the RMSE as the number of decision trees increases from 5 to 75, reaching a minimum at 75. However, increasing the number of trees to 75 does not significantly enhance accuracy, and it incurs longer training times. In consideration of the need for modeling at multiple grid points and balancing fitting quality with training time, a final decision was made to employ 55 trees to build the model.

Figure 2 Cross-validation RMSE for different numbers of decision trees.

2.3.2 Training the model

During the model training phase, we performed individual modeling at each grid point (1° × 1°) using ERA5 monthly average hourly PWV data at pressure levels ranging from 1000 to 225 hPa, within the 0 to 12 km altitude range, spanning the years 2008 to 2017. The ith PWV differences ($\mathrm{\Delta }{\mathrm{PWV}}_i={\mathrm{PWV}}_i-{\mathrm{PWV}}_{\mathrm{0}}$, where 0 stands for the bottom level) between the i level and the bottom level and the responding height differences ($\mathrm{\Delta }{\mathrm{GPH}}_i={\mathrm{GPH}}_i-{\mathrm{GPH}}_{\mathrm{0}}$) were all computed and utilized as the training dataset. In essence, each grid point contained 63 360 samples (22 levels × 24 h × 12 months × 10 years), and the region consisted of 2706 grid points (66 longitudes × 41 latitudes) at a 1° × 1° resolution. The model, denoted as the RF-PWV model, characterizes the relationship between ΔPWV and ΔGPH, as illustrated in Fig. 3. The input data included year, day of the year (doy is the first day of the corresponding month), hour of the day (hod), and ΔGPH; the output data were ΔPWV. The reason why “year” was included as an input variable to RF-PWV is that the PWV lapse rate has a significant periodic function with a year (Du et al., 2023; Huang et al., 2023). When users employ the model, they are only required to provide the geopotential height of the target point, the datum PWV, the time (year, doy, hod), and the height difference of the target point concerning the datum point (ΔGPH). They can obtain the corresponding ΔPWV and can then get the PWV of the target height by adding the datum PWV to the ΔPWV.

Figure 3 Network structure of the RF-PWV model based on the random forest algorithm.

In the application of the RF-PWV model, the four grid points surrounding the target point are determined based on the target point's geographical coordinates (latitude and longitude). Then, the ΔPWV at the corresponding height of the four selected points is calculated using the RF-PWV model. Finally, the ΔPWV at the target point's location is determined through bilinear interpolation. This process involves calculating the difference between the target point's GPH and the reference station's GPH₀ to get the ΔGPH. Next, the time information is input into the models for the four nearest grid points to the target point, yielding the ΔPWV at the corresponding height of these grid points. Finally, bilinear interpolation is employed to calculate the ΔPWV at the target point's location. This method offers the advantage of not requiring an exceptionally strong spatial generalization ability for a single model. It comprehensively considers the relationship between the target point and the four nearest grid points within the limited spatial context, resulting in enhanced consistency and higher accuracy at each grid point, ensuring the overall model's robustness.

3 Accuracy validation and analysis

To validate the RF-PWV model, we employed hourly ERA5 and RS pressure level data from the study area in 2018 as the test set, while also selecting a newly developed PWV vertical correction model (C-PWVC1) for comparison. Note that the authors of the C-PWVC model suggest using C-PWVC1 directly for PWV vertical correction in the study area, so C-PWVC2 is ignored. C-PWVC1 has been proven to be more accurate than the classical PWV vertical correction model (PWV lapse rate = −0.5 mm km⁻¹) in the study area (Huang et al., 2021). C-PWVC1 is a model using the exponential function to account for the height dependency of PWV. C-PWVC1 can be expressed as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathrm{PWV}}_{h_{\mathrm{1}}}={\mathrm{PWV}}_{h_{\mathrm{2}}}\cdot \exp \left(\mathit{\beta }\left(h_{\mathrm{1}}-h_{\mathrm{2}}\right)\right),\text{(7)}& {\displaystyle }\begin{array}{rl}\mathit{\beta }\left(\mathrm{doy}\right)& =-\mathrm{0.35}-\mathrm{0.026}\cos \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{2}\mathit{\pi }\right)-\mathrm{0.015}\sin \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{2}\mathit{\pi }\right)\\ & +\mathrm{0.008}\cos \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{4}\mathit{\pi }\right)+\mathrm{0.026}\sin \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{4}\mathit{\pi }\right).\end{array}\end{array}$$

Here, ${\mathrm{PWV}}_{h_{\mathrm{1}}}$ and ${\mathrm{PWV}}_{h_{\mathrm{2}}}$ denote the PWV at h₁ and h₂, respectively; β is the PWV lapse rate; and doy is the day of the year. C-PWVC1 requires inputs of datum height, datum PWV, target height, and time to provide the PWV correction value at the target height, but the model is unable to capture nonlinear variations in the vertical direction. Then, the accuracy metrics employed for evaluation are bias and RMSE, as outlined below:

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left(X_i-X_i^{\prime }\right),\text{(9)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(X_i-X_i^{\prime }{)}^{\mathrm{2}}}.\end{array}$$

Here, X^′ is the reference values, X denotes model outputs, and n is the number of samples.

3.1 Validation of RF-PWV using ERA5 PWV

The RF-PWV model and C-PWVC1 were applied to vertically correct the hourly ERA5 bottom-level PWV data (1° × 1°) for the year 2018 to other pressure levels within the 0–12 km altitude range, excluding the bottom level. The results were then compared with ERA5 data, and the overall bias and RMSE are presented in Table 1. RF-PWV exhibited a bias close to 0 mm, indicating minimal systematic bias between the interpolated PWV and ERA5 PWV. Moreover, it reduced bias by 1.42 mm compared with C-PWVC1, corresponding to a remarkable optimization of 99.84 %. The bias values for RF-PWV were observed to fluctuate slightly within the range of −0.01 to 0.01 mm. Additionally, the RF-PWV RMSE showed a substantial reduction of 63.40 % compared with C-PWVC1. Furthermore, the RMSE values for RF-PWV demonstrated a more stable fluctuation pattern with a considerably narrower range. Overall, RF-PWV exhibited significantly higher accuracy than C-PWVC1, with corrected results showing better agreement with the reference values.

Table 1 Validation results of the RF-PWV and C-PWVC1 models tested using ERA5 data.

To provide a spatial illustration of the models' consistency with respect to accuracy, Fig. 4 displays the bias and RMSE values for each grid point for both RF-PWV and C-PWVC1. Notably, C-PWVC1 exhibited a significant north–south difference in bias, with larger values in the south and smaller values in the north. In contrast, RF-PWV demonstrated a substantial reduction in bias across almost all grid points in the study area, approaching 0 mm, effectively eliminating the north–south discrepancy. The most noteworthy improvement in accuracy was observed in the Tibetan Plateau and low-latitude regions. Despite the challenging climate conditions in the Tibetan Plateau and the strong land–sea interactions in the study area's low latitudes, which contribute to complex PWV variations, RF-PWV still achieved a bias close to 0 mm. These results highlight RF-PWV's adaptability to diverse weather conditions and its wide applicability. Furthermore, C-PWVC1 displayed a north–south difference in the RMSE. Higher RMSE values were concentrated in the southwestern and southeastern regions, reaching a maximum of 4.25 mm. This phenomenon is mainly attributable to the proximity of these regions to the ocean, frequent water vapor exchange between land and sea, and the complexity of PWV variations. However, RF-PWV's RMSE in these regions was significantly smaller than that of C-PWVC1, consistently measuring below 2 mm. Overall, RF-PWV's RMSE was lower than that of C-PWVC1 across the study area. Furthermore, RF-PWV exhibited excellent agreement, with values mostly hovering around 0.75 mm, nearly independent of spatial variations. These outcomes underscore the higher accuracy and improved spatial accuracy consistency of RF-PWV across the study area.

Figure 4 Distributions of bias and RMSE for RF-PWV and C-PWVC1 with respect to the ERA5 data.

To further evaluate the models' performance across different seasons, we calculated the bias and RMSE values for four representative grid points using data from 2018. These grid points were selected to represent various regions: 40.00° N, 80.00° E in the northwestern region; 15.00° N, 95.00° E in the southwestern region; 25.00° N, 110.00° E in the southeastern region; and 45.00° N, 125.00° E in the northeastern region. Figure 5a, b, g, and h illustrate that C-PWVC1 exhibited the highest bias and RMSE values during June–September, reaching 5.41 and 6.23 mm at 40.00° N, 80.00° E and 6.85 and 7.75 mm at 45.00° N, 125.00° E, respectively. Conversely, the lowest bias and RMSE values were recorded during January–February and November–December, hovering around 0 mm, with discernible seasonal fluctuations. This pattern is primarily attributed to significant PWV variations during the wet and rainy northern summers, compared with relatively mild PWV variations during the cold and dry winters. In contrast, Fig. 5c, d, e, and f show that the seasonal differences in the bias and RMSE for C-PWVC1 were less pronounced in the southern regions than in the northern regions. At 25.00° N, 110.00° E, which experiences abundant PWV changes and heavy rainfall throughout the year, the model's accuracy was relatively lower, with no noticeable seasonal variations. Similarly, near the Equator (15.00° N, 95.00° E), overall bias and RMSE values were more significant, with minimal seasonal differences. Notably, RF-PWV achieved substantially lower bias and RMSE values than C-PWVC1 during the summer months. Throughout the year, RF-PWV's bias and RMSE exhibited relatively stable patterns, with minimal fluctuations around 0 mm. Conversely, Fig. 5c, d, e, and f reveal that RF-PWV maintained bias and RMSE values around 0 mm, offering greater accuracy compared with C-PWVC1 in the southern grid points. In summary, RF-PWV exhibited enhanced resistance to seasonal variations, maintaining stable and accurate performance throughout the year across the study area.

Figure 5 Time series of RF-PWV and C-PWVC1 bias and RMSE for four selected grid points.

Given that more than three-quarters of the water vapor is concentrated in the lower atmosphere, in practice, most of the vertical correction of PWV occurs in the lower atmosphere (Yang et al., 2020). The bias and RMSE for C-PWVC1 and RF-PWV are statistically determined based on height differences, divided into 12 sections ranging from 0 to 6 km with intervals of 500 m. This division helps assess the applicability of the two models across different height segments, as shown in Fig. 6. Notably, C-PWVC1 exhibits a positive bias in every height difference segment, with the bias increasing as the height difference rises from 0 to 2.5 km, ultimately stabilizing at around 2.0 mm. RF-PWV bias tends to approach 0 mm in all height difference segments and shows negative bias after the height difference exceeds 2.5 km, with the absolute value increasing and reaching a maximum value of less than 0.2 mm. In each height difference segment between 0 and 6 km, RF-PWV bias is closer to 0 mm than C-PWVC1 bias, indicating that the corrected value of RF-PWV is more consistent with the reference value across different height difference segments. Additionally, the RF-PWV RMSE is significantly smaller than that of C-PWVC1 in all height difference segments. The RMSE for C-PWVC1 exhibits the same increasing trend as bias, stabilizing around 3 mm after the height difference exceeds 2.5 km. In contrast, the RF-PWV RMSE is less than 1 mm in all height difference segments. These findings demonstrate that RF-PWV offers improved correction effectiveness and higher accuracy compared with C-PWVC1. This enhanced adaptability to height differences enables a finer-scale description of the vertical distribution of PWV.

Figure 6 Accuracy of RF-PWV and C-PWVC1 in each height difference segment with respect to ERA5 data.

3.2 Validation of RF-PWV using RS PWV

To further validate the applicability of RF-PWV, the PWV data for all pressure levels within the 0–12 km altitude range from 148 RS stations in 2018 were used to assess the accuracy of RF-PWV and C-PWVC1. As the sounding stratified data are not uniformly distributed vertically, the variation in PWV with elevation was fitted using an exponential function based on the 2018 PWV data from each sounding station. Using the fitting results, the PWV data of neighboring levels were interpolated using inverse distance weighting (IDW) to generate a sequence of PWV values within the range of 0–12 km with intervals of 500 m. This sequence of PWV values served as reference values. The datum PWV is the PWV corresponding to the surface height of the RS station. For each RS station, the four grid points (1° × 1°) in proximity were selected, and the ΔPWV_i (i=1, 2, 3, 4, where i denotes the four proximity grid points) of the target height relative to the datum height was computed based on the RF model of each grid point. Then ΔPWV_i values were bilinearly interpolated to the corresponding location of the RS station to obtain the RF-PWV result. To account for systematic bias between modeling data and reference values data, the average difference between the corrected RF-PWV value and the corresponding reference value was computed as the systematic bias for each level at the RS station. Finally, the statistical accuracy of RF-PWV and C-PWVC1, after eliminating the systematic bias, is presented in Table 2.

Table 2 Validation results of the RF-PWV and C-PWVC1 models tested using RS data.

Table 2 reveals that the accuracy of C-PWVC1 is significantly lower than that of RF-PWV. C-PWVC1 exhibits a bias of 1.36 mm, ranging from −6.62 to 3.46 mm, whereas the RF-PWV bias is only 0.05 mm, reduced by 1.31 mm and improved by 96.36 %. The range of variation is notably reduced to −0.25 to 0.33 mm. Moreover, the RF-PWV RMSE is considerably smaller and more stable, with the RMSE reduced to 2.59 mm, ranging from 0.49 to 4.89 mm, corresponding to a decline rate of approximately 5 % compared with C-PWVC1. Consequently, RF-PWV demonstrates superior accuracy and stability in vertical PWV correction at 148 RS stations in the study area. Moreover, these results show that the accuracy analyzed by RS data is slightly lower than that estimated by ERA5 data. This is because of the significant systematic bias between ERA5 and RSs (Zhu et al., 2022; Sun et al., 2019), but such accuracy can still meet the meteorological requirements for PWV accuracy.

The bias and RMSE values for each RS station are also computed to further illustrate the application capabilities of the two models, as shown in Fig. 7. As depicted in Fig. 7a and c, C-PWVC1 exhibits a positive bias for almost all stations except for the RS stations in the Yunnan–Guizhou Plateau, where the bias is less pronounced. In contrast, the RF-PWV bias is consistently less than 0.5 mm and closer to 0 mm. Compared with C-PWVC1, the absolute value of the RF-PWV bias is effectively reduced in the Yunnan–Guizhou Plateau region, with the most significant reduction reaching 3.13 mm. Meanwhile, the positive bias in other areas is also reduced to varying degrees. Figure 7b and d demonstrate that the RF-PWV RMSE exhibits a certain degree of reduction compared with C-PWVC1, with the most substantial decline occurring in the sites located in the Yunnan–Guizhou Plateau. Given the complex terrain and significant undulations in the Yunnan–Guizhou Plateau, where the difference in height between the target point and the reference grid can be up to 1–2 km (Chen et al., 2011), the RF-PWV RMSE at all RS stations is less than 5 mm, with a maximum RMSE reduction of 11.65 mm. Therefore, RF-PWV demonstrates superior performance and more stable accuracy compared with C-PWVC1 across the entire study area. This advantage is particularly pronounced in regions with significant variations in height.

Figure 7 Distributions of the bias and RMSE for RF-PWV and C-PWVC1 with respect to the RS data.

The bias and RMSE values of RF-PWV and C-PWVC1 were also statistically analyzed for each month to assess the models' performance under different seasonal conditions. These results are presented in Fig. 8. The RF-PWV bias demonstrates improvement in every month compared with C-PWVC1. Both models exhibit seasonal variation characteristics, with lower accuracy during summer and higher precision in winter. This seasonal variation is attributed to the warm and humid weather with abundant rainfall in summer, leading to significant PWV fluctuations. Nevertheless, RF-PWV still shows notable bias optimization compared with C-PWVC1. Winters are typically drier and experience less rainfall, resulting in relatively smoother PWV changes. Consequently, both models can accurately capture PWV variations during this period, with RF-PWV having a distinct bias advantage. Furthermore, the RMSE of RF-PWV and C-PWVC1 exhibits similar variations to bias. While the RF-PWV RMSE is slightly larger than that of C-PWVC1 in late autumn and winter, it is smaller than C-PWVC1 in other months, particularly during summer and early autumn. RF-PWV's advantage becomes more pronounced when dealing with spatiotemporal PWV variations that are more drastic. It is important to note that differences between validation results based on radiosonde and ERA5 data may be attributed to certain systematic deviations between radiosonde and ERA5 data. In summary, RF-PWV demonstrates superior performance in vertical PWV correction under various seasonal conditions.

Figure 8 The bias and RMSE each month from RF-PWV and C-PWVC1.

To further evaluate the models' application in the vertical direction, the bias and RMSE of RF-PWV and C-PWVC1 in different height difference segments were examined, and the results are depicted in Fig. 9. In Fig. 9a and c, for C-PWVC1, when the height difference is less than 0 km, the bias and RMSE are −1.81 and 2.89 mm, respectively. As the height difference increases, both the bias and RMSE increase as well. When the height difference exceeds 2.5 km, the bias stabilizes at 2–2.5 mm, while the RMSE remains around 3.5 mm. RF-PWV demonstrates higher accuracy and stability across all height difference segments, with the bias approaching 0 mm and the RMSE being smaller than that of C-PWVC1. Figure 9b and d depict the improvement rates of the absolute values of bias and RMSE for RF-PWV compared with C-PWVC1 (positive values indicate improvement). The absolute value of the bias exhibits an improvement rate of over 80 %, with the maximum value approaching 100 %. Meanwhile, the improvement rate of the RMSE is significantly larger when the height difference is less than 3.5 km; it decreases slightly when the height difference exceeds 3.5 km but still remains around 5 %. In summary, RF-PWV offers higher vertical correction accuracy and improved stability across various height differences, demonstrating its strong applicability at different elevations.

Figure 9 Variation in the bias and RMSE with height differences (a, c) and improvement rates of the absolute values of bias and RMSE (b, d). It is noted that, in order to obtain PWV values that are uniformly distributed in the height direction, we obtained PWV values with heights in the range of 0–12 km; when the surface heights of some RS stations are greater than 0 km, their height differences relative to a height of 0 km are less than 0 km.

4 Conclusions and outlooks

Modeling accurate PWV vertical corrections benefits PWV fusion and provides detailed PWV vertical distribution information for meteorological studies. The complex terrain in China, characterized by varying climates and frequent water vapor exchanges, makes it challenging to accurately capture PWV variations at different heights. Consequently, this paper aims to develop a high-precision vertical PWV correction grid model. The primary contributions of this research can be summarized as follows:

1. We establish a PWV vertically corrected grid model (RF-PWV) with a resolution of 1° × 1° by integrating RF and monthly averaged hourly PWV data. This model utilizes RF to estimate the vertical variation in PWV at each grid point and demonstrates excellent applicability within a 6 km height difference. It effectively approximates PWV vertical changes. Validation against ERA5 data reveals that RF-PWV reduces the bias and RMSE by 99.84 % and 63.40 %, respectively, compared with C-PWVC1. RS validation also shows a reduction of 96.36 % for the bias and of 5 % for the RMSE compared with C-PWVC1. Furthermore, RF-PWV exhibits robust resistance to seasonal and height difference interference.

2. RF is employed to model each grid point (1° × 1°), with the grid serving to decompose spatial variations and confine RF within the corresponding grid point. This simplifies the features of training samples for each grid point RF, potentially reducing the likelihood of RF getting stuck in a local optimum. Simultaneously, during training, issues with a particular grid will not impact the models of other grid points, thereby enhancing modeling efficiency. This approach also eliminates concerns about spatial generalization ability and ensures relatively stable accuracy across all grid points, contributing to the model's robustness.

Comprehensive validation demonstrates that RF-PWV can more accurately provide PWV vertical corrections in China and its surrounding areas. This model holds great potential for PWV vertical correction and is well suited for delivering detailed PWV vertical distribution information for multisource water vapor fusion and meteorological research. Consequently, this method can be used to develop a globally applicable vertical correction model with higher accuracy, benefiting a wider range of users.