- Email: [email protected]
Nowcasting of 1-h precipitation using radar and NWP data
Journal of Hydrology (2006) 328, 200– 211
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jhydrol
Nowcasting of 1-h precipitation using radar and NWP data Zbynek Sokol The Academy of Sciences of the Czech Republic, Inst. of Atmospheric Physics, Bocni II, 1401, 141 31 Prague, Czech Republic Received 23 March 2005; received in revised form 4 October 2005; accepted 17 December 2005
KEYWORDS
Summary A multiple linear regression model complemented by a correction procedure (REG) is applied to nowcasting of 1-h precipitation for warm seasons. The model aims at predicting the mean area precipitation in the 9 km · 9 km squares for the territory of the Czech Republic. Precipitation amounts are calculated from the radar reflectivity measured by two C band radars, which are operated by the Czech Hydrometeorological Institute. The model predictors are obtained from: (i) radar-derived precipitation in the squares; (ii) radar-derived precipitation advected by wind fields at the 700 hPa level produced by ALADIN/LACE NWP model forecast; (iii) variables derived from forecasts of the ALADIN/LACE NWP model. Accuracy of the REG forecast is evaluated and compared with the forecast obtained by the advection of precipitation fields by the NWP wind at the 700 hPa level (ADV). The ADV forecast serves as one of the predictors (ii). Results show that REG yields apparently better forecasts than ADV in terms of RMSE, correlation coefficients and categorical measures. Also local precipitation maxima have better positions with respect to the real ones. However, if ADV did not indicate any precipitation, REG was not capable to forecast a convective precipitation development. This fact deteriorates the REG accuracy during afternoon hours. We may say that the REG forecast is more accurate in night and morning hours since the situation mentioned occurs quite scarcely during those time periods. c 2006 Elsevier B.V. All rights reserved.
Precipitation forecast; Regression models; Nowcasting; Radar
Introduction
For some hydrological applications, especially the flash flood warning and urban drainage management, make frequently demands of the 0–6 h forecasting (nowcasting) of heavy precipitation. At present, many nowcasting methods employ extrapolations of radar echo. Methods based on a E-mail address: [email protected].
simple extrapolation technique only (e.g. Dixon and Wiener, 1993; Johnson et al., 1993; Mecklenburg et al., 2000) showed a rapidly decreasing accuracy with the increasing lead time. As they do not consider a storm initiation, growth and dissipation the simple extrapolation is usually applicable up to 30–60 min. Radar data including an extrapolated radar echo together with other data provided by satellites, lightning and numerical weather prediction (NWP) models were used
0022-1694/$ - see front matter c 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2005.12.023
Nowcasting of 1-h precipitation using radar and NWP data as predictors in a statistical model by Kitzmiller (1996). Probabilities of exceedence a defined threshold were forecasted by multivariable linear regression models for individual areas of 40 km · 40 km. The regression models were developed on the basis of real data by the application of a technique very similar to the model output statistics (MOS; Glahn and Lowry, 1972). Models for different time intervals during a day were developed separately in order to respect the diurnal cycle of precipitation processes. Another approach to the nowcasting of heavy precipitation is represented by the model NIMROD (Golding, 1998), which utilises both radar echo extrapolation and precipitation forecasts of a NWP model. For shorter time periods the major weight is given to the extrapolation of an existing precipitating field, which is derived from satellite and radar data. With increasing forecast time the weight shifts to the numerical model and becomes almost totally dependent on the NWP model for the 6-h forecast. A different approach is presented by an automated nowcasting system of convective precipitation GANDOLF (Pierce et al., 2000) that attempts to simulate a storm development. It contains a precipitation model, which incorporates a conceptual model of the life cycle of shower clouds. With help of satellite and radar data and a variety of forecast products from a NWP model, convective cells are identified and their possible development is predicted. In this paper, a nowcasting technique based on the application of a linear regression model, similar to that used by Kitzmiller (1996), is applied to the forecast of mean 1-h precipitation in the 9 km · 9 km squares. The forecast is focused on large precipitation amounts, which are important for hydrological applications. Main differences between the applied method and Kitzmiller’s approach are the forecast variables and a significantly higher horizontal resolution. The paper is divided into sections as follows. The next section describes data and a method calculating precipitation amounts from the radar data. Also both forecast and data regions and predictands are defined there. The ensuing section show types of the predictors and an advection algorithm used to derive them. A regression model and a posterior correction of model outputs are given next. The method of model verification and results and discussion are also given. The conclusions are outlined in the final section.
Data and forecast region The data used in this study are from the 2002 warm season and include radar-derived precipitation and prognostic fields of the NWP model ALADIN/LACE (Aire Limite ´e Adaptation Dynamique De ´veloppement InterNational/Limited Area Modelling in Central Europe), which is routinely applied by the Czech Hydrometeorological Institute (CHMI). Fig. 1 shows the data and forecast regions. The radar measurements collected from the data region are used. This region is larger than the forecast one because the radar data are advected by the NWP model wind. Both regions are well covered by radar and NWP model data. Precipitation amounts are derived from radar reflectivity measurements. Radar-derived precipitation cannot be adjusted by ground measurements as there are no available hourly measurements from rain gauges reliable enough. Ra-
201 dar data are recorded by the Czech radar network consisting of two weather radars Brdy and Skalky. Both radars measure in C band and the interval of their scans is 10 min (Fig. 1). Precipitation amounts are calculated by integration of rain rates in time. The rain rates are derived from the radar reflectivity at the elevation of 2 km above the sea level (CAPPI 2 km) with help of a standard algorithm operationally used in CHMI (Zacharov et al., 2004; Michelson et al., 2004). The bright band is not considered here, as it occurs very rarely at or below the elevation of 2 km in warm seasons. Rain rates are calculated by using the Z–R relationship Z = 200R1.6 in the range of 7 dBZ 6Z 6 55 dBZ. For Z < 7 dBZ and Z > 55 dBZ the rates are R = 0 mm/h and R = 99.85 mm/ h, respectively. For both radars the maximum range of rainfall estimation is 256 km and the azimuth resolution is 1. Data from Skalky radar are recorded in 20 elevation angles and the ray bin resolution is 1 km. Radar Brdy uses 14 elevation angles and ray bin resolution of 0.5 km. The raw radar coordinates are transformed into the gnomonic projection and rain rates are interpolated/extrapolated separately for both radars into 512 · 512 pixels with the horizontal resolution of 1 km · 1 km. The rain rates in pixels outside the 256 km range are set to zero. The pixel values are smoothed by using the median from 3 · 3 nearby pixels, which improves radar precipitation when compared with gauge measurements (Kra cmar et al., 1998). Then 20-min precipitation ´ˇ amounts are calculated in the pixels and the maximum of two radar fields is applied in the overlapping pixels to obtain precipitation estimates in the data region. The integration starts at 00, 20 and 40 min for each hour. The data region is covered by squares consisting of 9 · 9 radar pixels. For these squares the mean pixel values of 1-h precipitation amounts are calculated as the sums of corresponding 20-min amounts. The square size selected allows the forecast of precipitation even for quite small river basins, the area of which is about 100 km2. The 9 km · 9 km squares in the forecast region are a subset of the squares covering the data region. The forecast variables (predictands) are 1-h precipitation amounts in squares of the forecast region, consisting of 43 · 25 squares and covering a large part of the territory of the Czech Republic. For simplicity reasons, we will be using the term ‘‘precipitation amount’’ instead of mean area precipitation amount. Outputs from the NWP model ALADIN/LACE include predicted temperature, relative humidity, geopotential, horizontal wind and vertical velocity at levels of 925, 850, 700, 500 and 300 hPa, and temperature, relative humidity, wind at the surface, and the sea level pressure. ALADIN/ LACE NWP model is a hydrostatic limited area model with the horizontal resolution of 12 km. The model integration starts every 12 h and basic prognostic fields are recorded every 6 h. In order to obtain the prognostic values in other times the linear interpolation between two consecutive 6hourly terms is used.
Predictors All predictors are available no later than at time T and the forecast 1-h precipitation is related to the interval [T, T + 60 (min)]. A multiple linear regression employs four
202
Z. Sokol
Figure 1 Data and forecast regions. The coordinates are in radar pixels, the size of which is 1 km · 1 km. The squares show the radar locations, and circle parts show the maximum radar range 256 km.
types of predictors. The first type consists of 1-h precipitation amounts and instantaneous rain rates. Precipitation predictors are determined in the fixed squares as well as in the squares obtained by the advection of precipitation fields. Besides the squares of 9 · 9 pixels, which correspond to the predictand squares, also squares of 17 · 17 and 25 · 25, which have the same centres as the 9 · 9 squares, are considered. The advection makes use of a prognostic wind field at the 700 hPa level and a backward scheme with the time step of 20 min. We considered also a wind field at 500 hPa level and a mean wind field in the layer of 700– 500 hPa. We forecasted 1-h precipitation amounts by advection of the actual precipitation field by various wind fields and compared results. We found out that the forecast by the wind at the 700 hPa level was slightly more accurate than forecasts employing other wind fields. The advection starts at centres of the 9 · 9 squares in the forecast domain (Fig. 1). After each advection time step, 20min precipitation is determined in the squares centred in the pixels resulting from the advection steps. In order to determine accumulated precipitation over interval [t 20, t] we apply the prognostic wind corresponding to the t 10 time. The predictor value of 1-h precipitation is the sum of 20-min precipitation in corresponding squares. The application of three 20-min steps turned out to be slightly more accurate than one 60-min step. The following predictors are determined for times t = T, T 20, T 40 and T 60:
• Current mean and maximum 1-h precipitation in the squares of 9 · 9, 17 · 17 and 25 · 25. • Current mean and maximum instantaneous precipitation rates in the squares of 9 · 9, 17 · 17 and 25 · 25. • Forecasted mean and maximum 1-h precipitation in the squares of 9 · 9, 17 · 17 and 25 · 25 that are produced by the advection algorithm. Mean precipitation in the squares of 9 · 9 for t = T represents the forecast by advection only. • Forecasted mean and maximum instantaneous precipitation rates in the squares of 9 · 9, 17 · 17 and 25 · 25 obtained by the application of the advection algorithm. • Mean vertically integrated liquid water content (VIL) in the squares of 9 · 9. • Mean top of radar echoes (ECHOTOP) in the squares of 9 · 9. • Mean VIL and ECHOTOP in the squares of 9 · 9 obtained by the advection algorithm. The predictors of the second type are prognostic NWP model outputs, which are interpolated in space and time to get predictor values for the individual squares and forecast times. In addition, the following derived parameters are also used: K-index, vertical gradient of equivalent potential temperature in the layers of 700–500 and 850– 700 hPa, the mean vertical relative humidity and vertical velocity (in the layer from 925 to 500 hPa). The list of
Nowcasting of 1-h precipitation using radar and NWP data predictors is complemented by the tendencies of the NWP model predictors over 6-h intervals containing time T. The third type of predictors contains linearized predictors calculated in a way similar to that used by Charba (1998). The idea is to linearize relationships between the predictors and predictands. Predictor values are sorted in ascending order and divided into 10 groups with an approximately equal number of members. In each group mean predictor values, (Xi, i = 1, . . . , 10) and mean predictand values, (Yi, i = 1, . . . , 10) are calculated and Yi are treated as a piecewise linear function of Xi. Linearized predictor values Y(x) are calculated as follows: YðxÞ ¼ Y i
x X iþ1 X iþ1 x þ Y iþ1 ; X iþ1 X i X iþ1 X i
for X i 6 x 6 X iþ1 ; ð1Þ
YðxÞ ¼ X 1 ;
for x 6 X 1 ;
ð2Þ
203 YðxÞ ¼ X 10 ;
for X 10 6 x;
ð3Þ
where x is the observed predictor value. The grouping mentioned above follows from results of a number of preliminary tests. The fourth type of predictors has the form of (a Æ b)0.5, where a and b are predictors of the first type.
Regression models We use a linear model as many others have done in spite of the fact that the relationships between predictors and predictands are supposed to be non-linear. Those relationships seem to be complex and predictors probably do not contain sufficient information about precipitation processes to allow an accurate forecast. Moreover, values of predictors are supposed to be influenced by various types of errors. Therefore, we use a simple linear model, robust and easily applicable, as an approximation of the relationships.
Figure 2 Scatter plots (c) and (d) illustrate effects of the correction procedure for M = 5 and M = 10 values, respectively, on the standard regression model output (a) M = 0. The plot (b) shows the areas from which the pi and yi values are used in calculating corrected output pcor,i. For M = 10 the data from the darkest area are used. For M = 5 are used data from the union of the darkest and light grey areas. The data concerns a calibration subset for 22–24 UTC.
204
Z. Sokol
Table 1 Mean predictand values (Mean) and probabilities (%) that predictand values will reach or exceed thresholds P = 1, 3, 5 and 10 mm for 8 considered daytimes (UTC) UTC
Mean (mm)
P = 1 mm
P = 3 mm
P = 5 mm
P = 10 mm
01–03 04–06 07–09 10–12 13–15 16–18 19–21 22–24
0.85 0.93 0.72 0.68 0.93 0.94 1.16 0.81
24.2 29.0 22.9 18.6 22.5 22.6 27.5 24.7
8.5 9.7 6.4 6.1 8.4 8.8 11.8 7.2
3.4 3.2 2.1 2.6 3.9 4.5 6.0 2.8
0.6 0.2 0.2 0.5 1.3 1.4 1.7 0.4
The predictand is 1-h precipitation amount in the forecast squares 9 km · 9 km.
Following the scheme suggested by Kitzmiller (1996) we have developed separately 8 multiple linear regression models for the forecast times 01–03 UTC, 04–06
UTC, . . . , and 21–24 UTC. The models take a linear form such as y ¼ a0 þ a1 x 1 þ a2 x 2 þ þ aN x N ;
ð4Þ
where a1, . . . , aN are model parameters, y is the predictand value and xj, j = 1, . . . , N, are predictor values. The parameters are found by minimizing the sum S of squared differences between observed yi and predicted values pi: S¼
nc X ðy i pi Þ2 ;
ð5Þ
i¼1
where pi values are obtained with help of formula (4), and they are pi ¼ a0 þ a1 z1;i þ a2 z2;i þ þ aN zN;i ;
ð6Þ
z1,i, . . . , zN,i are predictor values corresponding to the ith predictand value and nc is the number of predictands used to derive the model parameters. The solution of the minimization problem can be found in standard regression texts e.g. Wilks (1995).
Figure 3 Scatter plots (a)–(c) show REG forecasts in dependence on M. The forecast is calculated for 16–18 UTC. In (d) the solid lines show average forecast values for the defined observation categories and dashed lines show average observation values for the defined forecast categories in dependence on M. See the text for details.
Nowcasting of 1-h precipitation using radar and NWP data The predictors are selected by the forward-selection screening algorithm (Wilks, 1995). The selection stops when the decrease of mean-square-error (MSE) by the inclusion of a new predictor is smaller than 0.2% of the previous MSE. The maximum number of predictors was restricted to 16 as the tests showed that additional predictors did not improve the forecast. The choice of the 0.2% threshold was subjectively made with respect to the preliminary tests. The model outputs are smoother than the observed precipitation and significantly underestimate heavy precipitation due to its scarcity within the dataset (for the confirmation see Fig. 2(a). The underestimation of heavy precipitation can be corrected by a posterior correction of model outputs pi by the linear transformation: pcor;i ¼ api þ b;
i ¼ 1; . . . ; nc ;
ð7Þ
where a and b are parameters. Their meaning will be elucidated in the following text.
205 Parameter a is found to minimize the function X FðaÞ ¼ ðy i api Þ2 .
ð8Þ
i;y i ;pi >M
The sum in (8) represents all data i from the calibration data set, whenever yi (observations) or pi (model outputs) are greater than a given value of M. An example of the areas from which the data are taken for M = 5 and M = 10 are shown in Fig. 2(b). In order to determine a for M = 10, only predictor values and corresponding observations from the darkest area are used. In case of M = 5 the light grey area is added to the darkest one. These areas contain mainly those forecasts, in which the model significantly underestimates observations. Therefore, the parameter a is greater than 1 and consequently, api > pi. Parameter b aims at maintaining an approximate equality between the sum of observations and corrected model outputs pcor,i, which yields, when (5) and (7) are employed:
Table 2 RMSE and CC for REG with M = 0, M = 5 and M = 10 applied to 16–18 UTC in dependence on observed precipitation amounts (PP) PP P 0 mm
M=0 M=5 M = 10
PP P 1 mm
PP P 5 mm
PP P 10 mm
BIAS
RMSE
CC
RMSE
CC
RMSE
CC
RMSE
CC
0.96 0.97 0.99
1.70 1.71 1.82
0.68 0.68 0.68
3.33 3.32 3.51
0.54 0.55 0.55
6.83 6.57 6.45
0.34 0.34 0.34
10.55 9.96 9.41
0.24 0.24 0.24
BIAS is calculated only for all data.
Table 3 UTC
RMSE and CC for ADV and REG methods in dependence on daytimes (UTC) and observed precipitation amounts (PP) Method
PP P 0 mm
PP P 1 mm
PP P 5 mm
BIAS
RMSE
CC
RMSE
CC
RMSE
PP P 10 mm CC
RMSE
CC
01–03
ADV REG
0.94 1.01
1.49 1.12
0.61 0.77
2.75 2.06
0.40 0.61
5.64 4.29
0.14 0.35
10.03 7.70
0.07 0.34
04–06
ADV REG
0.99 1.02
1.28 0.91
0.70 0.83
2.12 1.53
0.51 0.66
3.84 2.91
0.23 0.30
6.66 5.87
0.24 0.14
07–09
ADV REG
1.11 0.99
1.17 0.85
0.65 0.80
2.05 1.55
0.43 0.63
4.37 3.41
0.18 0.25
8.30 7.07
0.16 0.24
10–12
ADV REG
0.93 0.94
1.44 1.10
0.57 0.74
3.08 2.35
0.35 0.60
6.63 5.14
0.08 0.33
11.35 8.97
0.06 0.16
13–15
ADV REG
0.96 0.99
2.58 1.76
0.41 0.70
5.02 3.44
0.24 0.60
10.26 7.45
0.05 0.40
15.72 11.72
0.03 0.19
16–18
ADV REG
0.98 0.97
2.69 1.71
0.37 0.68
5.09 3.32
0.18 0.55
9.26 6.54
0.02 0.34
13.70 9.96
0.04 0.24
19–21
ADV REG
0.98 0.99
2.53 1.64
0.50 0.76
4.40 2.89
0.32 0.64
7.73 5.27
0.06 0.43
11.61 8.15
0.10 0.29
22–24
ADV REG
1.12 1.01
1.47 1.03
0.62 0.77
2.56 1.86
0.46 0.62
5.07 3.98
0.25 0.31
7.80 6.89
0.06 0.10
BIAS is calculated only for all data.
206 bð0Þ ¼
Z. Sokol nc 1 X ðy api Þ. nc i¼1 i
ð9Þ
The calculated values of a and b(0) may lead to negative values of pcor,i. Thus, the following iterative process (k = 1, . . . , kmax) is used in calculating b: ðkÞ ð10Þ pcor;i ¼ max api þ bðkÞ ; 0 ; i ¼ 1; . . . ; nc ; bðkþ1Þ ¼ bðkÞ þ
nc 1 X ðkÞ y i pcor;i ; nc i¼1
ð11Þ
It converges in several steps. The test showed that kmax = 5 is sufficient and b ¼ bðkmax Þ . Even after iterations pcor,i can be negative. In this case pcor,i is set to zero. Regression model (4) with the correction algorithm ((7)– (11)) is referred to as REG. The correction (7) is used for all regression outputs pi . The influence of the M value on the model outputs is illustrated by Fig. 2. It shows that the correction improves the forecast of heavy precipitation. On the other hand, it may overestimate low precipitation amounts.
The M values used in this study were derived at from the test results.
Calibration and verification of the models For each of 8 daytimes, 39 terms with large local precipitation are used for model calibration and verification. It means that 41,925 (39 · 43 · 25) predictand values are available for each model. Basic statistical characteristics of the predictand values are shown in Table 1. The majority of large local precipitation, frequently occurring in afternoons in the warm season, is of convective origin, and the forecast of the convective precipitation is very difficult. For each daytime the data are divided into 3 subsets, each containing 13 terms. The data are placed in the subsets by time sequence to maintain effects of synoptic regimes. Two subsets are used as the calibration data and the third one serves for the verification. A cross-validation technique (Elstner and Schmertmann, 1994), choosing different combinations of the subsets, is applied to get the independent
Figure 4 Scatter plots of forecast and observed values for 13–15 UTC. (a)–(c) show REG forecasts: (a) shows all the values, (b) shows only pairs when ADV forecasts positive precipitation, and (c) depicts events when ADV is equal to 0. (d) shows ADV forecast for all data.
Nowcasting of 1-h precipitation using radar and NWP data forecasts for all data. The independent forecasts are used in the following verification tests. The arrangement of the data into the subsets puts close terms, which may reflect similar development of precipitation, in the same subset in most cases. Therefore, they do not have a substantial influence on verification results. The accuracy of the forecast is evaluated by the rootmean-square-error (RMSE), bias (BIAS) and correlation coefficient (CC): sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 1X RMSE ¼ ðy pi Þ2 ; n i¼1 i
ð12Þ
Pn pi ; BIAS ¼ Pi¼1 n i¼1 y i
ð13Þ
Pn
P P y i n1 nj¼1 y j pi n1 nj¼1 pj CC ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 P 2ffi ; Pn Pn Pn n 1 1 y y p p i i i¼1 j¼1 j i¼1 j¼1 j n n i¼1
Figure 5
ð14Þ
207 where pi, yi are forecast and observed values of predictands. As the forecast is primarily aimed at large precipitation amounts, RMSE and CC are calculated with respect to the values of actual precipitation yi. The model results are compared with the forecast obtained by the advection of the precipitation field by the NWP wind at 700 hPa (ADV), which is one of the predictors.
Results and discussion The selected predictors and a number of them differ for daytimes and among the particular calibration data sets. Large variability of the selected predictors results from the fact that the majority of the predictors are highly correlated and therefore, there are various groups of selected predictors that produce comparable results. The screening algorithm yields numbers of predictors in the range from 12 to 16. Even higher number of predictors does not cause an overfitting but also does not improve the model accuracy.
The same as Fig. 3 for 04–06 UTC.
208
Figure 6 Relations between the forecast and observed values by REG and ADV methods for 13–15 UTC (a) and for 04–06 UTC (b). The solid lines show average forecast values for defined observation categories and the dashed lines show average observation values for defined forecast categories for M = 5. See the text for details.
The most important predictors belong to the first and third types. Predictors derived from the NWP model are contained almost in each selection but they have little impact on the forecast accuracy. The same is true for the fourth type of predictors. Although the majority of selected predictors are valid at time T, when the forecast is issued, predictors from earlier times are also frequently included. There is a difference between the most frequently selected predictors for afternoon and other daytimes. For the afternoon (13–15 UTC, 16–18 UTC) the most frequent predictors are: the mean precipitation rates in 9 · 9 pixels at time T and the maximum precipitation rates in 25 · 25 pixels at time T. Other predictors are related to precipitation or precipitation rates and squares 9 · 9 pixels. Some linearized predictors are also selected but not so often, e.g. advected VIL. For the remaining daytimes two predictors: the advected mean precipitation rates in 25 · 25 pixels
Z. Sokol at time T and the linearized advected ECHOTOP at time T are more frequently chosen than others. Other often selected predictors are the maximum precipitation rate in 9 · 9 pixels and advected VIL. The influence of the value of M (8) is demonstrated by Fig. 3 and Table 2 for the interval 16–18 UTC. Three scatter plots of observed and forecast values for M = 0 (standard regression), M = 5 and M = 10 show that REG is able to forecast higher values for larger M. Table 2 shows that with increasing M the accuracy of the forecast increases for higher measured precipitation. On the contrary, the large value of M is reducing accuracy for small precipitation amounts. It is worth mentioning that BIAS is small for all values of M. It proves that the iterative algorithm ((9)–(11)) is helpful in reducing BIAS values. Fig. 3(d) compares the forecasts for different M in a more transparent way. The horizontal coordinates of marks on full lines are calculated as the mean of observed values from intervals [0, 0.5], [0.5, 3.5], [3.5, 6.5], . . . , [21.5, 24.5] and P24.5 mm. The vertical coordinates are the means of forecast values corresponding to observed data from the relevant intervals. The dashed lines contain the marks calculated in a similar way, but at first, forecast values are averaged over the intervals. The standard regression should yield a dashed line close to the line representing the perfect fit. As it follows from Table 2 and Fig. 3 the REG model with M = 5 apparently improves the forecast of large precipitation in comparison with M = 0 and only negligibly decreases the accuracy for small precipitation amounts. REG with M = 10 even more improves the forecast of large precipitation; however, it deteriorates forecasts of small precipitation. Consequently, REG with M = 5 is used in the following calculations. The quality of the forecast depends essentially on the daytime. The forecast values are apparently less accurate (Table 3) and more scattered around the perfect fit for the terms 13–15 and 16–18 UTC than for the remaining time intervals. It is illustrated by Figs. 4–6 for 04–06 UTC and 13–15 UTC. The comparison of Figs. 4(a) and 5(a) with Figs. 4(d) and 5(d) show that REG significantly improves the forecast of ADV. For afternoon the simple advection ADV is definitely unsuitable for the precipitation forecasting. From Fig. 4(c), which shows the forecast for those events in which the forecast by ADV is equal to 0 mm, follows that REG is not capable to produce new convective cells. When these cases are excluded from the data then the REG forecast gives better results (Fig. 4(b)). The mentioned drawback of REG is not very important for morning and night terms (Fig. 5) as a development of new storm cells is rather scarce. Fig. 6 compares REG and ADV forecasts by using the means similar to those in Fig. 3(d). Also Fig. 7, where the quality of the forecast for 13–15 UTC and 04–06 UTC is compared in terms of categorical measures POD (probability of detection), FAR (false alarm rate) and ETS (equitable thread score; Wilks, 1995), confirms a lower accuracy of afternoon forecasts except for the 99% percentile. The continuous forecast is transformed into a binary forecast using thresholds corresponding to the 50%, 70%, 80%, 90%, 95% and 99% percentiles of the observed precipitation. Another drawback of the forecasts common for all terms is the tendency to smooth maximum and minimum values. This is due to minimization of mean-square-error in the regression. Histograms in Table 4 show that REG forecasts
Nowcasting of 1-h precipitation using radar and NWP data both low (<1 mm) and high values with lower frequencies than they are observed. This behaviour of REG can be noted in Fig. 8, where examples of three consecutive forecasts of extremely heavy precipitation are shown. The REG fields are smoother than the observed or advected precipitation fields; however, REG describes shapes of areas with precipitation better than ADV. Heavy precipitation forecasts provided by ADV have significant location errors, which lead to a larger RMSE and lower CC, as shown in Table 5. We tried to find out how well REG can localize large precipitation amounts. Therefore, we studied differences between the positions of observed local maxima and positions of local maxima in forecasts produced by REG and ADV. We consider all inner squares within the forecast domain with observed precipitation at least 10 mm under
209 the condition that 8 neighbouring values are lower (local maximum). For these squares we look for the nearest local maxima in the forecast fields and evaluate distances (in terms of square numbers) between the observed and forecast local maxima. Table 6 shows that the accuracy of the forecast by REG and ADV differs for individual daytimes. For REG about 12% of local maxima are localized well, 55% are localized correctly or in neighbouring squares and in 86% of all cases distances between the observed and forecast local maxima do not exceed the distance of 2 squares (18 km). For ADV the localization is less accurate: 8%, 48% and 78%. The good results of REG for afternoon terms 13– 15 UTC and 16–18 UTC are worth noting. Smooth stratiform precipitation fields, where local maxima are not well pronounced, might influence the evaluation for night terms.
Figure 7 Evaluation of REG forecasts for 13–15 UTC and 04–06 UTC by POD, FAR and ETS. The horizontal axis shows the 50%, 70%, 80%, 90%, 95% and 99% percentiles of observed precipitation. The corresponding values of precipitation amounts are 0.04, 0.58, 1.19, 2.62, 4.29, 11.84 mm and 0.12, 0.94, 1.70, 2.95, 4.19, 7.31 mm for 13–16 UTC and 4–6 UTC, respectively.
Table 4 Frequency in % of observed (OBS) and forecast (REG) values for the given precipitation intervals in mm and for daytimes 13–15 UTC and 04–06 UTC UTC
<0.1
0.1–1.0
1–3
3–5
5–10
10–15
15–20
P20
13–15
OBS REG
54.30 48.47
23.22 25.65
14.05 17.77
4.58 4.57
2.58 2.65
0.68 0.64
0.28 0.17
0.33 0.07
04–06
OBS REG
49.15 46.14
21.88 24.20
19.24 20.04
6.53 6.71
2.97 2.85
0.21 0.05
0.02 0.00
0.00 0.00
210
Z. Sokol
Figure 8 Example of three 1-h forecasts by REG (label FOR) and ADV for the 13th July, 2002, 16, 17 and 18 UTC. OBS shows the observed precipitation fields. The size of each square is 9 km · 9 km. Table 5
Evaluation of 1-h forecasts by RMSE, BIAS and CC for the given terms in dependence on observed precipitation amounts (PP)
Date
PP P 0 mm
PP P 1 mm
PP P 5 mm
PP P 10 mm
BIAS
RMSE
CC
RMSE
CC
RMSE
CC
RMSE
CC
02071316
ADV REG
1.05 1.27
5.27 3.00
0.29 0.66
11.01 5.96
0.05 0.51
13.20 8.63
0.01 0.33
16.15 11.33
0.11 0.37
02071317
ADV REG
0.87 1.13
5.26 3.45
0.22 0.59
8.70 6.48
0.04 0.31
11.96 9.11
0.02 0.11
15.68 12.10
0.12 0.12
02071318
ADV REG
1.02 1.06
4.27 2.81
0.46 0.68
7.44 4.73
0.20 0.46
8.51 6.15
0.07 0.26
11.70 9.15
0.18 0.14
Table 6 Number of cases in % when the distance D (in squares) between the local maximum in observed precipitation fields with value at least 10 mm and the nearest local maximum in forecast fields by REG and ADV model is less or equal 0, 1 and 2 UTC
01–03 04–06 07–09 10–12 13–15 16–18 19–21 21–24
REG
ADV
D=0
D61
D62
D=0
D61
D62
19.7 15.6 5.3 5.6 19.4 12.7 12.1 9.1
62.3 62.5 42.1 42.6 60.2 59.3 58.6 52.3
90.2 90.6 84.2 75.9 88.9 86.4 86.2 86.4
6.6 15.6 5.3 5.6 13.9 3.4 6.0 6.8
45.9 46.9 36.8 50.0 64.8 46.6 37.9 56.8
85.2 96.9 68.4 74.1 78.7 66.1 67.2 86.4
The distance D is calculated as the maximum difference in the coordinates x and y.
Nowcasting of 1-h precipitation using radar and NWP data
211
Conclusion
References
Simple advection of precipitation fields, which does not include evolution of precipitation itself, is not appropriate tool for forecasting of 1-h precipitation in the warm season. The proposed model REG, consisting of a regression model and the correction procedure, yields significantly better forecasts than the simple advection. The important part of the model is the correction procedure, which improves the forecast of large precipitation amounts while it does not significantly affect the forecast accuracy for small precipitation. REG is not able to forecast an initiation and development of precipitation in case that no precipitation is predicted by the simple advection forecast. Those cases can be observed quite often during afternoon hours, and therefore, the forecast accuracy significantly deteriorates. During nights and mornings that situation is less frequent, and is usually associated with rather small precipitation amounts. Therefore, the accuracy of the REG forecast is not significantly reduced. This study indicates that the radar and NWP model data combined here do not contain enough information to allow an accurate forecast of precipitation evolution. Therefore, we will turn our attention to a possibility of an inclusion of satellite data into the model. Also retrieved wind fields from tracking of radar echoes can provide a better guide to movements of convective cells which may contribute to the forecast quality.
Charba, J.P., 1998. The LAMP QPF product. Part I: Model development. Weather Forecasting 13, 934–965. Dixon, M., Wiener, G., 1993. TITAN: Thunderstorm identification, tracking, analysis and nowcasting—a radar-based methodology. J. Atmos. Ocean. Tech. 10, 785–797. Elstner, J.B., Schmertmann, C.P., 1994. Assessing forecast skill through cross validation. Weather Forecasting 9, 619–624. Golding, B.W., 1998. Nimrod: A system for generating automated very short range forecasts. Meterorol. Appl. 5, 1–16. Glahn, H.R., Lowry, D.A., 1972. The use of model output statistics (MOS) in objective weather forecasting. J. Appl. Meteor. 11, 1203–1211. Johnson, J.T., MacKeen, P.L., Witt, A., Mitchell, E.D., Stumpf, G.J., Eilts, M.D., Thomas, K.W., 1993. The storm cell identification and tracking algorithm: an enhanced WSR-88D algorithm. Weather Forecasting 13, 263–276. Kitzmiller, D.H., 1996. One-hour forecasts of radar-estimated rainfall by an extrapolative-statistical method. TDL Office Note 96-1, National Weather Service, NOAA, US Department of Commerce, 26 pp. Available from: NWS/TDL, W/OSD2, 1325 East West Highway, Silver Spring, MD 20910. Kra ´ˇ cmar, J., Joss, J., Nova ´k, P., Havra ´nek, P., ˇSa ´lek, M., 1998. First steps towards quantitative usage of data from Czech weather radar network. COST75 Advanced Weather Radar Systems. COST75 Final Seminar, Locarno, Switzerland, 23–27 March 1998, 91–101. Michelson, D., Einfalt, T., Holleman, I., Gjertsen, U., Friedrich, K., Haase, G., Lindskog, M., Szturc, J., 2004. Weather radar data quality in Europe: quality control and characterization. Document COST717. Available from:
Acknowledgement The work was supported by the GA ASCR under the grant S3042101 and by GACR under the grant 205/04/0114. The data were kindly provided by the Czech Hydrometeorological Institute.
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jhydrol
Nowcasting of 1-h precipitation using radar and NWP data Zbynek Sokol The Academy of Sciences of the Czech Republic, Inst. of Atmospheric Physics, Bocni II, 1401, 141 31 Prague, Czech Republic Received 23 March 2005; received in revised form 4 October 2005; accepted 17 December 2005
KEYWORDS
Summary A multiple linear regression model complemented by a correction procedure (REG) is applied to nowcasting of 1-h precipitation for warm seasons. The model aims at predicting the mean area precipitation in the 9 km · 9 km squares for the territory of the Czech Republic. Precipitation amounts are calculated from the radar reflectivity measured by two C band radars, which are operated by the Czech Hydrometeorological Institute. The model predictors are obtained from: (i) radar-derived precipitation in the squares; (ii) radar-derived precipitation advected by wind fields at the 700 hPa level produced by ALADIN/LACE NWP model forecast; (iii) variables derived from forecasts of the ALADIN/LACE NWP model. Accuracy of the REG forecast is evaluated and compared with the forecast obtained by the advection of precipitation fields by the NWP wind at the 700 hPa level (ADV). The ADV forecast serves as one of the predictors (ii). Results show that REG yields apparently better forecasts than ADV in terms of RMSE, correlation coefficients and categorical measures. Also local precipitation maxima have better positions with respect to the real ones. However, if ADV did not indicate any precipitation, REG was not capable to forecast a convective precipitation development. This fact deteriorates the REG accuracy during afternoon hours. We may say that the REG forecast is more accurate in night and morning hours since the situation mentioned occurs quite scarcely during those time periods. c 2006 Elsevier B.V. All rights reserved.
Precipitation forecast; Regression models; Nowcasting; Radar
Introduction
For some hydrological applications, especially the flash flood warning and urban drainage management, make frequently demands of the 0–6 h forecasting (nowcasting) of heavy precipitation. At present, many nowcasting methods employ extrapolations of radar echo. Methods based on a E-mail address: [email protected].
simple extrapolation technique only (e.g. Dixon and Wiener, 1993; Johnson et al., 1993; Mecklenburg et al., 2000) showed a rapidly decreasing accuracy with the increasing lead time. As they do not consider a storm initiation, growth and dissipation the simple extrapolation is usually applicable up to 30–60 min. Radar data including an extrapolated radar echo together with other data provided by satellites, lightning and numerical weather prediction (NWP) models were used
0022-1694/$ - see front matter c 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2005.12.023
Nowcasting of 1-h precipitation using radar and NWP data as predictors in a statistical model by Kitzmiller (1996). Probabilities of exceedence a defined threshold were forecasted by multivariable linear regression models for individual areas of 40 km · 40 km. The regression models were developed on the basis of real data by the application of a technique very similar to the model output statistics (MOS; Glahn and Lowry, 1972). Models for different time intervals during a day were developed separately in order to respect the diurnal cycle of precipitation processes. Another approach to the nowcasting of heavy precipitation is represented by the model NIMROD (Golding, 1998), which utilises both radar echo extrapolation and precipitation forecasts of a NWP model. For shorter time periods the major weight is given to the extrapolation of an existing precipitating field, which is derived from satellite and radar data. With increasing forecast time the weight shifts to the numerical model and becomes almost totally dependent on the NWP model for the 6-h forecast. A different approach is presented by an automated nowcasting system of convective precipitation GANDOLF (Pierce et al., 2000) that attempts to simulate a storm development. It contains a precipitation model, which incorporates a conceptual model of the life cycle of shower clouds. With help of satellite and radar data and a variety of forecast products from a NWP model, convective cells are identified and their possible development is predicted. In this paper, a nowcasting technique based on the application of a linear regression model, similar to that used by Kitzmiller (1996), is applied to the forecast of mean 1-h precipitation in the 9 km · 9 km squares. The forecast is focused on large precipitation amounts, which are important for hydrological applications. Main differences between the applied method and Kitzmiller’s approach are the forecast variables and a significantly higher horizontal resolution. The paper is divided into sections as follows. The next section describes data and a method calculating precipitation amounts from the radar data. Also both forecast and data regions and predictands are defined there. The ensuing section show types of the predictors and an advection algorithm used to derive them. A regression model and a posterior correction of model outputs are given next. The method of model verification and results and discussion are also given. The conclusions are outlined in the final section.
Data and forecast region The data used in this study are from the 2002 warm season and include radar-derived precipitation and prognostic fields of the NWP model ALADIN/LACE (Aire Limite ´e Adaptation Dynamique De ´veloppement InterNational/Limited Area Modelling in Central Europe), which is routinely applied by the Czech Hydrometeorological Institute (CHMI). Fig. 1 shows the data and forecast regions. The radar measurements collected from the data region are used. This region is larger than the forecast one because the radar data are advected by the NWP model wind. Both regions are well covered by radar and NWP model data. Precipitation amounts are derived from radar reflectivity measurements. Radar-derived precipitation cannot be adjusted by ground measurements as there are no available hourly measurements from rain gauges reliable enough. Ra-
201 dar data are recorded by the Czech radar network consisting of two weather radars Brdy and Skalky. Both radars measure in C band and the interval of their scans is 10 min (Fig. 1). Precipitation amounts are calculated by integration of rain rates in time. The rain rates are derived from the radar reflectivity at the elevation of 2 km above the sea level (CAPPI 2 km) with help of a standard algorithm operationally used in CHMI (Zacharov et al., 2004; Michelson et al., 2004). The bright band is not considered here, as it occurs very rarely at or below the elevation of 2 km in warm seasons. Rain rates are calculated by using the Z–R relationship Z = 200R1.6 in the range of 7 dBZ 6Z 6 55 dBZ. For Z < 7 dBZ and Z > 55 dBZ the rates are R = 0 mm/h and R = 99.85 mm/ h, respectively. For both radars the maximum range of rainfall estimation is 256 km and the azimuth resolution is 1. Data from Skalky radar are recorded in 20 elevation angles and the ray bin resolution is 1 km. Radar Brdy uses 14 elevation angles and ray bin resolution of 0.5 km. The raw radar coordinates are transformed into the gnomonic projection and rain rates are interpolated/extrapolated separately for both radars into 512 · 512 pixels with the horizontal resolution of 1 km · 1 km. The rain rates in pixels outside the 256 km range are set to zero. The pixel values are smoothed by using the median from 3 · 3 nearby pixels, which improves radar precipitation when compared with gauge measurements (Kra cmar et al., 1998). Then 20-min precipitation ´ˇ amounts are calculated in the pixels and the maximum of two radar fields is applied in the overlapping pixels to obtain precipitation estimates in the data region. The integration starts at 00, 20 and 40 min for each hour. The data region is covered by squares consisting of 9 · 9 radar pixels. For these squares the mean pixel values of 1-h precipitation amounts are calculated as the sums of corresponding 20-min amounts. The square size selected allows the forecast of precipitation even for quite small river basins, the area of which is about 100 km2. The 9 km · 9 km squares in the forecast region are a subset of the squares covering the data region. The forecast variables (predictands) are 1-h precipitation amounts in squares of the forecast region, consisting of 43 · 25 squares and covering a large part of the territory of the Czech Republic. For simplicity reasons, we will be using the term ‘‘precipitation amount’’ instead of mean area precipitation amount. Outputs from the NWP model ALADIN/LACE include predicted temperature, relative humidity, geopotential, horizontal wind and vertical velocity at levels of 925, 850, 700, 500 and 300 hPa, and temperature, relative humidity, wind at the surface, and the sea level pressure. ALADIN/ LACE NWP model is a hydrostatic limited area model with the horizontal resolution of 12 km. The model integration starts every 12 h and basic prognostic fields are recorded every 6 h. In order to obtain the prognostic values in other times the linear interpolation between two consecutive 6hourly terms is used.
Predictors All predictors are available no later than at time T and the forecast 1-h precipitation is related to the interval [T, T + 60 (min)]. A multiple linear regression employs four
202
Z. Sokol
Figure 1 Data and forecast regions. The coordinates are in radar pixels, the size of which is 1 km · 1 km. The squares show the radar locations, and circle parts show the maximum radar range 256 km.
types of predictors. The first type consists of 1-h precipitation amounts and instantaneous rain rates. Precipitation predictors are determined in the fixed squares as well as in the squares obtained by the advection of precipitation fields. Besides the squares of 9 · 9 pixels, which correspond to the predictand squares, also squares of 17 · 17 and 25 · 25, which have the same centres as the 9 · 9 squares, are considered. The advection makes use of a prognostic wind field at the 700 hPa level and a backward scheme with the time step of 20 min. We considered also a wind field at 500 hPa level and a mean wind field in the layer of 700– 500 hPa. We forecasted 1-h precipitation amounts by advection of the actual precipitation field by various wind fields and compared results. We found out that the forecast by the wind at the 700 hPa level was slightly more accurate than forecasts employing other wind fields. The advection starts at centres of the 9 · 9 squares in the forecast domain (Fig. 1). After each advection time step, 20min precipitation is determined in the squares centred in the pixels resulting from the advection steps. In order to determine accumulated precipitation over interval [t 20, t] we apply the prognostic wind corresponding to the t 10 time. The predictor value of 1-h precipitation is the sum of 20-min precipitation in corresponding squares. The application of three 20-min steps turned out to be slightly more accurate than one 60-min step. The following predictors are determined for times t = T, T 20, T 40 and T 60:
• Current mean and maximum 1-h precipitation in the squares of 9 · 9, 17 · 17 and 25 · 25. • Current mean and maximum instantaneous precipitation rates in the squares of 9 · 9, 17 · 17 and 25 · 25. • Forecasted mean and maximum 1-h precipitation in the squares of 9 · 9, 17 · 17 and 25 · 25 that are produced by the advection algorithm. Mean precipitation in the squares of 9 · 9 for t = T represents the forecast by advection only. • Forecasted mean and maximum instantaneous precipitation rates in the squares of 9 · 9, 17 · 17 and 25 · 25 obtained by the application of the advection algorithm. • Mean vertically integrated liquid water content (VIL) in the squares of 9 · 9. • Mean top of radar echoes (ECHOTOP) in the squares of 9 · 9. • Mean VIL and ECHOTOP in the squares of 9 · 9 obtained by the advection algorithm. The predictors of the second type are prognostic NWP model outputs, which are interpolated in space and time to get predictor values for the individual squares and forecast times. In addition, the following derived parameters are also used: K-index, vertical gradient of equivalent potential temperature in the layers of 700–500 and 850– 700 hPa, the mean vertical relative humidity and vertical velocity (in the layer from 925 to 500 hPa). The list of
Nowcasting of 1-h precipitation using radar and NWP data predictors is complemented by the tendencies of the NWP model predictors over 6-h intervals containing time T. The third type of predictors contains linearized predictors calculated in a way similar to that used by Charba (1998). The idea is to linearize relationships between the predictors and predictands. Predictor values are sorted in ascending order and divided into 10 groups with an approximately equal number of members. In each group mean predictor values, (Xi, i = 1, . . . , 10) and mean predictand values, (Yi, i = 1, . . . , 10) are calculated and Yi are treated as a piecewise linear function of Xi. Linearized predictor values Y(x) are calculated as follows: YðxÞ ¼ Y i
x X iþ1 X iþ1 x þ Y iþ1 ; X iþ1 X i X iþ1 X i
for X i 6 x 6 X iþ1 ; ð1Þ
YðxÞ ¼ X 1 ;
for x 6 X 1 ;
ð2Þ
203 YðxÞ ¼ X 10 ;
for X 10 6 x;
ð3Þ
where x is the observed predictor value. The grouping mentioned above follows from results of a number of preliminary tests. The fourth type of predictors has the form of (a Æ b)0.5, where a and b are predictors of the first type.
Regression models We use a linear model as many others have done in spite of the fact that the relationships between predictors and predictands are supposed to be non-linear. Those relationships seem to be complex and predictors probably do not contain sufficient information about precipitation processes to allow an accurate forecast. Moreover, values of predictors are supposed to be influenced by various types of errors. Therefore, we use a simple linear model, robust and easily applicable, as an approximation of the relationships.
Figure 2 Scatter plots (c) and (d) illustrate effects of the correction procedure for M = 5 and M = 10 values, respectively, on the standard regression model output (a) M = 0. The plot (b) shows the areas from which the pi and yi values are used in calculating corrected output pcor,i. For M = 10 the data from the darkest area are used. For M = 5 are used data from the union of the darkest and light grey areas. The data concerns a calibration subset for 22–24 UTC.
204
Z. Sokol
Table 1 Mean predictand values (Mean) and probabilities (%) that predictand values will reach or exceed thresholds P = 1, 3, 5 and 10 mm for 8 considered daytimes (UTC) UTC
Mean (mm)
P = 1 mm
P = 3 mm
P = 5 mm
P = 10 mm
01–03 04–06 07–09 10–12 13–15 16–18 19–21 22–24
0.85 0.93 0.72 0.68 0.93 0.94 1.16 0.81
24.2 29.0 22.9 18.6 22.5 22.6 27.5 24.7
8.5 9.7 6.4 6.1 8.4 8.8 11.8 7.2
3.4 3.2 2.1 2.6 3.9 4.5 6.0 2.8
0.6 0.2 0.2 0.5 1.3 1.4 1.7 0.4
The predictand is 1-h precipitation amount in the forecast squares 9 km · 9 km.
Following the scheme suggested by Kitzmiller (1996) we have developed separately 8 multiple linear regression models for the forecast times 01–03 UTC, 04–06
UTC, . . . , and 21–24 UTC. The models take a linear form such as y ¼ a0 þ a1 x 1 þ a2 x 2 þ þ aN x N ;
ð4Þ
where a1, . . . , aN are model parameters, y is the predictand value and xj, j = 1, . . . , N, are predictor values. The parameters are found by minimizing the sum S of squared differences between observed yi and predicted values pi: S¼
nc X ðy i pi Þ2 ;
ð5Þ
i¼1
where pi values are obtained with help of formula (4), and they are pi ¼ a0 þ a1 z1;i þ a2 z2;i þ þ aN zN;i ;
ð6Þ
z1,i, . . . , zN,i are predictor values corresponding to the ith predictand value and nc is the number of predictands used to derive the model parameters. The solution of the minimization problem can be found in standard regression texts e.g. Wilks (1995).
Figure 3 Scatter plots (a)–(c) show REG forecasts in dependence on M. The forecast is calculated for 16–18 UTC. In (d) the solid lines show average forecast values for the defined observation categories and dashed lines show average observation values for the defined forecast categories in dependence on M. See the text for details.
Nowcasting of 1-h precipitation using radar and NWP data The predictors are selected by the forward-selection screening algorithm (Wilks, 1995). The selection stops when the decrease of mean-square-error (MSE) by the inclusion of a new predictor is smaller than 0.2% of the previous MSE. The maximum number of predictors was restricted to 16 as the tests showed that additional predictors did not improve the forecast. The choice of the 0.2% threshold was subjectively made with respect to the preliminary tests. The model outputs are smoother than the observed precipitation and significantly underestimate heavy precipitation due to its scarcity within the dataset (for the confirmation see Fig. 2(a). The underestimation of heavy precipitation can be corrected by a posterior correction of model outputs pi by the linear transformation: pcor;i ¼ api þ b;
i ¼ 1; . . . ; nc ;
ð7Þ
where a and b are parameters. Their meaning will be elucidated in the following text.
205 Parameter a is found to minimize the function X FðaÞ ¼ ðy i api Þ2 .
ð8Þ
i;y i ;pi >M
The sum in (8) represents all data i from the calibration data set, whenever yi (observations) or pi (model outputs) are greater than a given value of M. An example of the areas from which the data are taken for M = 5 and M = 10 are shown in Fig. 2(b). In order to determine a for M = 10, only predictor values and corresponding observations from the darkest area are used. In case of M = 5 the light grey area is added to the darkest one. These areas contain mainly those forecasts, in which the model significantly underestimates observations. Therefore, the parameter a is greater than 1 and consequently, api > pi. Parameter b aims at maintaining an approximate equality between the sum of observations and corrected model outputs pcor,i, which yields, when (5) and (7) are employed:
Table 2 RMSE and CC for REG with M = 0, M = 5 and M = 10 applied to 16–18 UTC in dependence on observed precipitation amounts (PP) PP P 0 mm
M=0 M=5 M = 10
PP P 1 mm
PP P 5 mm
PP P 10 mm
BIAS
RMSE
CC
RMSE
CC
RMSE
CC
RMSE
CC
0.96 0.97 0.99
1.70 1.71 1.82
0.68 0.68 0.68
3.33 3.32 3.51
0.54 0.55 0.55
6.83 6.57 6.45
0.34 0.34 0.34
10.55 9.96 9.41
0.24 0.24 0.24
BIAS is calculated only for all data.
Table 3 UTC
RMSE and CC for ADV and REG methods in dependence on daytimes (UTC) and observed precipitation amounts (PP) Method
PP P 0 mm
PP P 1 mm
PP P 5 mm
BIAS
RMSE
CC
RMSE
CC
RMSE
PP P 10 mm CC
RMSE
CC
01–03
ADV REG
0.94 1.01
1.49 1.12
0.61 0.77
2.75 2.06
0.40 0.61
5.64 4.29
0.14 0.35
10.03 7.70
0.07 0.34
04–06
ADV REG
0.99 1.02
1.28 0.91
0.70 0.83
2.12 1.53
0.51 0.66
3.84 2.91
0.23 0.30
6.66 5.87
0.24 0.14
07–09
ADV REG
1.11 0.99
1.17 0.85
0.65 0.80
2.05 1.55
0.43 0.63
4.37 3.41
0.18 0.25
8.30 7.07
0.16 0.24
10–12
ADV REG
0.93 0.94
1.44 1.10
0.57 0.74
3.08 2.35
0.35 0.60
6.63 5.14
0.08 0.33
11.35 8.97
0.06 0.16
13–15
ADV REG
0.96 0.99
2.58 1.76
0.41 0.70
5.02 3.44
0.24 0.60
10.26 7.45
0.05 0.40
15.72 11.72
0.03 0.19
16–18
ADV REG
0.98 0.97
2.69 1.71
0.37 0.68
5.09 3.32
0.18 0.55
9.26 6.54
0.02 0.34
13.70 9.96
0.04 0.24
19–21
ADV REG
0.98 0.99
2.53 1.64
0.50 0.76
4.40 2.89
0.32 0.64
7.73 5.27
0.06 0.43
11.61 8.15
0.10 0.29
22–24
ADV REG
1.12 1.01
1.47 1.03
0.62 0.77
2.56 1.86
0.46 0.62
5.07 3.98
0.25 0.31
7.80 6.89
0.06 0.10
BIAS is calculated only for all data.
206 bð0Þ ¼
Z. Sokol nc 1 X ðy api Þ. nc i¼1 i
ð9Þ
The calculated values of a and b(0) may lead to negative values of pcor,i. Thus, the following iterative process (k = 1, . . . , kmax) is used in calculating b: ðkÞ ð10Þ pcor;i ¼ max api þ bðkÞ ; 0 ; i ¼ 1; . . . ; nc ; bðkþ1Þ ¼ bðkÞ þ
nc 1 X ðkÞ y i pcor;i ; nc i¼1
ð11Þ
It converges in several steps. The test showed that kmax = 5 is sufficient and b ¼ bðkmax Þ . Even after iterations pcor,i can be negative. In this case pcor,i is set to zero. Regression model (4) with the correction algorithm ((7)– (11)) is referred to as REG. The correction (7) is used for all regression outputs pi . The influence of the M value on the model outputs is illustrated by Fig. 2. It shows that the correction improves the forecast of heavy precipitation. On the other hand, it may overestimate low precipitation amounts.
The M values used in this study were derived at from the test results.
Calibration and verification of the models For each of 8 daytimes, 39 terms with large local precipitation are used for model calibration and verification. It means that 41,925 (39 · 43 · 25) predictand values are available for each model. Basic statistical characteristics of the predictand values are shown in Table 1. The majority of large local precipitation, frequently occurring in afternoons in the warm season, is of convective origin, and the forecast of the convective precipitation is very difficult. For each daytime the data are divided into 3 subsets, each containing 13 terms. The data are placed in the subsets by time sequence to maintain effects of synoptic regimes. Two subsets are used as the calibration data and the third one serves for the verification. A cross-validation technique (Elstner and Schmertmann, 1994), choosing different combinations of the subsets, is applied to get the independent
Figure 4 Scatter plots of forecast and observed values for 13–15 UTC. (a)–(c) show REG forecasts: (a) shows all the values, (b) shows only pairs when ADV forecasts positive precipitation, and (c) depicts events when ADV is equal to 0. (d) shows ADV forecast for all data.
Nowcasting of 1-h precipitation using radar and NWP data forecasts for all data. The independent forecasts are used in the following verification tests. The arrangement of the data into the subsets puts close terms, which may reflect similar development of precipitation, in the same subset in most cases. Therefore, they do not have a substantial influence on verification results. The accuracy of the forecast is evaluated by the rootmean-square-error (RMSE), bias (BIAS) and correlation coefficient (CC): sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 1X RMSE ¼ ðy pi Þ2 ; n i¼1 i
ð12Þ
Pn pi ; BIAS ¼ Pi¼1 n i¼1 y i
ð13Þ
Pn
P P y i n1 nj¼1 y j pi n1 nj¼1 pj CC ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 P 2ffi ; Pn Pn Pn n 1 1 y y p p i i i¼1 j¼1 j i¼1 j¼1 j n n i¼1
Figure 5
ð14Þ
207 where pi, yi are forecast and observed values of predictands. As the forecast is primarily aimed at large precipitation amounts, RMSE and CC are calculated with respect to the values of actual precipitation yi. The model results are compared with the forecast obtained by the advection of the precipitation field by the NWP wind at 700 hPa (ADV), which is one of the predictors.
Results and discussion The selected predictors and a number of them differ for daytimes and among the particular calibration data sets. Large variability of the selected predictors results from the fact that the majority of the predictors are highly correlated and therefore, there are various groups of selected predictors that produce comparable results. The screening algorithm yields numbers of predictors in the range from 12 to 16. Even higher number of predictors does not cause an overfitting but also does not improve the model accuracy.
The same as Fig. 3 for 04–06 UTC.
208
Figure 6 Relations between the forecast and observed values by REG and ADV methods for 13–15 UTC (a) and for 04–06 UTC (b). The solid lines show average forecast values for defined observation categories and the dashed lines show average observation values for defined forecast categories for M = 5. See the text for details.
The most important predictors belong to the first and third types. Predictors derived from the NWP model are contained almost in each selection but they have little impact on the forecast accuracy. The same is true for the fourth type of predictors. Although the majority of selected predictors are valid at time T, when the forecast is issued, predictors from earlier times are also frequently included. There is a difference between the most frequently selected predictors for afternoon and other daytimes. For the afternoon (13–15 UTC, 16–18 UTC) the most frequent predictors are: the mean precipitation rates in 9 · 9 pixels at time T and the maximum precipitation rates in 25 · 25 pixels at time T. Other predictors are related to precipitation or precipitation rates and squares 9 · 9 pixels. Some linearized predictors are also selected but not so often, e.g. advected VIL. For the remaining daytimes two predictors: the advected mean precipitation rates in 25 · 25 pixels
Z. Sokol at time T and the linearized advected ECHOTOP at time T are more frequently chosen than others. Other often selected predictors are the maximum precipitation rate in 9 · 9 pixels and advected VIL. The influence of the value of M (8) is demonstrated by Fig. 3 and Table 2 for the interval 16–18 UTC. Three scatter plots of observed and forecast values for M = 0 (standard regression), M = 5 and M = 10 show that REG is able to forecast higher values for larger M. Table 2 shows that with increasing M the accuracy of the forecast increases for higher measured precipitation. On the contrary, the large value of M is reducing accuracy for small precipitation amounts. It is worth mentioning that BIAS is small for all values of M. It proves that the iterative algorithm ((9)–(11)) is helpful in reducing BIAS values. Fig. 3(d) compares the forecasts for different M in a more transparent way. The horizontal coordinates of marks on full lines are calculated as the mean of observed values from intervals [0, 0.5], [0.5, 3.5], [3.5, 6.5], . . . , [21.5, 24.5] and P24.5 mm. The vertical coordinates are the means of forecast values corresponding to observed data from the relevant intervals. The dashed lines contain the marks calculated in a similar way, but at first, forecast values are averaged over the intervals. The standard regression should yield a dashed line close to the line representing the perfect fit. As it follows from Table 2 and Fig. 3 the REG model with M = 5 apparently improves the forecast of large precipitation in comparison with M = 0 and only negligibly decreases the accuracy for small precipitation amounts. REG with M = 10 even more improves the forecast of large precipitation; however, it deteriorates forecasts of small precipitation. Consequently, REG with M = 5 is used in the following calculations. The quality of the forecast depends essentially on the daytime. The forecast values are apparently less accurate (Table 3) and more scattered around the perfect fit for the terms 13–15 and 16–18 UTC than for the remaining time intervals. It is illustrated by Figs. 4–6 for 04–06 UTC and 13–15 UTC. The comparison of Figs. 4(a) and 5(a) with Figs. 4(d) and 5(d) show that REG significantly improves the forecast of ADV. For afternoon the simple advection ADV is definitely unsuitable for the precipitation forecasting. From Fig. 4(c), which shows the forecast for those events in which the forecast by ADV is equal to 0 mm, follows that REG is not capable to produce new convective cells. When these cases are excluded from the data then the REG forecast gives better results (Fig. 4(b)). The mentioned drawback of REG is not very important for morning and night terms (Fig. 5) as a development of new storm cells is rather scarce. Fig. 6 compares REG and ADV forecasts by using the means similar to those in Fig. 3(d). Also Fig. 7, where the quality of the forecast for 13–15 UTC and 04–06 UTC is compared in terms of categorical measures POD (probability of detection), FAR (false alarm rate) and ETS (equitable thread score; Wilks, 1995), confirms a lower accuracy of afternoon forecasts except for the 99% percentile. The continuous forecast is transformed into a binary forecast using thresholds corresponding to the 50%, 70%, 80%, 90%, 95% and 99% percentiles of the observed precipitation. Another drawback of the forecasts common for all terms is the tendency to smooth maximum and minimum values. This is due to minimization of mean-square-error in the regression. Histograms in Table 4 show that REG forecasts
Nowcasting of 1-h precipitation using radar and NWP data both low (<1 mm) and high values with lower frequencies than they are observed. This behaviour of REG can be noted in Fig. 8, where examples of three consecutive forecasts of extremely heavy precipitation are shown. The REG fields are smoother than the observed or advected precipitation fields; however, REG describes shapes of areas with precipitation better than ADV. Heavy precipitation forecasts provided by ADV have significant location errors, which lead to a larger RMSE and lower CC, as shown in Table 5. We tried to find out how well REG can localize large precipitation amounts. Therefore, we studied differences between the positions of observed local maxima and positions of local maxima in forecasts produced by REG and ADV. We consider all inner squares within the forecast domain with observed precipitation at least 10 mm under
209 the condition that 8 neighbouring values are lower (local maximum). For these squares we look for the nearest local maxima in the forecast fields and evaluate distances (in terms of square numbers) between the observed and forecast local maxima. Table 6 shows that the accuracy of the forecast by REG and ADV differs for individual daytimes. For REG about 12% of local maxima are localized well, 55% are localized correctly or in neighbouring squares and in 86% of all cases distances between the observed and forecast local maxima do not exceed the distance of 2 squares (18 km). For ADV the localization is less accurate: 8%, 48% and 78%. The good results of REG for afternoon terms 13– 15 UTC and 16–18 UTC are worth noting. Smooth stratiform precipitation fields, where local maxima are not well pronounced, might influence the evaluation for night terms.
Figure 7 Evaluation of REG forecasts for 13–15 UTC and 04–06 UTC by POD, FAR and ETS. The horizontal axis shows the 50%, 70%, 80%, 90%, 95% and 99% percentiles of observed precipitation. The corresponding values of precipitation amounts are 0.04, 0.58, 1.19, 2.62, 4.29, 11.84 mm and 0.12, 0.94, 1.70, 2.95, 4.19, 7.31 mm for 13–16 UTC and 4–6 UTC, respectively.
Table 4 Frequency in % of observed (OBS) and forecast (REG) values for the given precipitation intervals in mm and for daytimes 13–15 UTC and 04–06 UTC UTC
<0.1
0.1–1.0
1–3
3–5
5–10
10–15
15–20
P20
13–15
OBS REG
54.30 48.47
23.22 25.65
14.05 17.77
4.58 4.57
2.58 2.65
0.68 0.64
0.28 0.17
0.33 0.07
04–06
OBS REG
49.15 46.14
21.88 24.20
19.24 20.04
6.53 6.71
2.97 2.85
0.21 0.05
0.02 0.00
0.00 0.00
210
Z. Sokol
Figure 8 Example of three 1-h forecasts by REG (label FOR) and ADV for the 13th July, 2002, 16, 17 and 18 UTC. OBS shows the observed precipitation fields. The size of each square is 9 km · 9 km. Table 5
Evaluation of 1-h forecasts by RMSE, BIAS and CC for the given terms in dependence on observed precipitation amounts (PP)
Date
PP P 0 mm
PP P 1 mm
PP P 5 mm
PP P 10 mm
BIAS
RMSE
CC
RMSE
CC
RMSE
CC
RMSE
CC
02071316
ADV REG
1.05 1.27
5.27 3.00
0.29 0.66
11.01 5.96
0.05 0.51
13.20 8.63
0.01 0.33
16.15 11.33
0.11 0.37
02071317
ADV REG
0.87 1.13
5.26 3.45
0.22 0.59
8.70 6.48
0.04 0.31
11.96 9.11
0.02 0.11
15.68 12.10
0.12 0.12
02071318
ADV REG
1.02 1.06
4.27 2.81
0.46 0.68
7.44 4.73
0.20 0.46
8.51 6.15
0.07 0.26
11.70 9.15
0.18 0.14
Table 6 Number of cases in % when the distance D (in squares) between the local maximum in observed precipitation fields with value at least 10 mm and the nearest local maximum in forecast fields by REG and ADV model is less or equal 0, 1 and 2 UTC
01–03 04–06 07–09 10–12 13–15 16–18 19–21 21–24
REG
ADV
D=0
D61
D62
D=0
D61
D62
19.7 15.6 5.3 5.6 19.4 12.7 12.1 9.1
62.3 62.5 42.1 42.6 60.2 59.3 58.6 52.3
90.2 90.6 84.2 75.9 88.9 86.4 86.2 86.4
6.6 15.6 5.3 5.6 13.9 3.4 6.0 6.8
45.9 46.9 36.8 50.0 64.8 46.6 37.9 56.8
85.2 96.9 68.4 74.1 78.7 66.1 67.2 86.4
The distance D is calculated as the maximum difference in the coordinates x and y.
Nowcasting of 1-h precipitation using radar and NWP data
211
Conclusion
References
Simple advection of precipitation fields, which does not include evolution of precipitation itself, is not appropriate tool for forecasting of 1-h precipitation in the warm season. The proposed model REG, consisting of a regression model and the correction procedure, yields significantly better forecasts than the simple advection. The important part of the model is the correction procedure, which improves the forecast of large precipitation amounts while it does not significantly affect the forecast accuracy for small precipitation. REG is not able to forecast an initiation and development of precipitation in case that no precipitation is predicted by the simple advection forecast. Those cases can be observed quite often during afternoon hours, and therefore, the forecast accuracy significantly deteriorates. During nights and mornings that situation is less frequent, and is usually associated with rather small precipitation amounts. Therefore, the accuracy of the REG forecast is not significantly reduced. This study indicates that the radar and NWP model data combined here do not contain enough information to allow an accurate forecast of precipitation evolution. Therefore, we will turn our attention to a possibility of an inclusion of satellite data into the model. Also retrieved wind fields from tracking of radar echoes can provide a better guide to movements of convective cells which may contribute to the forecast quality.
Charba, J.P., 1998. The LAMP QPF product. Part I: Model development. Weather Forecasting 13, 934–965. Dixon, M., Wiener, G., 1993. TITAN: Thunderstorm identification, tracking, analysis and nowcasting—a radar-based methodology. J. Atmos. Ocean. Tech. 10, 785–797. Elstner, J.B., Schmertmann, C.P., 1994. Assessing forecast skill through cross validation. Weather Forecasting 9, 619–624. Golding, B.W., 1998. Nimrod: A system for generating automated very short range forecasts. Meterorol. Appl. 5, 1–16. Glahn, H.R., Lowry, D.A., 1972. The use of model output statistics (MOS) in objective weather forecasting. J. Appl. Meteor. 11, 1203–1211. Johnson, J.T., MacKeen, P.L., Witt, A., Mitchell, E.D., Stumpf, G.J., Eilts, M.D., Thomas, K.W., 1993. The storm cell identification and tracking algorithm: an enhanced WSR-88D algorithm. Weather Forecasting 13, 263–276. Kitzmiller, D.H., 1996. One-hour forecasts of radar-estimated rainfall by an extrapolative-statistical method. TDL Office Note 96-1, National Weather Service, NOAA, US Department of Commerce, 26 pp. Available from: NWS/TDL, W/OSD2, 1325 East West Highway, Silver Spring, MD 20910. Kra ´ˇ cmar, J., Joss, J., Nova ´k, P., Havra ´nek, P., ˇSa ´lek, M., 1998. First steps towards quantitative usage of data from Czech weather radar network. COST75 Advanced Weather Radar Systems. COST75 Final Seminar, Locarno, Switzerland, 23–27 March 1998, 91–101. Michelson, D., Einfalt, T., Holleman, I., Gjertsen, U., Friedrich, K., Haase, G., Lindskog, M., Szturc, J., 2004. Weather radar data quality in Europe: quality control and characterization. Document COST717. Available from:
Acknowledgement The work was supported by the GA ASCR under the grant S3042101 and by GACR under the grant 205/04/0114. The data were kindly provided by the Czech Hydrometeorological Institute.