Spatial–Temporal Relationship Study between NWP PWV and Precipitation: A Case Study of ‘July 20’ Heavy Rainstorm in Zhengzhou

Abstract

In order to study and forecast extreme weather, a comprehensive and systematic analysis of the spatial and temporal relationship between Precipitable Water Vapor (PWV), predicted by Numerical Weather Predication (NWP) data, and precipitation, is necessary. The goal of this paper was to study the temporal and spatial relationship between PWV and precipitation during the so-called ‘July 20’ (18–21 July 2021) heavy rainstorm in Zhengzhou. Firstly, the PWV data provided by 120 radiosonde stations uniformly distributed throughout the world, and two IGS stations in China, in 2020, was used to evaluate the accuracy of PWV estimation by ERA5 and MERRA-2 data, and the factors affecting the accuracy of NWP PWV were explored. Secondly, ERA5 PWV and the precipitation data of six meteorological stations were used to qualitatively analyze the relationship between PWV and precipitation during the ‘July 20’ heavy rainstorm in Zhengzhou. Finally, a quantitative study was conducted by an eigenvalue matching method. The main experimental results were as follows. Compared with MERRA-2 PWV, the accuracy of ERA5 PWV was slightly higher. Latitude, altitude and season were the influencing factors of the NWP PWV estimation accuracy. The change trend of ERA5 PWV was consistent with both 24 h cumulative precipitation and surface precipitation during the ‘July 20’ heavy rainstorm in Zhengzhou. The average optimal matching degree and optimal matching time between NWP PWV and surface precipitation during the ‘July 20’ heavy rainstorm in Zhengzhou was 56.6% and 3.68 h, respectively. The maximum optimal matching degree was 80.3%. The spatial–temporal relationship between NWP PWV and surface precipitation was strong.

1. Introduction

The troposphere, as the atmosphere layer most closely related to human activities, contains almost all water vapor and 75% of the air quality of the whole atmosphere, which is an important part of global space. Although the water vapor content in the atmosphere is only 0–4%, its change has a direct impact on surface temperature, humidity and surface precipitation [1,2,3]. In addition, water vapor is an important greenhouse gas, and about 60% of greenhouse gases that cause climate change are affected by water vapor [4]. At present, global climate change is irregular: drought, flood disasters and extreme weather occur frequently, which is having a serious impact on human production and life. Accurately and systematically monitoring the distribution of atmospheric water vapor, and exploring its change rule, play an important role in studying the evolutionary processes of varieties of complex weather, and providing early warning of natural disasters [5].

Precipitable Water Vapor (PWV) is the most commonly used indicator to express the water vapor content in the atmosphere, which refers to the total water vapor content contained in the air column of a unit cross-sectional area from the surface to the top of the troposphere [6]. Traditional PWV observation methods include radiosonde (RS), microwave radiometer and satellite remote sensing [7,8]. Of the traditional methods, RS is one of the most accurate PWV detection methods at present, but its time resolution is low, the site distribution is sparse, and the detection cost is high, which cannot meet the needs of small-scale and medium-scale meteorological research [9]. Microwave radiometer and satellite remote sensing are mainly based on thermal infrared and near-infrared bands, to detect meteorological parameters in the atmosphere, which have high accuracy and large detection range. However, correction is required before use. The detection accuracy is easily affected by weather conditions, and the vertical profile of the water vapor cannot be obtained, which restricts its application in weather forecasting and meteorological research [10,11]. With the development of the Global Navigation Satellite System (GNSS), the concept of GNSS meteorology was first proposed in the 1990s [12]; the PWV monitoring method based on GNSS technology has been widely used since then. This method has the merits of high resolution, low cost, continuous operation and high precision. However, due to the uneven distribution of GNSS stations, it is difficult to obtain sufficient GNSS PWV data for areas with few GNSS stations, which limits the application of GNSS technology to invert PWV in climate research. With the increasing accuracy of reanalysis data, the estimation of PWV using Numerical Weather Prediction (NWP) forecasts or reanalysis data, and the study of the relationship between PWV and extreme weather, have become current hot topics [13].

At present, organizations including the European Centre for Medium-Range Weather Forecasts (ECMWF), the National Centers for Environmental Prediction (NCEP), the National Aeronautics and Space Administration (NASA), and the China Meteorological Administration (CMA), provide their users with the latest global atmospheric numerical prediction reanalysis information and forecast information. The provision, by these organizations, of accurate assessment of PWV predicted by the NWP data, is necessary. Bock, et al. [14], compared PWV calculations from 120 IGS stations worldwide with PWV calculated from ERA-Interim reanalysis data provided by ECMWF: the results showed that the daily standard deviation of PWV was usually less than 2 mm. Zhang, et al. [15], used GNSS PWV to evaluate the accuracy of PWV estimations retrieved from ERA5 and ERA-interim reanalysis data provided by ECMWF in China: the RMS values of PWV retrieved from ERA5 and ERA-interim were 1.8 mm and 2.1 mm, respectively. Vey, et al. [16], showed good consistency between GNSS PWV and PWV estimations retrieved from NCEP reanalysis data; however, in the Antarctic and tropical regions, PWV retrieved from NCEP reanalysis data were smaller than those retrieved from GNSS PWV. Chen, et al. [17], used the GNSS PWV to evaluate PWV retrieved from the Climate Forecast System Reanalysis (CFSR) dataset provided by NCEP, with an RMS value of 4.13 mm. Huang, et al. [18], used GNSS PWV to evaluate the accuracy of ERA5 PWV provided by ECMWF, and MERRA-2 PWV provided by NASA, on the Tibetan Plateau, with RMS values of 1.77 mm and 2.12 mm, respectively. The CRA40 reanalysis data, released by the CMA, showed a temporal resolution of 3 h, a spatial resolution of 43 km, and data starting time of 1979 [19]. There is scant literature evaluating the accuracy of the estimation of PWV from CRA40 reanalysis data.

Many scholars have studied the relationship between PWV and extreme weather. The occurrence of rainfall requires sufficient water vapor content in the atmosphere, and there is constant condensation of water vapor in the atmosphere a few hours, or sometimes even more, before rainfall occurs, resulting in a sustained growth trend in PWV values [20]. Duc, et al. [21], used the extended fractions skill score method to evaluate the performance of MF10 km and MF2 km ensemble forecast data for predicting rainfall. Experiments showed that MF2 km was more reliable than MF10 km in predicting moderate rainfall and rainstorms. In contrast, the MF10 km data was superior to the MF2 km data in predicting small rainfall. Champollion, et al. [22], studied the variation and distribution of PWV during a rainstorm in southwest France on 9 September 2002: the results showed that small-scale variation of PWV can also cause deep convection phenomena. The necessary condition for rainfall is the increase of PWV, but the increase of PWV does not necessarily lead to rainfall events. Yao, et al. [23], pointed out that PWV surged 2–6 h before the rainstorm. Heavy precipitation events occur after a steep rise in PWV [24]. Zhu, et al. [25], studied the variation of PWV during three typhoons in Hong Kong from 2013 to 2014, and the experiment showed that PWV near the surface (0~1.6 km) varied less during the typhoon, while PWV between 1.6~8.5 km in height varied more during the typhoon. Valjarevi, et al. [26], used cloud data from moderate resolution imaging spectroradiometer satellite for 30 years (1989–2019), to analyze the relationship between cloud cover and topography in Serbia. The results showed that the lower the cloud cover, the lower the precipitation in the eastern mountainous areas of Serbia.

Henan Province is located in central China, and its capital is Zhengzhou city. The terrain is high in the west and low in the east. The north, west, south are surrounded by Taihang Mountain, Funiu Mountain, Tongbai Mountain and Dabie Mountain, respectively. The central and eastern part is the Huang-Huai-Hai Plain. Affected by typhoon, topography and atmospheric circulation, a rare continuous heavy rainfall occurred in Zhengzhou from 18 to 21 July 2021. The whole city generally suffered heavy and extraordinarily heavy rains, and the cumulative average precipitation was 449 mm. This was the so-called ‘July 20’ heavy rainstorm in Zhengzhou. However, no research on PWV in the Zhengzhou rainstorm has been published in English.

To sum up, MERRA-2 and ERA5 are the latest generation of NWP datasets, and a comprehensive evaluation of the accuracy of the PWV estimation of both datasets is necessary for studying the PWV variation pattern [27,28]. On the other hand, the existing research on the ‘July 20’ heavy rainstorm in Zhengzhou mainly focuses on the temporal relationship between PWV and precipitation, and there are few studies on the temporal–spatial relationship between PWV and precipitation, and the change process of precipitation trajectories. A comprehensive and detailed study on the relationship between PWV and precipitation during the rainstorm will be of great significance for improving the early-warning ability of rainstorm prediction. Therefore, in this paper, the accuracy of PWV estimation from MERRA-2 and ERA5 reanalysis data was evaluated by using PWV data from 120 radiosonde stations uniformly distributed throughout the world, GNSS PWV data calculated by 2 IGS stations in China, and analysis of the influencing factors of PWV estimation accuracy from reanalysis data. Taking the ‘July 20’ heavy rainstorm in Zhengzhou in 2021 as an example, the change process of the precipitation trajectory during a rainstorm was qualitatively studied, based on the precipitation data of meteorological stations and ERA5 PWV, and the eigenvalue matching method was proposed, to quantitatively explore the spatial–temporal relationship between NWP PWV and surface precipitation.

The experimental data, the PWV calculation method and the accuracy index are described in the second section of this paper. The third section is the experimental results and analysis, which mainly includes the evaluation of the accuracy of NWP-estimated PWV, and the analysis of the spatial and temporal distribution relationship between PWV and precipitation during the ‘July 20’ heavy rainstorm in Zhengzhou. The fourth section summarizes the work of this paper.

2. Data and Methodology

2.1. Data Description

In order to analyze the accuracy and error distribution characteristics of PWV retrieved from MERRA-2 reanalysis data and ERA5 reanalysis data in the global area, and to study the spatial–temporal relationship between PWV and precipitation in Henan Province during the ‘July 20’ heavy rainstorm in Zhengzhou in 2021, this paper collected the PWV data from the RS, MERRA-2 reanalysis dataset, the ERA5 reanalysis dataset, the tropospheric delay products of the International GNSS Service (IGS) and the precipitation data of six meteorological stations around Zhengzhou city.

2.1.1. MERRA-2 Data

MERRA-2 data is published by NASA, and started in 1980 [27,29]. This paper collected the MERRA-2 M2I6NPANA dataset at UTC 12 every day from January 1 to 31 December 2020. The MERRA-2 reanalysis data type was griding, the time resolution was 6 h, the horizontal resolution was 0.625° × 0.5°, the vertical level was 42 layers, and the data download address was https://disc.gsfc.nasa.gov (accessed on 5 December 2021).

2.1.2. ERA5 Data

The ERA5 pressure-level data and single-level data were derived from the ERA5 reanalysis data provided by the ECMWF [30]. The ERA5 pressure-level data at UTC 12 daily in 2020 was experimentally collected for PWV accuracy evaluation, and the surface precipitation data was used to analyze the spatial–temporal relationship between PWV and surface precipitation. Both the ERA5 pressure-level data and the surface precipitation data were from the ERA5 dataset, and the horizontal resolution of this dataset was 0.625° × 0.5°. The ERA5 reanalysis data horizontal resolution was 0.25° × 0.25°, the time resolution was 1 h, and the vertical resolution was divided into 37 layers. The data download address was https://www.ecmwf.int (accessed on 26 November 2021).

Table 1. Reanalysis data information.

Dataset	Agency	Maximum Time Resolution	Maximum Horizontal Resolution	Vertical Resolution	Assimilation Method
CRA40	CMA	6 h	0.3125° × 0.3125°	47	4DVAR
MERRA-2	NASA	6 h	0.625° × 0.5°	42	GEOS-5
ERA5	ECMWF	1 h	0.25° × 0.25°	37	4DVAR

shows the CRA40, MERRA-2 and ERA5 reanalysis data information tables. It can be seen from Table 1 that ERA5 had the highest maximum spatio–temporal resolution. As CRA40 data could not be download freely, PWV estimated by MERAA-2 and ERA5 data was assessed.

2.1.3. RS Data, GNSS ZTD Data and Precipitation Data

The RS station data was collected from the Integrated Global Radiosonde Archive Version 2 (IGRA2) dataset generated by the National Climate Data Center. The location distribution of the selected 120 stations is shown in Figure 1. The experiment used the RS data of UTC 12 each day in 2020. The IGRA2 dataset included pressure, temperature, relative humidity, potential height, wind speed and other data. The time resolution of the data was 12 or 6 h. The data download address was ftp://ftp.ncdc.noaa.gov/pub/data/igra [31] (accessed on 5 November 2021).

In this paper, GNSS PWV estimations, recorded at UTC 12 every day between 1 January and 31 December, 2020, were used as a reference [32]. The selected IGS station location is shown in Figure 1.

The precipitation data of the meteorological stations from July 2021 were provided by the Meteorological Information Center of the China Meteorological Administration. The time resolution of this data was 1 h.

2.2. Methodology

In this paper, the PWV data provided by the RS stations and GNSS PWV were regarded as true values for evaluating the accuracy of NWP PWV. RMSE and Bias were used as precision evaluation indexes.

2.2.1. NWP PWV Estimation Method

The PWV estimation formula, using NWP reanalysis data, was [15]:

$$\mathrm{q}=\frac{0.622e}{p-0.378e}$$

(1)

$$PWV=\frac{1}{\mathrm{g}}{\displaystyle {\int }_0^{p_s}qdp}\approx \frac{1}{2\mathrm{g}}{\displaystyle \sum _{i=1}^n\left(q_i+q_{i+1}\right)·\left(p_i-p_{i+1}\right)}$$

(2)

where $q$ was specific humidity, $e$ was water vapor pressure, $p$ was atmospheric pressure, $p_s$ was surface pressure and $\mathrm{g}$ was a gravitational parameter. As the $\mathrm{g}$ values in different regions were different, Equation (3) was used to calculate the gravitational parameter in different regions in this paper. The formula was as follows ($\phi$, $H$ was the latitude and geodetic height of the station [33]):

$$\mathrm{g}\left(\phi ,H\right)=9.80616\left(1-2.59\times {10}^{-3}\cos 2\phi \right)·\left(1-3.14\times {10}^{-7}H\right)$$

(3)

The elevation of MERRA-2 and ERA5 reanalysis data was geopotential height, the elevation of RS data was orthometric height, and the elevation of IGS data was geodetic height. In order to eliminate the elevation error caused by different elevation systems, different elevations were converted to orthometric height. The geopotential height $H_{geo}$ calculate formula, using geopotential, was:

$$H_{geo}=\frac{Geo}{\mathrm{g}}$$

(4)

The $Geo$ was geopotential and $\mathrm{g}$ was a gravitational parameter.

To resolve the problem that the plane position and the elevation of the NWP reanalysis dataset grid points did not coincide with those of the IGRA2 RS stations, the PWV value of four grid points around the sounding station was used to interpolate the PWV value at the RS station. In order to weaken the errors caused by the multiple interpolations, this paper first calculated the PWV value of the grid points at the same elevation as the site, and then used the Kriging interpolation method for plane interpolation [34]. Since the reanalysis dataset did not provide meteorological data at the surface, Formulas (5) and (6) were used, in this paper, to extrapolate the pressure and other parameters [33], taking the pressure calculation at station height $H$ as an example:

$$P_H=P_{lower}\exp \left(-\frac{H-H_{lower}}{h_P}\right)$$

(5)

$$h_{\mathrm{p}}=\frac{H_{upper}-H_{lower}}{\ln \left(P_{lower}/P_{upper}\right)}$$

(6)

where, $H_{upper}$, $H_{lower}$, $P_{lower}$ and $P_{upper}$ were the upper and lower geopotential heights, and Pressure, $P_H$, was the pressure at height $H$.

2.2.2. GNSS PWV Calculation Method

The GNSS PWV calculation method was as follows: (1) the Zenith Hydrostatic Delay (ZHD) was calculated using the Saastamoinen model [35]; (2) the Zenith Wet Delay (ZWD) was stripped from ZTD provided by IGS using Equation (7); (3) the GNSS PWV was calculated by ZWD, using Equation (8) [36].

$$ZWD=ZTD-ZHD=ZTD-\frac{0.002277·P}{1-0.00266·\cos \left(2\phi \right)-0.00028·H}$$

(7)

$$PWV=\frac{{10}^6}{\left(K_2^{’}+K_3/T_{\mathrm{m}}\right)·R_V·\rho }·ZWD$$

(8)

where $H$, $P$ and $\phi$ were, respectively, the geodetic height, pressure and latitude at the station; $K_2^{’}$, $K_3$, and $R_V$ were constants, and their values were 16.48 K∙hPa⁻¹, (3.776 ± 0.014) × 105 K²∙hPa⁻¹ and 461∙J∙(Kg∙K)⁻¹, respectively. $\rho$ was the density of liquid water, and its value was 10³ Kg·m⁻³. $T_{\mathrm{m}}$ was the atmospheric weighted mean temperature, which was obtained by the Bevis formula, $T_{\mathrm{m}}$ = 70.2 + 0.72T. T was the atmospheric temperature of the station, in K. Research showed that the RMSE of $T_{\mathrm{m}}$, calculated by the Bevis formula in the middle latitude region, was 4.74 K [12]. P and T were provided by MERRA-2 or ERA5 reanalysis data.

2.2.3. Gross Error Detection of RS PWV

As the sounding balloon was easily affected by multiple factors in the measurement, there were outliers in the measured values. As a consequence, the RS PWV estimation needed to be pre-processed. In this paper, the Interquartile Range (IQR) method was used to eliminate the deviation in the RS PWV data. The IQR criterion assumed that the detected target sequence conformed to the standard normal distribution. If the target sequence was arranged from small to large, the 25th percentile was called the lower quartile, and the 75th percentile was called the upper quartile. The IQR value was the difference between the upper quartile and the lower quartile. The anomaly detection interval based on the IQR was [37]:

$$\left[Q1-1.5·IQR,Q3+1.5·IQR\right]$$

(9)

where, $Q1$ and $Q3$ denoted the lower quartile and the upper quartile, respectively, and the IQR was IQR = $Q3$ − $Q1$.

The specific steps of the IQR gross error detection method were as follows: (1) the periodic model of each station was established by Formula (10), and the coefficients of the periodic model were fitted by RS PWV data; (2) the PWV model value was calculated according to the periodic model, and the difference between the PWV model value and the RS PWV was calculated to obtain the residual sequence of each station; (3) the RS PWV that the residual offset out of the given gross error detection interval was removed, and replaced by the model value [38].

$$\begin{array}{r}PWV=A_0+A_1\cos \left(\frac{doy}{365.25}2\pi \right)+B_1\sin \left(\frac{doy}{365.25}2\pi \right)\\ +A_2\cos \left(\frac{doy}{365.25}4\pi \right)+B_2\sin \left(\frac{doy}{365.25}4\pi \right)\end{array}$$

(10)

where: $A_0$ was the mean value of $PWV$; $A_1$, $A_2$, $B_1$ and $B_2$ were the required parameters; $doy$ was the Day of the Year.

2.2.4. Eigenvalue Matching Method

In order to quantitatively analyze the relationship between PWV and surface precipitation during rainstorms in Henan Province, an eigenvalue matching method based on the upper quartile algorithm was proposed. At present, there are many different ways to determine the similarity level between spatial distribution images [37,39,40,41]. In this paper, an eigenvalue matching method was used. In the IQR criterion, the value at the 75th percentile was called the ‘larger quartile’, also known as the upper quartile [37]. In the proposed eigenvalue matching method, the upper quartile was used as the threshold value, and PWV and surface precipitation values greater than their upper quartiles were taken as eigenvalues. The selected experimental area was 110°–117.5°E and 30°–38°N, which contained the Henan Province and surrounding areas. The flow chart of the eigenvalue matching method is shown in Figure 2. Firstly, taking the UTC 00 18 July as an example, the PWV value of the selected area at this hour was calculated, and the PWV value was arranged in descending order to determine the upper quartile value. The PWV value greater than the upper quartile was used as the PWV eigenvalue of the selected area. Secondly, as the rainfall mainly occurred within 2–6 h after PWV climbed [24], the surface precipitation eigenvalue of the selected area at 2–6 h after UTC 00 was determined using the same method. Thirdly, studies have shown that the greater the PWV, the greater the rainfall intensity [20]. This indicated that PWV at this grid point was strongly correlated with surface precipitation when the grid points of the PWV eigenvalue and the surface precipitation eigenvalue were the same. In this step, the grid points with both the surface precipitation eigenvalues at each hour and the PWV eigenvalues were selected, and were named the matching grid points. These points expressed the spatial relationship between PWV and surface precipitation. Fourthly, the matching degree at each hour was calculated using Equation (11). The maximum matching degree at 2–6 h after UTC 00 was named as the optimal matching degree (OMD) between PWV and surface precipitation, and the hour with the OMD was called the optimal matching time (OMT). The surface precipitation at the MOT had the strongest correlation with the current PWV, and this OMT expressed the temporal relationship between PWV and surface precipitation.

$$MD=\frac{CGP}{GP}\times 100\%$$

(11)

where $CGP$ was the number of matched grid points, $GP$ was the number of PVW eigenvalues at the selected hour, and $MD$ was the matching degree.

2.3. Precision Evaluation Index

Bias and RMSE values were introduced to evaluate the accuracy of PWV estimated from the MERRA-2 and ERA5 reanalysis data. The expressions of Bias and RMSE were as follows:

$$B\mathrm{ias}=\frac{{\displaystyle {\sum }_{i=1}^N\left(X_{r\mathrm{e},i}-X_{r,i}\right)}}{N}$$

(12)

$$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^N{\left(X_{r\mathrm{e},i}-X_{r,i}\right)}^2}}{N}}$$

(13)

where, $N$ denoted the total number of samples, $X_{r\mathrm{e},i}$ was PWV calculated for reanalysis data, and $X_{r,i}$ was RS PWV or GNSS PWV.

3. Results and Analysis

3.1. Accuracy Evaluation and Accuracy Affecting Factor Analysis of NWP PWV

In order to compare the accuracy of PWV inversed by the MERRA-2 and ERA5 data, and to analyze the factors affecting the estimation accuracy of NWP PWV, MERRA-2 PWV and REA5 PWV from 1 January to 31 December 2020 were calculated, and the PWV calculated by the 120 RS stations and 2 IGS stations was used as the reference. The correlation between RS PWV and NWP PWV, Bias and RMSE of NWP PWV in the global region, and the relationship between latitude, altitude, season and NWP PWV accuracy is discussed in this section.

3.1.1. NWP PWV Accuracy Evaluation

Due to space limitation, this paper randomly selected an RS station to show its correlation between MERRA-2 PWV and RS PWV, and between ERA5 PWV and RS PWV, as shown in Figure 3. The selected station, CAM00071802, was located in the middle latitudes of the northern hemisphere. It can be seen from this figure that the correlation coefficients between MERRA-2 PWV and RS PWV, and between ERA5 PWV and RS PWV, were 0.9838 and 0.9891, respectively.

Table 2. The correlation coefficient of MERRA-2 PWV, ERA5 PWV and RS PWV at portion stations.

Station Name	Latitude	Longitude	Altitude (m)	Coefficient between MERRA-2 PWV and RS PWV	Coefficient between ERA5 PWV and RS PWV
GLM00004360	65.61°N	37.63°W	54.0	0.8451	0.9181
CAM00071802	47.51°N	52.78°W	112.4	0.9679	0.9784
CHM00058027	34.28°N	117.15°E	42.0	0.9751	0.9822
SAM00040430	24.55°N	39.70°E	654.0	0.7224	0.7326
IDM00097014	1.53°N	124.91°E	80.0	0.6434	0.7365
IDM00097072	0.68°S	119.73°E	6.0	0.5516	0.6223
IDM00097560	1.18°S	136.11°E	11.0	0.5820	0.6209
IDM00097180	5.06°S	119.55°E	14.0	0.4539	0.5477
ASM00094403	28.80°S	114.69°E	36.9	0.7354	0.7373
NZM00093844	46.41°S	168.31°E	2.0	0.7123	0.7428

shows the correlation coefficient table of MERRA-2 PWV, ERA5 PWV and RS PWV at portion stations selected randomly. It can be seen from Table 2 that the correlation of stations in Southeast Asia was lower than that in other regions.

The annual average Bias distribution of MERRA-2 PWV and ERA5 PWV in 2020, referred to RS PWV at 120 RS stations, is shown in Figure 4. From Figure 4, we can see that the Bias of MERRA-2 PWV and ERA5 PWV ranged between −2 mm to 5 mm. The average MERRA-2 PWV Bias of the 120 RS stations was 1.26 mm, which was 1.15 mm for ERA PWV. The Bias of MERRA-2 PWV was similar to that of ERA5 PWV. The average Bias of MERRA-2 PWV and ERA5 PWV in the northern hemisphere was 1.22 mm and 1.04 mm, respectively. The average Bias of MERRA-2 PWV and ERA5 PWV in the southern hemisphere was 1.35 mm and 1.44 mm, respectively. The PWV accuracy of these two data in the northern hemisphere was slightly higher than that in the southern hemisphere.

The annual RMSE distribution of MERRA-2 PWV and ERA5 PWV in 2020, referred to RS PWV at the 120 RS stations, is shown in Figure 5. It can be seen from Figure 5 that the RMSE of PWV estimated from these two data was less than 4 mm, except for individual stations in the equatorial region of Southeast Asia. The average RMSE of MERRA-2 PWV and ERA5 PWV was 3.76 mm and 3.20 mm, respectively. The average RMSE of MERRA-2 PWV in the northern hemisphere was 3.66 mm, and that of the southern hemisphere was 4.04 mm. The average RMSE of the ERA5 PWV in the northern hemisphere and the southern hemisphere was 3.07 mm and 3.53 mm, respectively. The RMSE of MERRA-2 PWV was higher than that of ERA5 PWV in the southern and northern hemispheres. In short, the accuracy of ERA5 PWV was slightly higher than that of MERRA-2 PWV referred to RS PWV.

It can be seen from Figure 4 and Figure 5 that in the equatorial region of Southeast Asia, the accuracy of PWV estimated from NWP reanalysis data was slightly lower. The equatorial region of Southeast Asia is a tropical rainforest climate. Affected by typhoons and other extreme weather, the atmospheric water vapor is abnormally active, resulting in low accuracy of NWP data [42]. As the NWP reanalysis dataset assimilated data from different data sources during its generation [43], there were more ground observation stations in the northern hemisphere, and the data sources were richer than those in the southern hemisphere (such as RS and surface meteorological station observations). Therefore, the accuracy of PWV estimated from the two NWP reanalysis data in the northern hemisphere was higher than that in the southern hemisphere.

In order to explore the accuracy of NWP PWV referred to GNSS PWV, PWV at the IGS station of URUM and JFNG was retrieved from the reanalysis data of MERRA-2 and ERA5 in 2020. Taking the GNSS PWV from IGS as reference, the Bias distribution of MERRA-2 PWV and ERA5 PWV is shown in Figure 6. The average Bias and RMSE of MERRA-2 PWV and ERA5 PWV at URUM station were 1.3 mm and 0.54 mm, and 2.82 mm and 2.31 mm, respectively. The average Bias and RMSE of MERRA-2 PWV and ERA5 PWV at JFNG station were 1.9 mm and 1.03 mm, and 3.79 mm and 2.43 mm respectively. The accuracy of PWV calculated by ERA5 data was higher than that of MERRA-2 referred to GNSS PWV. In addition, the PWV Bias in summer was higher than that in winter.

3.1.2. NWP PWV Accuracy Affecting Factor Analysis

To investigate the factors affecting the accuracy of NWP PWV, the accuracy of NWP PWV was studied from the aspects of latitude, altitude and season. The relationship between the Bias/RMSE and the latitude of the 120 RS stations is shown in Figure 7. The negative value of latitude represents the south latitude, and the positive value of latitude represents the north latitude. It can be seen from the figure that the accuracy of PWV estimation from the NWP reanalysis data showed a high correlation with latitude, and its accuracy increased with the increase of latitude. In the equatorial and low latitude regions, the accuracy of PWV was relatively low. The reason was that the atmospheric water vapor content changed rapidly due to the complex meteorological conditions in low latitude regions. In short, latitude was one of the main factors affecting the accuracy of PWV estimation from NWP reanalysis data.

The relationship between the PWV Bias/RMSE and the altitude of the 120 RS stations is shown in Figure 8. In Figure 8, the horizontal axis represents height, and the vertical axis represents Bias and RMSE, respectively. It can be seen from the figure that the PWV Bias/RMSE reduced with the increase of the altitude of the station. At high altitude, the water vapor content was relatively small, and the climatic conditions were stable. Therefore, the accuracy of the NWP PWV at high altitude was better than that at low altitude. In short, altitude was also one of the main factors affecting the accuracy of PWV estimation from NWP reanalysis data.

The correlation between the monthly mean Bias value (the red circle) of NWP PWV and the seasons is shown in Figure 9. CAM00071802 station and NZM00093844 station are located in the northern hemisphere and southern hemisphere, respectively. In summer, the PWV estimation accuracy of the two reanalysis data was lower than in the other seasons. This was related to the fact that the atmosphere in summer is more active than in other seasons, and the distribution and variation of water vapor are more complex. The trends of ERA5 PWV and MERRA-2 PWV were similar to one another. In short, ‘the season’ was one of the factors affecting the accuracy of PWV estimated from reanalysis data.

3.2. Spatial–Temporal Relationship between NWP PWV and Precipitation during the ‘July 20’ Heavy Rainstorm in Zhengzhou in 2021

A comprehensive and detailed study of the relationship between PWV and precipitation in the ‘July 20’ heavy rainstorm in Zhengzhou is of great significance to further improving early-warning ability for extreme weather, and reducing the loss of people’s lives and property. Therefore, this paper qualitatively and quantitatively analyzed the change process of the precipitation trajectory, and the spatial–temporal relationship between PWV and surface precipitation during the ‘July 20’ heavy rainstorm in Zhengzhou. Firstly, the variation process of precipitation trajectories during rainstorms, and the relationship between PWV and precipitation were qualitatively analyzed, based on the PWV and hourly precipitation of six meteorological stations (Figure 10) around Zhengzhou for each July from 2019 to 2021. Secondly, an eigenvalue matching method was proposed, to quantitatively study the spatial–temporal distribution relationship between surface precipitation and PWV during the ‘July 20’ heavy rainstorm in Zhengzhou in 2021. As the accuracy of ERA5 PWV was slightly higher than that of MERRA-2 PWV, according to the previous experiment, the ERA5 data was used for PWV calculation in subsequent experiments. As the heavy rain is mainly concentrated in the central and northern regions of Henan province, the meteorological stations of these regions were chosen. The distribution of the six selected meteorological stations is shown in Figure 10. The red circle represents the location of the meteorological station, and the blue line represents the river.

3.2.1. Qualitative Analysis of the Spatial–Temporal Relationship between NWP PWV and Precipitation

The average PWV of the six meteorological stations each July, from 2019 to 2021, is shown in

Table 3. Average PWV in July 2019 to 2021.

Meteorological Station	2019 (mm)	2020 (mm)	2021 (mm)
PingDingShan Station	43.36	44.35	47.22
JiaoZuo Station	49.62	50.77	53.48
ZhengZhou Station	44.92	45.98	48.88
XuChang Station	48.72	50.17	52.67
KaiFeng Station	49.30	51.87	54.13
XinXiang Station	49.76	51.44	53.82

. The trends of PWV in July 2019 to 2021, and precipitation per hour in July 2021, are shown in Figure 11. From Table 3 and Figure 11, we can see that: (1) the average PWV of the six meteorological stations in July 2021 was higher than that in 2020 and 2019; (2) during the rainstorm period (18–21 July 2021), each day’s PWV of six meteorological stations was higher than that in the previous three years; (3) precipitation was positively correlated with PWV, and there was a stage of PWV cumulative growth before each rainfall event (cyan ellipses, for instance).

The 24 h cumulative precipitation of the meteorological stations during heavy rainfall is shown in

Table 4. 24 h cumulative precipitation during rainstorms (the red values represent values of daily precipitation greater than 50 mm).

Meteorological Station	18 July 2021 (mm)	19 July 2021 (mm)	20 July 2021 (mm)	21 July 2021 (mm)
PingDingShan Station	55.0	209.5	20.2	5.3
JiaoZuo Station	59.5	66.4	209.2	235.8
ZhengZhou Station	37.0	228.0	376.3	77.4
XuChang Station	3.6	166.2	175.2	27.9
KaiFeng Station	0.3	63.7	83.1	22.1
XinXiang Station	34.1	42.5	242.7	258.6

, and the surface precipitation distribution in Henan and the surrounding areas at UTC 00, 06, 12 and18, from 18 to 21 July 2021, is shown in Figure 12. According to the definition of ‘rainstorm’ by the CMA, ‘rainstorm’ refers to 24 h cumulative precipitation greater than 50 mm. The World Meteorological Organization (WMO) divides precipitation with 1 h of precipitation greater than 10 mm into rainstorm grades. From Table 4, we can see that the rainstorm first occurred in the western mountainous area (Pingdingshan and Jiaozuo Station, Pingdinghsan/Jiaozuo, China), and then the precipitation trajectory gradually moved to the central region of Henan province (Zhengzhou and Xuchang Station, Zhengzhou/Xuchang, China). Finally, the rainstorm moved to the eastern area (Kaifeng Station, Kaifeng, China) and the north area (Xinxiang Station, Xinxiang, China). A similar trend of surface precipitation can also be seen in Figure 12. Taking the second column as an example, at the beginning of the rainstorm, the surface precipitation in the western and near-alpine regions was greater than that in other regions, and then gradually shifted from south to north and from west to east. The PWV distribution in Henan and the surrounding areas at UTC 00, 06, 12 and 18, from 18 to 21 July 2021, is shown in Figure 13. During the rainstorm period, the orange patches in Figure 13 (PWV bigger than 65 mm) moved from west to east and from south to north in Henan, and then dissipated after the rainstorm ended. This trend of PWV was similar to that of the 24 h cumulative precipitation and surface precipitation.

3.2.2. Quantitative Analysis of the Relationship between NWP PWV and Precipitation by Eigenvalue Matching Method

The OMT and OMD between PWV and surface precipitation during rainstorms was calculated every 6 h, and saved in

Table 5. Optimal matching degree and optimal matching time. TP refers to surface precipitation.

	Time of PWV	Time of TP	OMT	OMD		Time of PWV	Time of TP	OMT	OMD
18 July	00:00	02:00	2 h	31.0%	20 July	00:00	03:00	3 h	68.8%
	06:00	08:00	2 h	36.1%		06:00	08:00	2 h	54.1%
	12:00	18:00	6 h	56.0%		12:00	15:00	3 h	55.7%
	18:00	23:00	5 h	47.5%		18:00	00:00	6 h	75.4%
19 July	00:00	02:00	2 h	59.0%	21 July	00:00	02:00	2 h	54.0%
	06:00	10:00	4 h	65.6%		06:00	12:00	6 h	42.6%
	12:00	17:00	5 h	80.3%		12:00	15:00	3 h	50.5%
	18:00	22:00	4 h	60.6%		18:00	22:00	4 h	68.9%

. Due to space limitation, the PWV eigenvalue grid point at (a) UTC 12 on 19 July 2021 and (b) UTC 18 on 20 July 2021, and the corresponding matching grid point distribution, are shown in Figure 14. From this table and figure, we can see that the OMD between PWV and surface precipitation ranged between 31.0% and 80.3%, and the average OMD was 56.6%. Twelve groups of PWV and surface precipitation had OMD greater than 50%. The average OMD values corresponding to the OMT values of 2, 3, 4, 5 and 6 h were 46.8%, 58.3%, 65.0%, 63.9% and 58.0%, respectively. The OMT ranged between 2 and 6 h after the time of PWV, and 2 h was the most frequent OMT. The average OMT was about 3.68 h. The spatial–temporal relationship between PWV and surface precipitation was strong.

4. Discussion

Comprehensive and systematic monitoring of PWV changes is necessary for studying and forecasting extreme weather. In this paper, the accuracy of MERRA-2 PWV and ERA5 PWV were assessed by using RS PWV and GNSS PWV. The experimental results showed that, compared with RS PWV, the average RMSE of PWV estimated by MERRA-2 and ERA5 was 3.76 mm and 3.20 mm, respectively. Compared with GNSS PWV, the average RMSE was 3.305 mm and 2.37 mm, respectively. Consequently, the accuracy of the ERA5 PWV was slightly higher. This result was consistent with Huang’s research [18]. It can be seen from Table 2, Figure 4 and Figure 5 that NWP PWV had little correlation with RS PWV in the equatorial region of Southeast Asia, and low accuracy. The equatorial region of Southeast Asia is a tropical rainforest climate, and the atmospheric water vapor is extremely active, resulting in the low accuracy of NWP data here [42].

From Table 3 and Figure 11, it can be seen that there was a continuous growth stage of PWV before each rainfall. The more intense the PWV changes, the higher the probability of heavy rainfall. This was consistent with Zhao’s research [20]. Figure 14 and Table 5 show that the spatial distribution of PWV was consistent with the spatial distribution of precipitation, and that the average OMD and OMT values between PWV and precipitation were 56.63% and 3.68 h, respectively. Barindelli’ s research showed that the spatial distribution of PWV was correlated with the spatial distribution of rainfall [44]. Yao’s research showed that rainfall lagged behind PWV by 2–6 h [23]. The experiment results of this paper are similar to these research findings.

Based on Figure 10 and Figure 13, the reasons for the high PWV during the whole rainfall in western and southern Henan Province can be summarized. Firstly, the high PWV before rainfall was due to the continuous convergence of water vapor caused by the blocking of the western mountains. Secondly, the reason for the high PWV during rainfall was due to the influence of typhoons and the western Pacific subtropical highs, and to the warm and humid airflow from the southeast direction continuously recharging the water vapor in the atmosphere of Henan Province, resulting in the PWV being at a high level in the process of heavy rainfall. Thirdly, due to the precipitation in the eastern and northern regions, and the fact that the atmospheric water vapor could not be supplemented in time, the PWV in the eastern and northern regions was smaller than that in the western and southern regions after rainfall.

In addition, in this paper, NWP PWV, GNSS PWV and eigenvalue calculations were calculated by our own program. MATLAB software was used as statistical software. For the illustrations, we used the m_map toolbox.

5. Conclusions

A comprehensive and detailed study of NWP PWV accuracy and the relationship between PWV and precipitation is of great significance for improving early-warning ability for extreme weather, and reducing the loss of people’s lives and property. Consequently, in this paper, the accuracy of NWP PWV estimated from MERRA-2 and ERA5 reanalysis data, and the affecting factors of NWP accuracy, were evaluated using PWV from 120 RS stations uniformly distributed throughout the world, and 2 IGS stations in China. Furthermore, the spatial–temporal relationship between NWP PWV and surface precipitation was analyzed qualitatively and quantitatively by an eigenvalue matching method, using the ‘July 20’ heavy rainstorm in Zhengzhou for analysis. The following conclusions were drawn:

(1)

The PWV of both the MERRA-2 data and the ERA5 data had good consistency with RS PWV and GNSS PWV. Compared with MERRA-2 PWV, the accuracy of ERA5 PWV was slightly higher. Latitude, altitude and season were the influencing factors on the NWP PWV estimation accuracy.

(2)

The change trend of ERA5 PWV was consistent with both 24 h cumulative precipitation and surface precipitation during the ‘July 20’ heavy rainstorm in Zhengzhou. The average OMD and OMT between PWV and surface precipitation during the ‘July 20’ rainstorm in Zhengzhou were 56.63% and 3.68 h, respectively, and the maximum optimal matching degree was 80.3%. The spatial–temporal relationship between PWV and surface precipitation was strong.