Quantitative Precipitation Forecast Experiment Based on Basic NWP Variables Using Deep Learning

In 1950, numerical weather prediction (NWP) was first performed on an electronic computer at Princeton University (Charney et al., 1950). Since then, NWP has been widely adopted and has long played an increasingly important role in weather forecasting, gradually becoming established as a foundation of modern weather forecasting operations.

The Navier-Stokes equations (including the effects of the Earth's rotation), the first law of thermodynamics, and the ideal gas law represent the full set of prognostic equations with which spatiotemporal changes of wind speed and direction, pressure, density, and temperature in the atmosphere are conventionally described (Kalnay, 2002; Bauer et al., 2015). As a result, the prediction of basic meteorological variables such as air temperature, wind, geopotential height, and humidity are relatively accurate. However, some physical processes, such as convective activities, cloud microphysics, turbulence, and radiation, among others, cannot be resolved by these basic atmospheric equations, so parameterized schemes (PSs) are critical (Bauer et al., 2015).

Different NWP models may use various PSs containing some key characteristics of the weather process simulation and play decisive roles in the performance of precipitation forecasting methods. However, the physical PSs involve some inherent uncertainties, owing to either incompletely understood physical processes or the impacts of impossible or invalid scale resolutions on NWP grids (Bauer et al., 2015). Therefore, such PSs have become a vital bottleneck, restricting further improvements in the performance of NWP forecasting methods.

Many research efforts have demonstrated that deep learning (DL; Table A1 in Appendix) has powerful capabilities in extracting features. DL methods have performed well on tasks involving the automatic extraction of relevant features from a large number of historical observations or NWP reanalysis data and have enabled the construction of models with effective forecasting capabilities (Shi et al., 2015; Zhou et al., 2019; Sønderby et al., 2020). The DL method can extract the characteristics of convective evolution patterns from observational data to provide nowcasting or very short-range forecast results. Shi et al. (2015) established a DL model for radar echo extrapolation based on historical evolution characteristics of radar echoes, which demonstrated better forecasting capability than the traditional optical flow method. Zhou et al. (2020) constructed a DL semantic segmentation (Table A1 in Appendix) model to extract lightning occurrence and development characteristics from multisource observation data, such as radar echoes, satellite cloud images, and lightning density information. The model realized effective lightning forecasts approximately 0–1 h before strikes, with suitable performance forecasting lightning initiations. Sønderby et al. (2020) used satellite, radar, and precipitation data to construct a predictive DL model with a high spatiotemporal resolution for 0–8 h forecasts. The accuracy of precipitation predictions of 0.1 mm to 2 mm in depth exceeded the high-resolution numerical mode-HRRR (high-resolution rapid refresh).

DL has also demonstrated the ability to extract relevant features of synoptic systems from numerical modeling data. Zhou et al. (2019) used convolutional neural networks to extract the characteristics of an environmental field around a set of grid points and thereby realized effective forecasts of thunderstorms, short-term heavy precipitation (>20 mm h^–1), hail, and thunderstorms gusts, outperforming subjective forecasts by meteorologists. Lin et al. (2019) used WRF (weather research and forecasting) model data and the lightning observation data to construct a DL model based on ConvLSTM (convolutional long short-term memory) for extracting features from NWP and observational data. Considering the temporal and spatial variation of multisource data, their DL model produced lightning forecasts for both 6 h and 12 h periods. Because DL effectively integrates observational WRF data, these methods have demonstrated significantly better forecast performance than WRF convective PSs alone.

DL has been applied widely in precipitation nowcasting based on radar data or multisource observations (Shi et al., 2015; Sønderby et al., 2020). Although many machine learning methods have been used to correct NWP forecasts, using DL to extract the rainfall features from NWP to obtain QPF and replace the PS has rarely been explored.

Inspired by the abovementioned prior works, we aim to develop DL methods to go beyond the conventional understanding of physical laws and directly establish non-linear causal relationships between predictors and labels to avoid the impacts of the uncertainty of physical PSs on complex weather forecasting.

Specifically, the forecast performance of basic variables by DL methods has been shown to be significantly better than quantitative precipitation forecasts (QPFs) for most current NWP models. Because basic meteorological variables already include atmospheric movements with ascending and descending air masses, water vapor, and thermodynamic conditions that induce rainfall, it is possible to construct automated systems to establish non-linear mappings between basic meteorological variables and precipitation by DL. Such methods are expected to enable the realization of a better QPF.

Recent attempts have been made to drive NWP using DL directly. Dueben and Bauer (2018) designed a "toy model" to explore the possibilities of performing NWP with machine learning (ML). Arcomano et al. (2020) generated 20-day global forecasts with reservoir computing, a framework for computation derived from recurrent neural network, using ERA5 that were superior to climatology or persistence over three-day prediction periods. Compared with NWP models, ML demonstrates advantages in forecasting physical processes which involve complex parameterization. Rasp et al. (2020) also showed similar results. Although considerable further development is required to perform NWP by ML in operational applications, logistic regression and convolutional neural network models have shown considerable promise on tasks related to predicting variables such as Z500, T850, and T2m. Thus, using ML to improve or drive NWP has become possible in practice.

In this study, a semantic segmentation DL model was constructed to directly establish a mapping relationship between observed precipitation and basic meteorological variables from the reanalysis data (ERA5), and the basic variables including air temperature, wind, geopotential height, and humidity at various atmospheric levels. The basic forecast variables of NWP were fed to the trained model to derive quantitative precipitation prediction. To demonstrate the effectiveness of QPF with DL, we analyzed the importance of several variables and discuss the predictability of precipitation.

The remainder of this study is structured as follows. The data used and the proposed DL model are introduced in section 2. The performance of QPFNet is evaluated in section 3, and its limitations are discussed in section 4, along with the importance of meteorological variables in different pressure layers. Finally, section 5 presents our conclusions.

ERA5 is the fifth generation ECMWF atmospheric reanalysis of the global climate.Reanalysis combines modeling data with observations from across the world into a globally complete and consistent dataset that adheres to the laws of physics (Copernicus Climate Change Service, 2017). ERA5 provides hourly estimates for many atmospheric, ocean-wave, and land-surface quantities. The raw spatial resolution for surface and high-level areas is 0.125° × 0.125° and 0.250° × 0.250°, respectively. In the present work, the spatial resolution is unified to 0.250° × 0.250°, consistent with the QPFNet’s forecasts. In this work, the ERA5 data from June 2016 to June 2019 was used to build the training set.

Some predictors from ERA5 reanalysis data are listed in Table 1. Terrain altitude was included alongside the basic atmospheric variables of standard pressure levels and near the ground. The experimental forecast area spanned 18°–54°N, 72°–135°E (The topographic distribution can be seen in Fig. 1). In performing rainfall prediction, basic deterministic forecast variables were fed into QPFNet from ECMWF's highest-resolution model (HRES) with a spatiotemporal resolution of 3 h and 0.125° × 0.125°. The HRES data from July to September of 2019 was used to build the test set.

Figure 1.  The topography distribution (shaded; in units of km) of the study domain (18°–54°N, 72°–135°E.).

Observational precipitation data was collected from the Multi-source merged Precipitation Analysis System (CMPAS-V2.1) by the China Meteorological Administration (CMA). The data sources used by CMPAS include observational precipitation from more than 40 000 gauges, satellite-derived rainfall from FY2, CMORPH (CPC MORPHing technique), and radar-derived precipitation. The CMPAS-V2.1 data from June 2016 to September 2019 labeled the samples. An independent verification showed that its accuracy was higher than that of any of the three sources of precipitation (Xie and Xiong, 2011; Pan et al., 2015). The spatiotemporal resolution of the data was 0.01° × 0.01° and 1 h.

In consideration of the requirements of operational weather forecasting, the accumulated precipitation of three hours (R_3h) was classified into 102 classes, corresponding to precipitation depths of {0, 0.1, 1, 2, 3……99, ≥100} mm. Within the proposed network architecture, a softmax classifier (Table A1 in Appendix) is used to obtain the probability for each category.

Consistent with the physical laws that govern precipitation, the QPFNet DL forecasting model comprised three-dimensional (3D) convolutional layers, residual connection layers, pooling layers, upsampling layers, and attention layers, among others, as shown in Fig. 2. The predictors of QPFNet are listed in Table 1, while the labels were taken from observed precipitation information from CMPAS-V2.1. Based on the abovementioned research on weather forecasting tasks and DL methods, we explain the role of the proposed neural network architecture in performing QPF.

Figure 2.  Semantic segmentation model (QPFNet) based on an encoder-decoder architecture with a residual mechanism, a convolutional block attention module (CBAM), and skip connections (Copy & Crop, black arrows).

Synoptic systems present a typical 3D structure. The different vertical distributions of variables, such as forward sloping troughs and temperature inversion layers, etc., can be useful in assessing weather system development from a meteorological perspective (Sun and Tao, 2012). Therefore, further development of 3D DL models is considered necessary to extract generation, development, and dissipation features of synoptic systems.

ERA5 reanalysis data provides 3D atmospheric structures with various fundamental elements from high altitudes to the surface. Based on ERA5 reanalysis data, QPFNet was designed to extract 3D features of atmospheric evolution and thus consists of 3D convolutional layers, 3D max-pooling layers, 3D upsampling layers, etc.

Meteorologists must analyze various sources of information to predict the evolution of weather systems and often seek to extract key information from various dynamical, thermodynamical and environmental features to support weather prediction. Similarly, an encoder-decoder architecture was adopted to construct an overall U-shape network architecture. During the encoding process, the encoder continuously extracts useful precipitation information while the feature maps (Table A1 in Appendix) are continuously compressed. Finally, the feature map size is reduced to 1/16 of the input size.

The detailed characteristics of the feature map are partly lost by the compression encoding, so the decoding process is necessary to restore the lost information. During the decoding process, the QPFNet model deduces the forecasting results in each grid. The decoder has the same number of convolutions and the same number of blocks. Instead of pooling, the decoder performs upsampling using upsampling layers that replace the pooling layers in the network architecture. A mirrored decoder in the decoder network connects their corresponding encoder feature map(s) by skip connections (copy and crop, black arrows in Fig. 2).

Residual connections are a type of skip connection that learns residual functions by referencing layer inputs instead of learning unreferenced functions, which were first proposed by He et al. (2016). The depth of a DL network is considered crucial (Simonyan and Zisserman, 2015; Szegedy et al., 2015). However, increasing network depth is associated with many challenges, such as gradient explosion and vanishing. Due to the increase in the number of layers, the gradient may become unstable during the backpropagation process. Residual connections can automatically learn identity mappings to accelerate training. As QPFNet consists of 124 layers, we apply residual connection in the proposed model.

To investigate precipitation processes, meteorologists analyze synoptic-scale systems such as cold vortices, upper-level troughs, typhoons, etc., as well as mesoscale systems such as lower-level convergence lines and cold pools to perform a comprehensive analysis. Hence, a good QPF must comprehensively extract global and local information from weather conditions.

Convolutional Block Attention Modules (CBAM) were proposed by Woo et al. (2018). In the proposed approach, CBAM was used to calculate correlation coefficients between feature maps and a label matrix so that a DL network could better extract global features. As a result, the DL network was able to acquire a better understanding of weather systems of various scales and then forecast occurrences of convective weather.

Average and maximum pooling (Table A1 in Appendix) aggregate the spatial information of the feature map. Then, it is sent to a shared multilayer perceptron network (MLP) to compress the spatial dimensions of the input feature map. The channel’s attention map M_ch was generated by summing pixel by pixel.

where F is the feature map, ${F}_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{ch}}$ and ${F}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{ch}}$ are feature maps with channel information after global average pooling and maximum global pool. W₀ and W₁ are the weight parameters of two layers in MLP, and $ \sigma $ is the sigmoid activation function (Table A1 in Appendix).

The average pooling and maximum pooling processes were also applied to compress the input feature map at the channels. We connect these layers with 7 × 7 convolutions ($ {\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{v}}^{7\times 7} $) to generate a spatial attention map M_s as given below.

$ {F}_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{s}} $ and $ {F}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{s}} $ are the feature maps with spatial information after global average pooling and maximum global pool.

In general, QPF tasks involve continuous value prediction, which falls under the category of regression models. Instead of modeling the continuous values of the targets y, we discretize the variables into many small intervals that continuously cover a range of interest (Allwein et al., 2000; Sønderb et al., 2020).

The probability of all categories is the output of the softmax classifier. Suppose we assign the category with the highest probability to be the predicted category. In that case, the model is inclined to predict light rain because the larger is the rain rate, the smaller its climate probability.

Therefore, a new strategy was designed. Assuming that the probability of each category output by softmax is $ {p}_{c} $ and a probability threshold τ_c is set for each precipitation level c, the value of τ_c decreases as c increases, consistent with the climate probability. For each c, the algorithm calculates the cumulative probability ${p}'_{c}$ of all categories that exceed c. If ${p}'_{c}$ is greater than τ_c, then c is taken as the predicted category. If multiple values of ${p}'_{c}$ are greater than τ_c, the largest c is selected.

where argmax is the operation that finds the argument that gives the maximum value from the target function,

and max is the maximum function, y_i is the predicted category value, m is the total number of categories. τ_c is calculated according to the evaluation results from the validation set. For each c, the cumulative probability ${p}'_{c}$ is assigned the highest threat score (TS; see section 4.2) so long as the BIAS remains less than V_BIAS, and then ${p}'_{c}$ was confirmed to be τ_c. V_BIAS could be set as required, and we set it to be 1.5 in this study.

QPFNet is designed to pay more attention to heavy rainfall areas, which improves the QPF, especially for heavy precipitation events. Nonetheless, QPFNet is also designed to avoid paying excessive attention to the heavy rainfall grid in the training process. Doing so may lead to failure to predict occurrences of weak rainfall events. Moreover, QPFNet is designed to focus on difficult samples to enhance its forecast capability in such cases.

Therefore, the multiclass focal loss [Eq. (4), Lin et al., 2017] is applied as given below.

where we apply $ \mathrm{\alpha }\;\mathrm{a}\mathrm{n}\mathrm{d}\;\mathrm{\gamma } $ as weights whose default configurations are 0.25 and 2.0, respectively. The parameter $ {y}_{i,c} $ is the label for sample i, c is the category, $ {p}_{i,c} $ is the prediction probability of sample i for category c, and m is the total number of categories.

The predictors used by QPFNet are listed in Table 1, labeled by the observational precipitation data from CMPAS2.1. Predictors of time T were labeled by the precipitation of [T–(T + 3)] h. Finally, 6542 samples were constructed, spanning June 2016 to June 2019. A total of 1308 samples (about 20%) were randomly selected to generate the validation set, and the remaining 5234 samples were selected as the training set. The test set contained 1200 samples from July to September of 2019.

The forecasting area was 18°–54°N, 72°–135°E. Consistent with a spatial resolution of 0.25°, the width and length were 145 and 253 grids, respectively. The size of the samples was N_layer × 145 × 253 × N_predictor, where N_layer is the number of layers (from the surface to 100 hPa, a total of 16 layers), and N_predictor is the number of predictors, including air temperature, humidity, pressure or geopotential height, U wind, V wind, and vertical velocity (w). Because there was no win the surface layer, we replaced it with topographic altitude.

The Adam algorithm (Kingma and Ba, 2014) was utilized as an optimizer, and the learning rate was set to 10⁻⁴. Other parameters were held at their default values (Perol et al., 2018). The model was trained for 100 epochs, with a batch size of 2.

An early stopping strategy was used during the training to avoid overfitting. When the loss on the validation set was no longer reduced in 10 epochs, the training process was terminated, and the model weight with the minimum loss on the validation set was saved.

The NVidia CUDA (Compute Unified Device Architecture) library and an NVidia Tesla graphic processing unit (GPU) were used to perform the training and forecasting processes of QPFNet.

Using the GPU, the training time was roughly 27 hours. Moreover, 0–72 h forecasts (at 3-h intervals) at a 0.25° × 0.25° resolution over a study area of (18°–54°N, 72°–135°E) could be completed in 5 min, which is suitable for practical forecasting operations.

After training, the optimal weights of QPFNet were obtained. The basic forecast variables from HRES were fed into QPFNet to obtain the QPF of corresponding forecast times. A comprehensive evaluation and a case evaluation of QPFNet and HRES are given below.

The evaluation results of the QPF performed by QPFNet and HRES are listed in Table 2. Classical evaluation measures, including the probability of detection (POD), false alarm ratio (FAR), threat score (TS), bias, F1-score, and equitable threat score (ETS), were used to evaluate the forecasts generated by the two models.

The overall evaluation scores from July to September 2019, with all leading times from 0 to 72 hours, are listed in Table 2. It was observed that the performance of QPFNet’s QPF was better than that of HRES in the following three aspects.

(1) In predicting drizzle (R_3h ≥ 0.1 mm), QPFNet improved the POD of drizzle forecasts while reducing the FAR, resulting in an obvious improvement in ETS and TS. The TS of QPFNet improved by 19.7% compared to the HRES, and the BIAS was relatively closer to 1.

(2) In predicting heavy rain (R_3h ≥10 mm and R_3h ≥ 20 mm), the PODs of QPFNet obviously increased, while the FARs decreased slightly compared to the HRES, while the POD of QPF of 10 mm of precipitation increased from 0.132 to 0.264, while that of 20 mm depth increased from 0.051 to 0.159. As a result, TS and ETS were obviously improved. The TS of 10 mm (3 h)^–1 and 20 mm (3 h)^–1 precipitation increased by 43.2% and 87.1 %, respectively.

(3) In predicting heavy rain (R_3h≥ 10 mm and R_3h≥20 mm), the QPFs produced by HRES were relatively conservative predictions. HRES’s BIAS was less than 0.76, indicating the forecasting area was sub-optimally small. In contrast, the forecast area of QPFNet was relatively large, with a BIAS larger than 1.5.

Some valuable information may be extracted from Fig. 3. First, for varying precipitation intensities and forecasting lead times, the TS values of QPFNet were better than those of HRES. The greater the precipitation intensity, the larger the improvement in TS. Second, with increasing forecast lead time, the performance of the QPFNet and HRES showed a clear downward trend, consistent with the conventional understanding of weather predictability.

Figure 3.  Threat score (TS) distribution of 12–72 h precipitation forecasts from QPFNet and HRES for 0.1, 3, 10, and 20 mm (3 h)^–1. (a) Evaluation with grid observation from CMPAS2.1; (b) evaluation based on observations from 2411 national meteorological observation stations in China.

Precipitation was recorded over large areas in China from 17–20 July 2019. The precipitation intensity ranged from 0 to more than 50 mm (3 h)^–1, providing a good case to evaluate the performance of the QPFNet. The QPFs of both QPFNet and HRES are shown in Fig. 4. The initial time was 0000 UTC on 17 July 2019.

Figure 4.  Precipitation forecasts of QPFNet and HRES initiated at 0000 UTC on 17 July 2019 for the next 12–72 h (the tables list the corresponding TS of QPFNet and HRES of each forecast).

The following facts may be observed from Fig. 4.

(1) For light rain (0.1 mm ≤ R_3h < 3 mm) events, the QPF of QPFNet was more conservative than that of HRES. The BIAS of QPFNet was closer to 1, consistent with the results in Table 2. In the experiments conducted, QPFNet demonstrated the capability to provide more accurate predictions for rainfall occurrences. The TS of QPFNet forecasts of varying lead times was also better than that of HRES.

(2) For moderate or heavier precipitation (R_3h ≥ 3 mm), the QPFNet's forecast area was less conservative. While HRES missed many heavy precipitation grids, QPFNet hit an obviously larger number of heavy precipitation grids, although QPFNet also gave more false alarms. This is also consistent with the evaluation results in Table 2. Obviously, the forecasts of QPFNet would be more informative in practice.

(3) The precipitation area produced by QPFNet was smoother than that of HRES, indicating that there would be less information in the forecasts produced QPFNet. This may result from the QPFNet seeking the global minimum loss during training, which would lead to some details being smoothed.

In general, QPFNet showed better QPF capability than HRES, both for drizzle and heavy precipitation. However, the distribution of QPFNet’s QPF was smoother, and some detailed information may have been lost, which may confuse meteorologists in operational applications.

QPFNet was trained with ERA5 reanalysis data, and then basic forecast variables from HRES were fed to QPFNet to perform QPF. Therefore, the performance of the ERA5 reanalysis data ultimately determined the performance limitations of QPFNet.

Some heavy precipitation events produced by convective systems are beyond the predictive capability of global NWP models (Hoskins, 2013). Moreover, the performance of QPFNet was also dependent on the accuracy of the basic meteorological variables produced by HRES and the predictability of rainfall events themselves. Therefore, it is necessary to discuss further the performance ceiling of QPF performed by QPFNet.

Corresponding to the various rain rates, the difference in the performance of QPFNet with the reanalysis field (P_a) and the forecast field (P_f) can be found in Fig. 5.

Figure 5.  Performance of QPFNet's QPF by feeding HRES’s basic variables forecast data (circles, 0–72 h, 3-h interval), and the performance of QPFNet simulating precipitation by ERA5 data at corresponding times (squares) (the triangle shows the QPF of HRES; red, green, blue, black, and yellow represent a R_3h of 0.1 mm, 3 mm, 10 mm, 20 mm, and 30 mm, respectively. The evaluation period was from July to September of 2019).

(1) The difference between P_a and P_f was negligible in forecasting light rain (0.1 mm) or the lack thereof. In general, most of the regional drizzle and light rain events were usually produced by synoptic-scale systems and with high predictability (Palmer and Hagedorn, 2006). This indicates that HRES already exhibited a good forecasting ability for synoptic-scale systems.

(2) For moderate rain [3–10 mm (3 h)^–1], the more intense the precipitation, the greater the difference between P_a and P_f. Therefore, we considered that the forecast errors of HRES on basic meteorological variables induced moderate rain events to be missed or incorrectly forecasted.

(3) For heavy precipitation of over 10 mm (3 h)^–1, the convective precipitation is the prevailing mode of rainfall, and synoptic-scale NWP models with sub-grid PSs generally perform poorly (Bauer et al., 2015). Both P_a and P_f were very poor. The negative feedback mechanism between the grid-scale and sub-grid scale in NWP may lead to a sharp decline in heavy rainfall forecast performance.

The permutation importance (PI) method was first proposed by Breiman (2001) for the random forest algorithm. Fortunately, PI can also be used in traditional machine learning and DL models.

By randomly shuffling the predictor x, we obtained a prediction ${y}'$. We then compared the performance difference (ΔP) between ${y}'$ and the normal prediction y. The larger ΔP, the more critical x is indicated to be. For very small values of ΔP, x may not be essential, or there could be information redundancy, such as a linear correlation, between x and other predictors. Because Threat Score (TS) is a crucial criterion in evaluating forecast productions, TS was selected to evaluate the relative importance of each predictor. Ten random shuffling operations were performed for each predictor, and the average result was taken to reduce random errors in the calculation process. The relative importance of each predictor x was defined by $ {\mathrm{T}\mathrm{S}}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e},x}^{} $ as formulated below

Here, $ {\mathrm{T}\mathrm{S}}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{l},x}^{} $ is the TS of predictor x with the original arrangement, and $ {\mathrm{T}\mathrm{S}}_{\mathrm{s}\mathrm{h}\mathrm{u}\mathrm{f}\mathrm{f}\mathrm{l}\mathrm{e}\mathrm{d},x}^{} $ is the TS corresponding to predictors x with a random shuffling.

To further analyze the reasons for QPFNet’s improved QPF performance, the impact on forecast performance of different combinations of various variables was investigated with PI. In addition, the sensitivity of QPFNet’s predictors was also analyzed.

The investigations focused on 0–12 h forecasts (3-h intervals) to minimize the impact of increased forecast errors with increasing lead times. Furthermore, because the topography of the west of China is very irregular, rainfall events, which occurred over the middle and eastern areas (20°–40°N, 110°–125°E), were chosen to estimate the sensitivity of the predictors. Figure 6 shows the impact on precipitation of six basic meteorological variables, taken as a whole, in different pressure layers.

Figure 6.  The relative permutation importance of meteorological variables in different pressure layers (0–12 h forecasts).

(1) The meteorological variables of high layers play an important role in generating precipitation of various intensities—the more intense the precipitation, the more important the high-level atmosphere. The meteorological variables at high levels from 100 hPa to 300 hPa can obviously impact various rainfall intensities. Specifically, the influence from the 200 hPa level is considered the most significant. The high-level jet stream could induce frontogenesis in lower layers. Prior works have indicated that an accelerating high-level jet stream would be conducive to underlying frontogenesis and tend to result in heavy rainfall (Gao and Tao, 1991). The impacts of the meteorological variables above 200 hPa are relatively lower.

(2) Meteorological variables in the middle-level troposphere are crucial for precipitation. Interestingly, the basic variables at 400 hPa were the most critical for various precipitation intensities, while the basic variables at 500 hPa, which are usually more focused on by forecasters, were less important.

(3) Notably, for precipitation events with intensity ≥ 3 mm (3 h)^–1, the predictors in 900 hPa were most important under middle level (500 hPa–ground surface). There was a sharp decrease of relative importance from 900 hPa to 925 hPa, and the basic variables had even lower sensitivity to various precipitation intensities in levels from 925 hPa to the surface, indicating that the atmospheric evolution at the surface or near-surface layers was less important for precipitation processes in QPFNet. This finding seems to conflict with conventional knowledge, especially for convective rainfall. The traditional forecasting wisdom emphasizes the importance of atmospheric conditions at the surface or on near-surface layers during nowcasting or the forecasting of precipitation systems. For example, convergence lines in the boundary layer, dry lines or dewpoint fronts, surface frontal, outflow boundaries, cold pools, etc., are closely related to conditions triggering or maintaining convective activity (Corfidi, 2003; Schumacher and Johnson, 2005; Sun et al., 2015). The reduction in the relative importance of the variables in the boundary layer in QPFNet may be attributed to the poor simulation capability by NWP in the boundary layer (Warner, 2010; Garcia-Carreras et al., 2013). In the eastern plains of China, 900 hPa approximates the upper limit of the boundary layer, coincident with the level that governs the exchange of heat, momentum, vapor, etc., between free atmospheric layers and the boundary layer. Therefore, it is reasonable that the predictors of 900 hPa were most important among the lower layers.

The above analysis shows that the variables in the 400 hPa and 900 hPa layers were more important than the others in the middle and lower troposphere, respectively. Figure 7 compares the relative importance of the variables in 400 hPa (Fig. 7a) and 900 hPa (Fig. 7b).

Figure 7.  The relative permutation importance of different basic meteorological variables for different precipitation intensities at (a) 400 hPa and (b) 900 hPa (T stands for temperature, H for geopotential height, Q for relative humidity, U for zonal wind, V for meridional wind, and w for vertical wind).

(1) The vertical velocity (w) in 400 hPa layer was the most sensitive predictor among all QPF variables except for cases involving extremely heavy rain. When the precipitation intensity was less than 10 mm (3 h)^–1, the relative importance of the 400 hPa vertical velocity increased along with increasing precipitation intensity. In contrast, when the precipitation intensity increased further, its relative sensitivity began to decrease. This phenomenon is related to the relationship between w and precipitation intensity in the ECMWF HRES. The w in the mid-troposphere was directly related to the precipitation intensity of stratiform clouds, consistent with prior observations. For example, based on satellite observations, Liu et al. (2017) found that the w on the 400 hPa layer is closely related to the precipitation rate. However, when the precipitation is driven by violent convective activity, the vertical movement derived from the grid-scale diagnostic and grid PSs was unreliable. As a result, the greater the convective precipitation intensity, the lower the sensitivity of the w from NWP.

(2) The relative sensitivity of air temperature in the 400 hPa layer was similar to w when the precipitation intensity was weaker than 10 mm (3 h)^­–1, showing the influence of the frontal intensity or the cold/warm air activities in the mid-troposphere on the non-convective precipitation. When the rainfall intensity continued to grow, the relative importance of air temperature declined slightly. However, when the precipitation intensity surpassed 20 mm (3 h)^–1, the relative importance of the air temperature at 400 hPa increased again. This may imply a complex influence of convective activity on mid-tropospheric precipitation processes, partly confirmed by the impact of geopotential height on precipitation intensity. When the precipitation intensity was weaker than 10 mm (3 h)^–1, the geopotential height in 400 hPa was of lower sensitivity and remained nearly invariant with the precipitation intensity. However, when the rainfall intensity reached 20 mm (3 h)^–1, its sensitivity obviously increased with precipitation intensity, and its relative importance even surpassed W at 30 mm (3 h)^–1. These findings indicate that convective precipitation intensity was impacted by the combined effects of the geopotential stability and the mid-tropospheric evolution of the synoptic system.

(3) In the lower troposphere (below 900 hPa), the most important variables were air temperature and w. With the increase of precipitation intensity, the relative PI of air temperature showed a unimodal distribution. Meanwhile, w, as well as the horizontal wind (U, V), fluctuated little. Generally, moderate or weaker rainfall [the intensity ≤10 mm (3 h)^–1] is related to frontal system forcing (Zhu et al., 1981), so the relative importance of air temperature increases with rainfall intensities. For precipitation events dominated by convective activities, the sensitivity of the air temperature decreases as the precipitation intensity increases, which may be attributed to the NWP air temperature forecast errors that often occur when resolving convective activities. The synoptic-scale NWP was almost impossible to accurately resolve due to the complex exchange processes between the boundary layer and free atmosphere and the abrupt and potentially violent changes in various meteorological elements in lower layers resulting from convective activity as the cold pool effect, outflow boundaries, etc.

(4) When the precipitation intensity did not exceed 10 mm (3 h)^–1, the rainfall intensity forecast by the QPFNet was not sensitive to vapor variation both at the 400 hPa and 900 hPa layer, and its relative importance was minor. When precipitation intensity exceeded 10 mm (3 h)^–1, its relative importance increased. This may imply the vertical gradient ($ \partial q/\partial z $) or vertical flux ($ \partial wq/\partial z $) of humidity in lower layers was much more sensitive to the intensity of convective precipitation than that of non-convective rainfall (Sun, 2015).

Based on basic meteorological variables from ERA5, as well as terrain data and the observational precipitation data from CMPAS-V2.1, a semantic segmentation DL model named QPFNet has been proposed, incorporating basic meteorological variables including air temperature, geopotential height or sea-level pressure, humidity, and wind in 16 pressure layers. A residual mechanism, an attention mechanism, and a multi-classification method were applied in QPFNet to fit the non-linear relationship between basic meteorological variables and precipitation. The trained model was employed to perform QPF with inputs from the basic variables forecast of the HRES. Forecast evaluation results show that QPFNet outperformed HRES on various precipitation intensities.

(1) The QPFNet achieved a better precipitation forecast than HRES, showing obvious improvement in the various evaluation indices. For the 0–72 h forecast of precipitation accumulated over 3 h, the TS at 0.1 mm, 3 mm, 10 mm, and 20 mm precipitation depths improved by 19.7 %, 15.2 %, 43.2 %, and 87.1 %, respectively.

(2) In forecasting drizzle (0.1 mm), QPFNet was more precise than HRES in terms of higher POD, lower FAR, and having a BIAS closer to 1. For convective precipitation [the intensity≥ 20 mm (3 h)^–1], QPFNet could better predict the central location of heavy precipitation and obviously improve TS.

The above results show that QPFNet can extract precipitation characteristics from three-dimensional basic meteorological variables to realize an effective QPF. Especially for intensive precipitation, the forecast performance of QPFNet obviously improved upon that of the conventional model. Therefore, using DL to extract precipitation features from basic meteorological variables can provide an important reference for QPF, and avoid the uncertainties of PSs.

Based on the proposed DL model, the predictability of precipitation has been discussed. QPFNet’s performance on precipitation simulations on the ERA5 analysis and HRES forecast fields has been compared. Furthermore, as physical laws of cumulus clouds, lightning, heavy fog, and other weather phenomena are not fully understood, the present work can provide insight and a feasible solution for forecasting these complex weather phenomena.

The experimental results show that the forecast errors of the basic variables from HRES have less impact on forecasting drizzle than intensive rainfall. The greater the precipitation intensity, the more obvious the impact is. However, the error effect induced by basic variables was tiny on the convective precipitation process [≥30 mm (3 h)^–1]. This indicates that the basic variables, either from reanalysis or forecast, could not effectively capture the relevant information on violent convection.

Machine Learning (ML) is often criticized by forecasters as being a “black box” because of a perceived inability to understand on what physical basis ML makes its predictions (Mcgovern et al., 2019), when in reality it may reveal further insights. In this study, the permutation importance is used to analyze the sensitivity of the different pressure levels and meteorological variables to precipitation forecasting. This method is beneficial in investigating the complex mechanisms relevant to different rainfall intensities. Some results are consistent with conventional knowledge, such as the variables on 200 hPa playing an important role in various rainfall intensities. However, some results differed from conventional understanding and may be valuable for understanding rainfall processes and forecasting errors. Some examples are given below.

(1) Out of the entire troposphere, the evolution of the basic variables at 400 hPa, instead of 500 hPa, is most important in rainfall forecasting. More specifically, the vertical velocity (w) in 400 hPa level was the most sensitive factor among all variables on forecasting precipitation intensity, except for extremely heavy rain. If the precipitation was led by violent convective activity, the vertical movements derived from the grid-scale diagnostic and grid PSs were unreliable. As a result, the greater the convective precipitation intensity, the lower the relative sensitivities to w from NWP. The relative sensitivity to air temperature at 400 hPa is beneficial to help understand the complex joint impacts of cold/warm air activity, frontal action, and evolution of stratification instability in the mid-troposphere on precipitation led by different dynamical processes.

(2) The 900 hPa level also had an important influence on precipitation. Located near the top of the boundary layer, the 900 hPa level is an exchange layer for heat, momentum, and water vapor between the free atmosphere and the boundary layer. Therefore, the atmospheric conditions in this layer exert an important influence on precipitation. In contrast, the influence of the boundary layer on precipitation is "relatively unimportant", which may be caused by the fact that the simulation capability of NWP in the boundary layer was much lower than that of the free atmosphere.

The model interpretation and visualization (MIV) of DL provides a new perspective to understand precipitation mechanisms. Specifically, some of the analytical results, such as the critical role of the 400 hPa layer on precipitation, can guide us to better analyze mechanisms related to precipitation. Other analytical results, such as the "relatively unimportant" effects of boundary layer conditions on precipitation, can guide us to better understand the limitations of NWP as well.

Future research should be conducted to improve QPFNet to allow it to become operationally feasible. For example, QPFNet's forecast is smoother than that of ECMWF, indicating a loss of detail. This may be caused by our training strategy and loss function, which sought the global minimum loss value. Moreover, the QPFNet’s forecast map is not a physical field like ECMWF forecast, and it is not drawn from the multivariate distribution of rain maps, leading to difficulties in obtaining realistic maps.

Therefore, a customized loss function and generative adversarial network may be used to enrich the details of the precipitation forecast. In addition, state-of-the-art operational approaches such as non-homogeneous regression or quantile random forests, etc., should be also be evaluated and compared under operational conditions.

Acknowledgements. The authors would like to acknowledge the financial support of the National Key Research and Development Program (Grant No. 2017YFC1502000) and the National Natural Science Foundation of China (Key Program, 91937301), and thank the Copernicus Climate Change Service (C3S) for providing ERA5 reanalysis data, and National Meteorological Information Center for providing the CMPAS-V2.1 product.

APPENDIX