Rainfall Nowcasting by Blending of Radar Data and Numerical Weather Prediction Rainfall Nowcasting by Blending of Radar Data and Numerical Weather Prediction

In order to improve conventional rainfall nowcasting, radar extrapolation and high- resolution numerical weather prediction (NWP) were blended to get a 6-h quantitative precipitation forecast (QPF) over the Yangtze River Delta region of China. Modifications and calibrations were done to both the extrapolation and NWP in order to get an integrated result from the two, which mainly included the extension for the extrapolation time and region, intensity and position calibration for the NWP, weighted blending of extrapolation and NWP based on scale and time, and a final real-time Z-R relation conversion. Forecast experiments were done, and results show that the blending technique could effectively extend forecast time compared with conventional radar extrapolation, mean- while applying a positive calibration to the NWP. The overall CSI score of 0 – 6 h reflectivity forecast was better than either single forecast.


Introduction
Practice has shown that deterministic short-term extrapolation forecasts have its predictive capability for forecasting precipitation within a few hours.Using modern high-resolution observation, combined with simple extrapolation model, a satisfactory precipitation prediction result with high time-space resolution can be obtained.When forecast time increases to 2-6 h, optimal prediction result can be obtained by combining the results of short-term extrapolations and numerical models together, and the blended product will be superior to any of the two [1][2][3].
As to various kinds of extrapolation methods, the main difference often lies in the acquisition of extrapolation vectors, and there are mainly three of them: The first one is to determine the position of the same echo at different times by using correlation method.The second one is to use echo recognition tracking technology.By decomposing the precipitation echo into different independent echo units, moving vectors of each unit can be obtained by matching units of adjacent times.The third one is to use observation wind field or NWP forecast wind field directly to carry on the echo extrapolation.All the three methods have shown their advantages in operation, such as the models of TITAN [4], TREC [5,6], and optical flow [7,8].At present, both TITAN and TREC are in operation here at Shanghai Central Meteorological Observatory (SCMO), and of which the modified COTREC method is more commonly used by forecasters [9].
As to forecasts of high-resolution regional numerical models, when they are applied directly to short-term prediction, several aspects need to be considered: The first one is the computation time, as NWP model takes time to do data assimilation and parallel integral computation; the result of the first 1-2 h of integral is often not available in practice.When time increases, its forecast will inevitably deviate from the current reality, and its results will need to be corrected.The second one is the resolution, mainly spatial resolution.The common resolution of regional numerical model is about 3-10 km, which is enough to reflect the small and mesoscale weather system, but it is still rough compared to the radar extrapolation of 1 km, particularly in the case of complex urban boundaries; high-impact disaster weather tends to occur within a few dozen square kilometers, and NWP is generally considered not able to effectively distinguish a system with a size less than five times of its resolution.Therefore, it is necessary to apply calibration to conventional NWP products.
Finally, as to the short-term quantitative precipitation forecast, the prevailing view is that to use the combination of conventional short-term forecasting techniques (mainly radar extrapolation) and numerical forecasts is a direct way to increase forecast to more than 2 h.The Nimrod system in the UK, for example, makes short-term strong precipitation forecasts by giving extrapolation and NWP different weights according to forecast time [10].The NIWOT system from NCAR, on the other hand, tries to modify the size of extrapolation using information from NWP [11].In the Swirl System of the Hong Kong Observatory (HKO), the prediction of NWP is revised by means of phase correction and strength correction, and a hyperbolic function is used to determine the weight of extrapolation and NWP, thus effectively improving the prediction effect of 0-6 h [12].Due to the relatively mature swirl system and good openness, HKO has now collaborated with several advanced local bureaus in China such as the Beijing Meteorological Service and the Zhejiang Meteorological Bureau on the task of blending nowcasting [13,14].
In this study, the authors also try to explore a set of short-term forecasting blending methods for the Yangtze River Delta region by using advanced radar network observation in East China and combining high-resolution numerical model forecast.In the following part of this chapter, we first briefly introduce the extrapolation method, and then we will focus on how to set the system configuration and calibrate the final forecast of NWP to obtain a better blending with the extrapolation.Conclusion and discussion are given in the end.

Blending scheme
The aim of this study is to establish a blended quantitative precipitation product for 0-6 h short-term forecast, covering the Yangtze River Delta region with spatial resolution of 3 km and time resolution of 10 min.Using high-resolution radar network observation extrapolation and specially configured numerical model, blended product is then carried out by real-time model evaluation and calibration, so as to improve the timeliness and accuracy of prediction (Figure 1).

Extended 0-6 h radar extrapolation
In this chapter, an improved TREC method, namely, the COTREC [9], is used for the extrapolation.First, a level-II real-time data quality control is done to remove non-meteorological echoes such as super refraction and geophysical echo.Non-meteorological echoes are mostly concentrated in low elevation angle; thus, it can be eliminated by doing comparative analysis of different elevation data, which in this study is 0.5 /1.5 VCP scan.Then, cross-correlation method (TREC) is used to find the TREC wind field.The wind field is first smoothed out by adopting a 9-point smoothing to avoid sudden change.As wind field retrieved by TREC method is only of value in the region where there is observed reflectivity, so the nonconvergence wind field obtained by TREC is usually weakened to a certain extent, which is particularly the case when it comes to isolated echoes or lined echoes.By introducing wind field from numerical forecast as the guide flow field of TREC wind field, the weakening of TREC wind field can be compensated.In this chapter, a 3-hourly forecast wind field of GFS, a global forecast model provided by the National Centers for Environmental Prediction (NCEP) of the USA, is used to supply the TREC wind field.The missing values of TREC wind is compensated by doing vertical averaging and horizontal interpolation of the GFS wind.Thus, a new TREC wind field (nu0 (i, j), nv0 (i, j)) is obtained.The final wind field of COTREC (u (i, j), v (i, j)) is obtained by applying non-convergence algorithm to the TREC field.Finally, by applying the COTREC wind field to the observation radar reflectivity of t2, the extrapolation at t2 + Δt is completed (Figure 2).
Traditional COTREC extrapolation is effective in about 2 h for a single radar observation.In order to extend the extrapolation time to more than 2 h, on the one hand, it is necessary to enlarge the radar data and extend the original single radar observation to the radar network over Yangtze River Delta region surrounding Shanghai.Thanks to the East China regional data sharing line established during the 2010 Shanghai World Expo, observations of 11 Doppler radars around Shanghai are available for the needs of longer time quality control and extrapolation (Table 1).Furthermore, as theoretically the vector field calculated by the TREC method is only effective for the current time, TREC method tends to lose validity due to the evolution of atmosphere, so when forecast time comes to over 2 h, the average wind field of numerical forecast is gradually taking part to replace the TREC wind field so as to indicate the evolution of the wind field.

NWP configuration for blending
For the NWP system, how to configure the numerical model system to make it meet the needs of short-term precipitation forecast as well as the afterward reflectivity blending is a prior task throughout the system configuration.Some of the key issues here include (1) the need to run a regional NWP model with high resolution to resolve weather processes involving small to mesoscale convections, (2) configuration of parametrization schemes in convection scale to keep strength and development of simulated echo system in basically the same range of the radar observation to be assimilated, (3) reduction of spin-up time of NWP to about 1 h to meet the needs for rapid updated short-term prediction, and (4), finally, as traditional short-term precipitation forecast relies heavily on radar observations, it is a common thought that better radar reflectivity assimilation and prediction are the keys to a better NWP for the blending.In the following part of this section, we focus on the NWP system configuration and radar data assimilation to describe their influence on the prediction of radar reflectivity.

The regional NWP system
Research and practice have shown that model spin-up can be effectively reduced by using Rapid Update Cycle [15][16][17][18].By doing RUC is basically to update numerical forecast of the previous time by continuously adding the latest observation data assimilation, thus to reduce spin-up time, and improve the short-term forecast.In this chapter, the Weather Research and Forecasting Model (WRF-ARW, V3.9) and Community Gridpoint Statistical Interpolation (GSI, V3.6) system were chosen as the numerical model and data assimilation system.In order to integrate with the radar extrapolation, the model configuration was focused on the forecast of reflectivity.Current 12 h forecast configuration includes a model grid of 700*700*51, horizontal resolution of 3 km, Thompson's microphysics parameterization, Noah land surface model, RRTM longwave radiation scheme, and Dudhia shortwave radiation scheme.Referring to the short-term system Rap [19] of the National Oceanic and Atmospheric Administration (NOAA), the overall operation scheme was configured to a partially RUC CYCLE run accompanied with a 6 h catch-up CYCLE over the East China (Figure 3).GFS forecast from the National Centers for Environmental Prediction (NCEP) was used to provide initial and boundary conditions twice a day at 02/14 BT for the catch-up cycle.The ECFLOW software package developed by the European Centre for Medium-Range Weather Forecasts (ECMWF) is used for task scheduling and monitoring (https://software.ecmwf.int/wiki/display/ECFLOW/).By using such a cycle-updated configuration with conventional data assimilation, spin-up time was essentially eliminated even in the absence of radar data; well-structured precipitation was usually able to be constructed within 1 h.
Cumulus and boundary layer parameterization are two other important schemes for model precipitation.Cumulus parameterization directly influences the thermal environment and precipitation forecast of model [20,21].While being the place where most heat and water vapor concentrate in the atmosphere, eddies in the boundary layer were the main trigger mechanism for precipitation, especially in the daytime [22,23].The chosen of different boundary layer parameterization schemes would have important implications for forecasting precipitation [24][25][26][27].In the current system, the UW boundary layer scheme [28] and the GF cumulus parameterization scheme [29] are chosen to be the boundary layer and cumulus parameterization scheme.Cumulus convection and boundary layer parameterization schemes can help the establishment of rain bands within the NWP model, while for our experiment, the selection of parametrization schemes should also be considered along with the adaptability of radar assimilation and extrapolation.With the application of radar data and cycled start-up, precipitation should already form in model start for most cases.The maintenance and development of the initial rain system are more important here.Research has found that the UW scheme showed a stronger mixing of the boundary layer turbulence in the daytime [30] and thus has a good inhibition on the afternoon false alarm of local convection, which is usually the case for regional NWPs.On the cumulus parameterization, the GF scheme was chosen for two reasons: one is that it can be applied to high-resolution numerical models [31] and the other is for the optional shallow convection scheme with GF scheme, which can also help to suppress false convection by distinguishing local convection type.

Radar data assimilation for reflectivity
Currently, weather radar is the only way to a continuous, high-resolution, three-dimensional direct observation of convective scale precipitation.Conventional Doppler weather radar observes radial wind and reflectivity, if one wants to combine them into the initialization of NWP; the key point here is to make radar observation compatible with the thermal dynamic fields of model atmosphere, so that they are neither smoothed away as noises nor to induce false development to the model.By introducing radar radial wind and reflectivity data, Xiao et al. [32] used the MM5 3DVAR system to improve initial wind field and moisture distribution of model, thus improving short-term precipitation simulation.Xue et al. [33] assimilated radar reflectivity data within a time window.By using the ARPS assimilation system, they were able to retrieve the distribution of water vapor and precipitation particles in the area of radar observation and thus to improve the initial condition of model and afterward tornado simulation.Some studies found that [34][35][36] the assimilation of radial wind data was generally less important to precipitation than reflectivity data.
The assimilation of radar reflectivity at present was usually done through "diabetic initialization," which is to convert the observation of radar reflectivity to the model variable such as cloud and precipitation particles.At the same time, in-cloud atmosphere profile was adjusted according to latent heat release and water vapor saturation ratio [37].Researches have shown that these kinds of assimilation can effectively reduce NWP spin-up time and improve forecast from the model start to a few hours later.While in the case of radial wind, conventional variation method (such as 3D-/4D-VAR) was used to do the assimilation.It should be pointed out that for the assimilation of single radar radial wind, the improvement to the NWP's threedimensional wind field was often limited, and it was usually done along with the assimilation of reflectivity to get a better improvement.
In recent years, China made a great effort in the construction of various meteorological observation networks.With the construction of the China's next-generation weather radar (CINRAD) network, even trivial precipitation is now hard to miss from observation.Especially to the east part of China, where radar observations were full covered, the integrated radar reflectivity observation can basically be viewed as the real-time precipitation distribution.After data quality control and interpolation, the complex cloud analysis scheme in GSI system was carried out to assimilate the reflectivity data.Model profile and particle mass variables were thus moderated to contain the observed precipitation system.
As shown in Figure 4, the reflectivity of 51 Doppler radars in the middle and east part of China was used.The new-generation weather radar 3D mosaic system [38,39] of the Chinese Academy of Meteorological Sciences (CAMS) was used to integrate the raw data into a three-dimensional mosaic reflectivity.The system was used to apply simple quality control to single radar observation and convert it to three-dimensional Cartesian coordinates.Afterward, data from each radar in the region was integrated together to the three-dimensional reflectivity mosaic.As GSI needed mosaic data of 31 levels which was a bit different from CinradMosaic, another vertical interpolation to the mosaic data was applied.Finally, the mosaic data was converted into Bufr formats onto the horizontal grid of NWP using the GSI radar data interface, so to be used in the cloud analysis [40,41].
Due to the fine coverage of the radar network observation, the distribution of precipitation can be basically indicated at model initial time through assimilation.In order to perform a better blending with the radar extrapolation afterward, the NWP system also needs to retain as much information as possible from the initial radar assimilation.If the difference between observation and prediction was so large that one can hardly relate one to the other, trying to blend them would certainly be difficult.Therefore, before the assimilation of reflectivity, all the initial precipitation particles from the original background were removed, so that the simulated reflectivity of model after assimilation was consistent with the observation (Figure 4).Although eliminating previous particles could reduce false precipitation to a certain extent, due to the application of hourly update cycle, radar observations were continuously added into the NWP.Adjustments done by cloud analysis could still easily lead to overloading of the NWP's hydrometer fields, which further caused the overdevelopment of precipitation.One way to prevent this problem in RAP was to reduce the overstatement of precipitation by setting up a threshold to limit the hydrometer change in every cycle.Similarly in this experiment, assimilation of reflectivity was restrained in two aspects.Firstly, the assimilation reflectivity was performed only in the region where reflectivity was above 15dBZ.Secondly, the retrieved hydrometer mass field was restricted to the limit of 3 g/kg in high-reflectivity area and 1 g/kg in the low-reflectivity area.

Calibration to the NWP forecasted reflectivity
While it often took the short-term NWP forecast a delay of about 2 h to do data collection and time integral, so even by the time forecasters got the NWP forecast, there would already be errors on intensity and position.In order to integrate the NWP forecast with radar extrapolation, it is necessary to calibrate the NWP prediction according to the current observation before it can be integrated with the extrapolation.The main calibrations here include: 1. Calibration of the intensity.A Weibull function was used here to fit the probability distribution of forecast and observation reflectivity separately.The probability density function is: where k is the shape parameter and λ is the scale parameter; they were used to determine the resemblance of forecast and observation.If the two sets of parameters deviated from above a certain threshold, calibration was carried out by doing a proportional scaling to the forecast reflectivity, so that the distribution of k and λ was closer to the observation ones.Recursive scaling was carried out until the final parameters meet the threshold (Figure 6).The same scaling value was then preserved and applied to the forecasts afterward but should be slightly weakened as forecast hour increased.
2. Calibration of the position.Position error of NWP was evaluated based on a target recognition algorithm of SCMO.Firstly, a low-pass filtering was carried out to remove noise and singularity for both forecast and observation.Target recognition was then carried out by identifying different blocks of distinct echoes.After target recognition, the predicted echoes were associated with the observation target by target by comparing the deviation characteristics such as center of gravity, scale, shape, period, intensity, and so on.Finally, when two targets from forecast and observation were identified as the same echo, the displacements of their centers of gravity were used to adjust the position of the original forecast reflectivity (Figure 7).

Weighted blending and Z-R relation determination
Different weight coefficients were applied to the extrapolation and calibrated NWP prediction in the stage of reflectivity blending.Coefficients were designed basically to make the extrapolation dominant in the early time (0-2 h) and make the NWP prediction smoothly take control later on (4-6 h).The weight function here used was based on the same equation used by Chen et al. [13]:  where t0 is the time when the extrapolation weight took the value of 0.5.At present the value of t0 is adjusted within 3-4 according to the scale of the system empirically.The larger the extrapolation echo scale was, the greater the t0 value would be.
Finally, the blended radar reflectivity is converted to precipitation rate through a Z-R relation function: There were five Z-R functions right now in operation at SCMO (Table 2).In order to decide which one to use in the final conversion, a real-time evaluation for each function was applied.
Precipitations from the ground AWS observation along with their observed reflectivity in the latest 3 h were used as statistical samples, mean absolute errors of these samples for each Z-R function was calculated, and the one with the smallest error was selected to use (Figure 8).The final operation run was a 6-h blending QPF forecast updated every 30 min with 3 km/10 min space-time resolution.

Forecast case and evaluation
The forecast to a rainfall event occurred in the late evening of July 5, 2018, is exhibited below for demonstration.This event happened during the Meiyu season in the southern part of Jiangsu and Anhui provinces.The rainfall was produced by a linear mesoscale convective system (MCS) which occurred on the low-level vortex of the Meiyu front.Early convections formed in the afternoon of July 5 in the middle of Anhui province and then moved southeast.
The precipitation lasted to until 18:00 BT in the evening, and new convections started to form to the southwest side of the former convection in the southwest of Anhui.The convections then organized into a quasi-linear convection zone, which lead to the heavy rainfall in the wide of south Anhui.Figures 9, 10 provide composite reflectivity evolution to the process of both the blending forecast and radar mosaic observation.Forecast started at 21:00 BT July 5 when the linear structure of the convection system was basically formed and about to enhance. Figure 9 shows that the first 3 h forecasts of the blending are mainly the results of extrapolation.The echo in south Anhui moves southeast while sustaining the initial intensity.Compared with the observation, the shape of the forecast echo was still satisfying, but the intensity basically keeps the same as initial.The influence of NWP is gradually strengthened after 3 h (Figure 10).New convection develops to the bordering area of Zhejiang, Jiangxi, and Anhui while echoes to the bordering area of Jiangsu and Zhejiang.While the moving direction of original extrapolation echo is mainly southeast, the 6 h blending forecast overall shows a rainfall system moving south and developing, which is consistent with the observation.Therefore, the blending forecast shows a better effect compared to that of extrapolation and NWP prediction at 0-6 h.
There is also precipitation in the east of Zhejiang at the beginning of the forecast, which was mainly generated by convections in the mountainous area of Zhejiang merging together.It gradually disperses while moving east after late afternoon, which is also the case in the blending forecast.
To objectively illustrate the effects of blending forecast, experiments were carried out for four rainfall events in the later part of Meiyu period from July 1-9, 2018, for forecast examination.With a total of 94 forecasts, the critical success index (CSI) for the reflectivity threshold of 20 dBZ is shown in Figure 11.The forecast score of extrapolation overrides that of the NWP in a longer time than common sense.The NWP only takes over at forecast time after 4 h.Although it may be case dependent, the extrapolation of the mosaic reflectivity does have its validity in longer forecast times in some occasion.And due to the differential weight coefficient adjustment, the blending forecast is able to keep resembling to that of extrapolation within 4 h overmatching the NWP.When the extrapolation score falls below the NWP after 4 h, weight of the NWP forecast begins to dominate.The fact that blending forecast scores override the NWP scores in all times shows the positive effect of NWP calibration to the forecasted reflectivity.

Conclusions and discussions
This chapter introduced the 0-6 h blending QPF forecast in SCMO, especially the configuration and calibration of NWP in it.The forecast is based on an extended extrapolation of 11 Doppler radars over the Yangtze River Delta and reflectivity forecast of a high-resolution rapid updated NWP system.Intensity correction based on Weibull transform, position correction based on target recognition and matching, weight coefficient adjustment based on precipitation scale, and Z-R relationship determination based on real-time AWS observation are used to calibrate and blend the NWP forecast.
1. Extrapolation time can be extended when combining observations of multiple radars in the region and using NWP wind forecast to do the extrapolation when forecast time increases.

2.
As NWP is still not perfect for all the forecast elements, configuration of the NWP system is better done according to one's special purpose of use.For short-term efficiency of the reflectivity forecast, attention needs to be paid to decrease model spin-up, improving radar data assimilation and various parameterization schemes.Based on the choice of the UW boundary layer scheme, the GF cumulus scheme with shallow convection, complex cloud analysis of three-dimensional reflectivity, and rapid update cycle along with catch-up configuration, the NWP system is able to sustain initial observation echo and forecast its development in the following hours and, at the same time, keep a good reduction to false alarms of model convections.
3. Calibrations to the NWP forecasts are carried out based on real-time error analysis before they can be integrated with the extrapolation.Spectrum distribution of forecast reflectivity is adjusted by the standards made by a Weibull distribution fitting.Echo positions are adjusted by target identification and matching according to parameters such as weight, scale, and strength.The final QPF is achieved by real-time Z-R relation determination.Forecast experiments were carried out, and results show that blending extrapolation and the NWP prediction can get a better short-term prediction of 0-6 h with the advantages of both.
4. Some of the technical scheme used in this experiment is still immature.For example, for the radar radial wind assimilation, although in some cases it might have got some improvement to the movement and intensity of precipitation, but due to the data coverless and lack of quality control, it was common that there were severe boundary effects which contaminated forecast later that we had to exclude it in the system.For the extrapolation, some state-of-the-art machine learning algorithms can already do forecasts with echo intensity development.And for the areas without radar coverage, how to use satellite data to get a 3D precipitation observation as a substitution is sure to be key to get a better shortterm forecast there.

Figure 1 .
Figure 1.Flowchart of the blending forecast.

Figure 3 .
Figure 3. Model domain and configuration of the rapid update cycle (blue arrows) and the catch-up cycle (green arrows).

Figure 5
provides a comparison of the 3-9 h radar composite reflectivity forecast and observation for a rainfall event in Central China on April 23, 2018.As shown in the figure, the system could maintain radar reflectivity observation from assimilation and continue to forecast its development up to about 9 h.Meanwhile, the overall intensity and position of the simulated reflectivity are kept consistent with the observation.In this way, both the prediction and observation had a good similarity in intensity and shape; then, it is possible to calibrate the forecast echoes based on the comparison of the two.

Figure 4 .
Figure 4. Distribution of radar stations (purple triangles) used in assimilation and retrieved initial composite reflectivity (shaded) from assimilation at 11:00 on June 5, 2018.

Figure 6 .
Figure 6.Intensity adjustment for a 3 hour forecast started at 11:00 BT 20th April, 2016.(a) NWP prediction (original), (b) NWP forecast (after calibration), and (c) Weibull fit for the forecast before/after calibration and the observation.

Figure 7 .
Figure 7. Position adjustment for a 3 hour forecast started at 11:00 BT 20th April, 2016.(a) target identification for observation (after low pass filtering), (b) target identification for NWP prediction (after low-pass filtering), and (c) target identification for NWP prediction (after target adjustment).

Figure 11 .
Figure 11.Critical success index for 94 forecast experiments of four rainfall events in July 2018 (threshold: 20dBZ).

Table 1 .
Radar sites around Shanghai used for extrapolation.

Table 2 .
Coefficients of the five Z-R relation functions in operation at SCMO.