Precipitation information with an accurate uncertainty estimate is highly relevant for research areas with special interest in areal precipitation information. The uncertainty information is needed for example to generate precipitation ensembles as input for probabilistic hydrological modelling (e.g. Krajewski and Ciach, 2006; Germann et al., 2009; Rico-Ramirez et al., 2015; Park et al., 2016) or nowcasting (e.g. Atencia and Zawadzki, 2015; Dai et al., 2015). Furthermore, precipitation observations are a valuable source of information for data assimilation in numerical weather prediction (e.g. Dowell et al., 2011; Chang et al., 2014; Bick et al., 2016). Recent studies stress the importance of correctly specifying errors in this context (e.g. Waller et al., 2016, 2017).

To assess the reliability of any observation, it is crucial to know its uncertainty. For precipitation observations, accurate measurements as well as a reliable estimate of their uncertainty are an ongoing topic of research. Observations from weather radars are of particular interest because they provide valuable areal precipitation information with high spatial and temporal resolution, in contrast to in situ observations. However, radar measurements are affected by numerous sources of error that diminish the accuracy of the provided precipitation quantification, e.g. calibration errors, noise, interferences, clutter, and attenuation. Despite the use of correction algorithms and filters on radar measurements, residual error remains inherent to the data. In addition, the estimation of rainfall intensity from radar data requires a relation between measured reflectivity and rain rate which is empirical and uncertain. Therefore, accurate rainfall precipitation quantification is not possible with radar measurements only. An extensive review of uncertainty sources in single polarisation radar measurements is given in Villarini and Krajewski (2010) and citations therein.

Improvements in the precipitation estimate are achieved by combining measurements from different sources, e.g. radar and rain gauge data, to take advantage of as much information as possible. Since no observation is error-free, they must be weighted with their respective uncertainty in order to obtain a statistically sound combined precipitation product. Many studies investigate the statistical combination of precipitation data. They rely on computing the optimal estimate of precipitation by minimising its error variance. The majority of the applied techniques focus on a static, i.e. time-independent, merging approach, where precipitation measurements are combined for every available time step independently. Most merging approaches use kriging or cokriging schemes to spatially merge radar and rain gauge data (e.g. Krajewski, 1987; Creutin et al., 1988; Haberlandt, 2007; Berndt et al., 2014). Kriging methods acknowledge the fact that precipitation observation errors are spatially correlated by the specification of variograms. These variograms are modelled based on theoretical considerations, and do not consider the temporal structure of the errors. Less commonly used are variational data assimilation methods that describe the correlation between observation errors through covariance matrices. These methods can easily be extended to include further observations. For example, Bianchi et al. (2013) use a variational approach to estimate precipitation from radar, gauge and microwave link data.

Some studies additionally take into account the temporal correlation of precipitation measurement errors while merging radar and gauge data. The temporal structure of the errors can be considered either by modelling spatio-temporal variograms in the kriging approach (e.g. Sideris et al., 2014, Pulkkinen et al., 2016) or by using Kalman filtering that combines knowledge on error variances from previous time steps for the current weighting of observations (e.g. Smith and Krajewski, 1991; Chumchean et al., 2006). The latter approach is also used by Grum et al. (2005) with the addition of microwave link data. Even though these methods consider temporal error structures in the merging process, the continuous evolution of the precipitation system is not represented.

Very few studies include a forecast model that evolves the system’s state to take advantage of previous knowledge of the system. The coupling of a forecasting model can be done using Kalman filtering or variational data assimilation methods. In addition to the merging of observations, the forecasting component allows for the extrapolation of information to data-void regions. Advection is commonly used as an approximation for the evolution of precipitation systems. Zinevich et al. (2009) and Mercier et al. (2015) present frameworks to estimate areal precipitation from microwave and television satellite links, respectively. Fielding et al. (2014) retrieve three-dimensional cloud properties from cloud radar and radiance observations. These studies integrate the temporal evolution of the system into the statistical merging of observations with a focus on retrieving the optimal combined estimate of the system’s state. However, they do not focus on an assessment of the uncertainty associated with the resulting precipitation field.

Considering radar precipitation data, uncertainty estimation is extensively studied outside the context of combined precipitation products. Numerous studies address the quantification of precipitation data uncertainty to overcome the deterministic view on radar precipitation measurements. The basis for the description of areal precipitation measurement errors mostly is an empirical, statistical study of radar and reference surface measurements. Some methods only provide static information (Krajewski and Ciach, 2006), while more recent approaches also allow for the description of the spatial and temporal structure of the errors through the description of the error covariance (e.g. Ciach et al., 2007; Germann et al., 2009; Villarini and Krajewski, 2009; Dai et al., 2014). Such spatial uncertainty description allows for a probabilistic assessment of precipitation information suitable for applications. However, the resulting error description is neither dynamical nor flow-dependent.

Here, we present a method that connects both aspects: precipitation data merging and probabilistic uncertainty assessment. The method allows for the combination of precipitation observations considering their respective uncertainty and takes advantage of the additional information provided by the temporal evolution of the system. At the same time, the method yields an areal uncertainty estimate for the resulting precipitation product. Because of the included temporal evolution, the uncertainty estimate is variable both in space and time and is flow-dependent. Thus, this method aims at providing both an accurate precipitation product and an improved areal and dynamical uncertainty estimate.

Our approach uses data assimilation as a tool to merge precipitation observations. This study considers radar data with different spatial resolution, but the method can be extended to incorporate any additional source of precipitation observation. Data assimilation techniques allow for statistically combining all available information considering its uncertainty within a temporal evolution model. The Local Ensemble Transform Kalman Filter (LETKF, Bishop et al., 2001; Ott et al., 2004; Hunt et al., 2007) used here is especially suited for the purpose of this study. Its main feature is the description of uncertainty through covariance matrices estimated by an ensemble approach. The additional use of localisation in the LETKF reduces the required number of ensemble members to adequately represent the covariances. Through both the use of an ensemble and the localisation, the LETKF is computationally much more affordable for our problem size than the standard Kalman filter approach (Kalman, 1960; Kalman and Bucy, 1961; Talagrand, 1997). Furthermore, the ensemble approach can handle non-linear forecast models. In the presented setup, the LETKF is coupled to an ensemble nowcasting model to provide information about the displacement of precipitation over time.

The added value of the presented method is demonstrated in four experiments with observations from an X-band radar network. The method’s functioning, including the LETKF and the coupled nowcasting scheme, is presented in Section 2, Section 3 briefly presents the data. The performed experiments are introduced in Section 4 and the obtained precipitation product is evaluated in Section 5.



The combination of a data assimilation method with a nowcasting method introduces flow dependency to the areal uncertainty estimate, as known uncertainty is propagated in time and space. In this way, the method yields a situation-dependent uncertainty estimate whose structure follows the evolution of the system. The overall functioning of the method is outlined in the flowchart in Fig. 1. The ensemble data assimilation cycle allows for obtaining a combined precipitation product and a corresponding uncertainty estimate from the ensemble mean and spread, respectively. The data assimilation scheme (LETKF) and the coupled nowcasting scheme are introduced in the following.

Fig. 1.  

Flowchart illustrating the data assimilation cycle of the method for uncertainty estimation. Combined precipitation product and uncertainty estimate are obtained from the analysis ensemble mean and spread, respectively.



Data assimilation aims at combining different sources of information about an observed system, weighted by their respective uncertainty, in order to get the best possible estimate of the system’s true state. The data assimilation method used here is the LETKF that allows to work with an ensemble of moderate size through a localisation approach. It was first introduced by Hunt et al. (2007) based on previous work by Bishop et al. (2001) and Ott et al. (2004). The complete set of LETKF equations (see Hunt et al., 2007 for a complete derivation), consisting of the analysis and forecast steps, are given as

((1) )
Analysis:x̂ka=x̂bk+Kk[ykoHkx̂bk]Zka=ZkbTkLETKFTkLETKF=Ck[In+Dk]1/2CkTxk,ia=x̂ka+l1Zk,iaForecast:xk+1,ib=M(xk,ia)with Kk=P˜kbHkT[HkP˜kbHkT+Rk]1and Zka,b=1l1(xk,1a,bx̂ka,bxk,2a,bx̂ka,bxk,la,bx̂ka,b),
where Ck is the matrix of eigenvectors of ZkbTHkTRk1HkZkb and Dk is a diagonal matrix with the corresponding eigenvalues. In order to get the analysis ensemble mean after the update, x̂ka, the background ensemble mean before the update, x̂bk, is corrected by a term formed of the gain matrix Kk and the innovation vector ykoHkx̂bk. The gain matrix Kk weights the information by taking into account the Gaussian observation and background model error covariance matrices, Rk and P˜kb, respectively. In the context of the LETKF, P˜kb is the sample covariance matrix obtained from the background ensemble perturbations, i.e. the deviations from the ensemble mean. The innovation vector holds the differences between observations yko used for assimilation in the considered time step and the ensemble mean x̂bk in observation space. Hk is the observation operator mapping model variables to observation space for comparison. The scaled analysis ensemble perturbations Zka are obtained from the background ensemble perturbations Zkb and the transformation matrix TkLETKF. The derivation of TkLETKF is comprehensively described in Hunt et al. (2007) and Livings et al. (2008). The analysis model ensemble after the update, xk,ia with i=1,,l denoting the member, is obtained by adding the new analysis perturbations, Zk,ia, to the analysis ensemble mean. The performed steps are summarised in Fig. 1. The analysis is performed independently at each grid point to increase the degree of freedom of the analysis otherwise limited by the low-rank approximation of the ensemble. After performing the analysis step, the data assimilation cycle is completed by running a forecast until the next update time step. The forecast is performed by applying the forecast model M on every model ensemble member separately, which allows for non-linear forecast processes. In the case of the presented study, the system to be described is areal precipitation information, the forecast model used to predict the areal precipitation evolution is the ensemble nowcasting scheme introduced in the following section.


Ensemble nowcasting

The nowcasting scheme used for the proof of concept study is a simple approach based on template-matching between time-lagged precipitation images. The original algorithm, the Automated Precipitation Extrapolator (APEX), was developed by van Horne (2003) and is adopted here with some modifications.

The algorithm considers two consecutive precipitation images and computes an estimate of the precipitation displacement from the first to the second image. A correlation is calculated between the second image and spatially shifted versions of the first image. The spatial shift is performed in every direction and with different amplitudes and the shift yielding the best correlation allows for estimating the displacement of the precipitation structure between both images. The resulting displacement vector field indicates the precipitation motion direction and speed, as demonstrated in Fig. 2. The computed displacement between both images is then used to extrapolate the precipitation structure to its position at future time steps.

Fig. 2.  

Example of input precipitation images for template-matching for the 3 July 2013. (a) 15:23:00 UTC and (b) 15:26:00 UTC and resulting displacement vector field, thinned out to every tenth pixel in each direction in this representation.

The complete displacement vector field (one vector per pixel in the precipitation field) is used to compute the precipitation forecast. Displacement vector lengths are scaled to the desired forecast time step and the future position of the precipitation field in the second image is predicted accordingly. In order to get a smooth, continuous precipitation field after the advection, a 3 × 3 window of pixels around each pixel is shifted to its new position. If more than one precipitation value is attributed to a pixel of the forecasted precipitation field, values are averaged, as illustrated in Fig. 3. The APEX algorithm computes a hierarchical correlation analysis considering first the whole domain of the images and then smaller subregions. Due to the computation of local displacement vectors at each pixel covered with precipitation, some divergence and rotation is allowed in the vector field. The original precipitation can therefore be distorted during the forecasting process, accounting for some internal variability in the structure of the precipitation field. Precipitation intensity can change compared to the initial field due to the averaging in the forecasting step, but for the sake of simplicity no cell growth and decay model is integrated in the nowcasting scheme. The prediction of the internal evolution of precipitation systems is complex. More elaborated probabilistic nowcasting methods which explicitly try to account for the uncertainty arising from growth and decay of precipitation cells exist (e.g. Bowler et al., 2006; Atencia and Zawadzki, 2014). In a further step subsequent to this study, it is possible to extend the presented nowcasting method to include this aspect, or even to replace it with a method of choice.

Fig. 3.  

Schematic representation of the forecasting step for three pixels, exemplary. 3 × 3 windows of pixels (initial position left) are shifted according to the local motion vectors calculated for the middle pixel. If these pixel windows overlap in the forecast (new position right), pixel values (precipitation intensity or motion velocity) are averaged to get the forecasted field. In this example, light grey pixels are averages over two values, dark grey over three values.

To create an ensemble forecast from the deterministic nowcasting scheme, we generate an ensemble of displacement vector fields from the deterministic field by randomly perturbing its x- and y-components. The perturbations are multiplicative and drawn from a Gaussian distribution with mean 1 and variance 0.4. The aim of the ensemble generation process is to produce a smooth perturbation of the displacement vector field. Therefore, the random noise from which the perturbations are built must be spatially correlated to avoid fuzzy and inconsistent displacement vectors. To create the pattern of the spatial correlation, areas with high wind speeds are assumed to behave similarly, i.e. perturbations applied to vectors with large displacement magnitude are strongly correlated. In this way, the structure of the precipitation field is reflected in the correlation of the applied random noise. The ensemble size for this study is 50 (results with smaller or larger ensemble size do not differ qualitatively from the results presented here). The detailed procedure as well as a sensitivity study for the chosen variance of the Gaussian noise is described in Merker (2017).



This study uses radar reflectivity data from a research radar network. Both the network and the data are thoroughly described in Lengfeld et al. (2014), together with all applied post-processing algorithms (including calibration and attenuation correction). The radar network comprises four X-band radars with a range of 20 km. Data are available as single radar product with 1° azimuthal, 60 m radial and 30 s temporal resolution, and in the form of a network composite product. The composite data have the same temporal resolution as the single radar data and a spatial resolution of 250 m × 250 m on a quasi-Cartesian grid, with a total of 297 × 213 pixels. The quasi-Cartesian grid is computed on rotated spherical coordinates with rotated pole at 36.0625° E, 170.415° S (refer to Lengfeld et al., 2014 for details on the composite generation). The location of the four radars, the domain covered by the radar composite and the position of additional micro rain radars (MRR) and gauges for reference purposes are indicated in Fig. 4.

Fig. 4.  

Map of the radar network. X-band radar locations and maximum range in black, MRR locations in orange. X-band radars and associated MRRs are installed at BKM, HWT, MOD and QNS sites. Three MRRs are installed between X-band radars, at sites WST, MST, and OST, from west to east respectively, together with reference rain gauges.

The Cartesian composite radar data are used to initialise the nowcasting method presented in Section 2.2. For simplicity and because the coordinate system is irrelevant to the presented study, all figures will be presented in rotated coordinates.

Observations considered for assimilation in this study are taken from single radar measurements (X-band radar MOD in Lengfeld et al., 2014, lower centre circle in Fig. 4). Observations are thinned to a 5000 m × 5000 m grid (20 × 20 grid boxes in the nowcasting domain) within the reach of X-band radar MOD, resulting in 52 observation locations (Fig. 5). Each reflectivity observation is selected from the radar’s polar grid by using the nearest neighbour method. These observations are illustrative for any additional information about the system and could be extended by e.g. data from a rain gauge network in a further development. In the setup of this study, observations are treated similarly to radar data, which is most commonly thinned or superobbed on a spaced, regular grid for direct radar data assimilation. This is done to reduce computational costs and representativity disparity and especially to avoid correlation of observation errors (e.g. Alpert and Kumar, 2007; Simonin et al., 2014; Bick et al., 2016; Waller et al., 2016). In addition, thinning of dense observations is needed to obtain an analysis ensemble not completely dominated by the amount of observed values, and to prevent collapsing of the ensemble spread. Ground observations are not used here to avoid issues induced by point-to-area differences and the empirical relation between precipitation intensity and radar reflectivity. Both would add additional uncertainty to the study and increase the complexity of results to a level that is beyond the scope of the intended proof of concept.

Fig. 5.  

Locations of thinned precipitation observations created for assimilation from data of the X-band radar at MOD site (Lengfeld et al., 2014) on a 5000 m × 5000 m grid (dots) and locations for verification on a 5000 m × 5000 m grid shifted by 2500 m north and east (crosses).

Observation errors on the thinned 5000 m × 5000 m grid are assumed to be uncorrelated and with constant error variance. This is a simplifying assumption, especially for radar data. The estimation of the full radar observation error statistics is a complex and ongoing topic of research (e.g. Xu et al., 2007; Berenguer and Zawadzki, 2008; Waller et al., 2016). Recent studies show the benefits of more sophisticated error models, which allow for less thinning or superobbing of the data and can significantly improve the quality of the analysis as shown for example in Jacques and Zawadzki (2014, 2015). However, the use of uncorrelated error statistics still is a state-of-the-art procedure (e.g. Chung et al., 2009; Simonin et al., 2014; Bick et al., 2016) and will also be applied here. The error variance of the radar data is estimated by comparing data from the network presented in Lengfeld et al. (2014) to reference measurements from MRR in order to get an error variance specific to the radar data used. The error variance is found to be 11.32 dB2.

We focus on reflectivity data to avoid additional uncertainty coming from the empirical conversion from reflectivity to rain rate. The logarithmic reflectivity unit dBZ is used throughout the study. For reference, a difference of 1 dB translates approximately to a factor 1.25, and 2 dB to a factor 2 difference in non-logarithmic reflectivity Z. For rain rate in mm h–1, these differences correspond to roughly a factor of 1.15 and 1.54, respectively. Following the example of Bick et al. (2016), the ‘no rain’ value is set to 5 dBZ, approximately 0.07 mm h–1, and all reflectivity data are thresholded accordingly throughout the study.



To assess the uncertainty estimate obtained by the method introduced in Section 2, we run coupled nowcasting-assimilation experiments using reflectivity data presented in Section 3. The experiment setup and the detailed evaluation of the method in the following Section 5 is first presented based on a precipitation event monitored by the radar network on 3 July 2013 between 15:30:00 UTC and 16:00:00 UTC (case 1), in which a non-stratiform cell passes through the radar network area. Thereafter, results from three additional precipitation events are included in the evaluation to cover relevant regimes. Stratiform rain dominates on the 11 August 2013 (18:30:00 UTC–19:00:00 UTC, case 2), convective precipitation characterise the event on the 13 August 2013 (09:30:00 UTC–10:00:00 UTC, case 3) and the last case on the 17 August 2013 (03:30:00 UTC–04:00:00 UTC, case 4) features a frontal passage. An impression of these precipitation events is given in Fig. 6.

Fig. 6.  

X-band radar network composite reflectivity data (Lengfeld et al., 2014) for the four cases presented in the study: Case 1 on the 3 July 2013 15:23:00 UTC (a) and 15:26:00 UTC (b), case 2 on the 11 August 2013 18:30:00 UTC (c), case 3 on the 13 August 2013 09:30:00 UTC (d) and case 4 on the 17 August 2013 03:30:00 UTC (e).

As presented in Section 3, the ensemble forecast for the data assimilation experiment is initialised with radar composite data. Considering the first case of the study, forecast start is 15:26:00 UTC. The displacement vector field for the forecast is computed between composite radar images at 15:23:00 UTC (Fig. 6a) and 15:26:00 UTC (Fig. 6b). The determined global displacement is three pixels towards the east and four pixels towards the north, i.e. u= 4.2 m s–1 and v= 5.6 m s–1 respectively. The speed of the precipitation cell is thus estimated to be 7 m s–1. The forecast is performed for 50 ensemble members until 16:00:00 UTC in time steps of 2 min.

Observations used for assimilation are selected from single radar data and thinned to a 5000 m × 5000 m grid, as indicated in Section 3. We assimilate observations at the 52 locations resulting from the radar data thinning every 4 min, i.e. every two time steps, from 15:30:00 UTC to 16:00:00 UTC. This assimilation frequency is similar to the temporal resolution of 5 min of the C-band radar data made available operationally by the German Meteorological Service (DWD) and the assimilation frequency used in Bick et al. (2016). The data assimilation cycle ends with a new ensemble nowcast started after the analysis step. For this new nowcast, we compute displacement vector fields from the precipitation fields at the current and the previous analysis time for each ensemble member, without applying additional perturbations.

The LETKF further requires the specification of an observation influence radius. During the analysis step, model grid points are affected by all observations within the influence radius. Other observations are ignored. We use an influence radius of 4000 m, which ensures a complete coverage of grid points within X-band radar MOD reach. Therefore, every grid point in this area is affected by the update during data assimilation which favours a smooth precipitation analysis. An analysis of the sensitivity of the system with respect to different observation influence radii can be found in Merker (2017).

To illustrate the system’s behaviour, the performed data assimilation cycle is shown in Fig. 7 for one arbitrary observation point in the model domain, indicated by a cross in Fig. 8. The uncertainty of the ensemble mean is indicated by the ensemble spread, which is defined here as the ensemble standard deviation. During the data assimilation cycle, the analysis ensemble spread depends on the forecast uncertainty described by the background error covariance matrix, and on the observation uncertainty described by the observation error covariance matrix. The forecast uncertainty results from the errors associated with the extrapolation scheme used for the precipitation forecast. The ensemble spread is zero at the beginning of the forecast because all ensemble members are initialised with the same reflectivity field (Section 2.2). After one time step, ensemble members start to diverge and the standard deviation of reflectivity at the considered location increases. At time steps where observations are available, the assimilation pulls the ensemble mean towards the observation and reduces the ensemble spread at grid points influenced by the observations. The flow-dependent background error covariance matrix P˜b and the observation error covariance matrix R (Eq. 1) determine the effect of this new information on the ensemble at affected grid points. The observation error covariance used for assimilation is constant in this study and was derived from statistical data analysis (Section 3). After the assimilation step, the ensemble evolves according to the forecast model until the next assimilation is performed. The model ensemble and its spread are subject to the forecast and therefore the uncertainty information contained in the ensemble spread develops a flow dependency. In this way, the presented method allows for a physically consistent propagation of uncertainty in time.

Fig. 7.  

Reflectivity ensemble mean evolution (dark blue line) and uncertainty range (light blue envelope, ± one ensemble standard deviation) throughout the data assimilation cycle at observation grid point 0.041°E, – 0.053°N (indicated by a cross in Fig. 8) and observations (orange whiskers, ± observation error standard deviation).

Fig. 8.  

Spatial distribution of the reflectivity ensemble mean (with contour highlighting the region in which 80% of the ensemble members show precipitation above 5 dBZ, top row) and spread (bottom row), i.e. ensemble standard deviation, for the analysis at (a, c) 15:34:00 UTC and (b, d) 15:46:00 UTC on the 03 July 2013. Circles indicate observation influence radii, the cross indicate the location of the observation location used as an example in Fig. 7.



The spatial structure and temporal evolution of the uncertainty information is the main asset of the method presented here. The structure and evolution of the uncertainty is described by the ensemble spread at analysis time, which is shown together with the ensemble mean as an example for two time steps in Fig. 8. The potential of this uncertainty estimate is demonstrated in the following. The new uncertainty estimate is shown to be more accurate than a constant benchmark value, which assumes a constant uncertainty to be representative for the precipitation system. Prior to that, a short verification of first-order statistics is performed based on precipitation case 1.

Verification is performed at locations selected on a 5000 m × 5000 m grid, similar to the assimilation locations. The verification grid is shifted by 2500 m in each direction with respect to the assimilation grid in order to select independent observation locations for assimilation and verification (Fig. 5). The verification grid consists of 47 grid points within X-band radar MOD range. Reflectivity forecasts and observations are compared at verification grid points at each of the eight assimilation time steps, after the update step. Thus, the comparison comprises 376 data points in total.

The overall root mean square error (RMSE) of the forecast ensemble mean for case 1 amounts to 4.49 dB and the bias, defined here as the mean error of the ensemble mean compared to the observations, to 1.54 dB. These forecast statistics can be improved slightly using the probabilistic ensemble information and applying a 80% precipitation probability threshold to the forecast data, i.e. setting the ensemble mean for grid points with precipitation probability lower than 80% to the ‘no rain’ value of 5 dBZ. Then, RMSE and bias decrease to 4.44 dB and 1.14 dB, respectively, showing a slight improvement over using model ensemble mean without integrating probabilistic information. The bias indicates an overestimation of precipitation by the forecast model. This is an expected feature considering the simple extrapolation scheme used for the study. Cells can be spread over different regions among the ensemble members. Especially cell maxima then lead to an overestimation of the precipitation by the ensemble mean compared to colocated observations. The bias is small compared to the random uncertainty reflected by the RMSE and will have little impact on the results of the proceeding data assimilation experiments. The RMSE of 4.44 dB is approximately 1 dB larger than the estimated accuracy of the original radar data of 3.36 dB (Section 3). Considering the simplicity of the nowcasting scheme used here and the forecast range of 30 minutes, results of the forecast ensemble mean are found to be reasonable.

To characterise the quality of the areal uncertainty estimate of the combined precipitation, we analyse the statistical relation between forecasted uncertainty, i.e. the ensemble standard deviation, and the actual error, i.e. the deviation between forecasted and observed reflectivity. In a so-called spread-skill diagram, the perfect relation between forecast uncertainty and forecast error is a one-to-one relation. The uncertainty forecast perfectly describes the uncertainty of the system if, statistically, the actual forecast error equals its forecasted uncertainty (Leutbecher and Palmer, 2008). For all verification grid points, absolute forecast error and corresponding ensemble spread are compared in a spread-skill diagram (Fig. 9). Ensemble spread is binned into classes of 0.5 dB width from 0 to 8 dB. Statistics of the absolute forecast error distribution within these classes are described by the median and the first and third quartiles.

Fig. 9.  

Comparison of absolute precipitation product (model ensemble mean) error and ensemble spread (model ensemble standard deviation) at the available 47 verification grid points and eight analysis time steps for (a) case 1, (b) case 2, (c) case 3 and (d) case 4. Ensemble spread values are divided into bins of 0.5 dB width, boxes indicate the median (blue line) and the first and third quartiles. Data distribution is shown in frequency histograms, the solid grey lines show the cumulative distribution function.

Most data points show ensemble spread and absolute product error below 0.5 dB. This category is dominated by grid points without precipitation, both in the forecast and the observations. Overall, there is a clear correlation between ensemble spread and absolute product error, especially for cases 1 and 3, hinting at a good uncertainty representation by the ensemble spread. In the range between 0 and 3 dB ensemble spread, the forecast model ensemble is overdispersive because ensemble spread is substantially larger than the median of the absolute product error distributions (median below the one-to-one line). This is due to non-zero model forecast values at grid points outside the observed (true) location of the precipitation cell. Most of the ensemble members correctly predict no precipitation for those grid points, but some have higher reflectivity values because of the differently predicted displacements of the cell. The ensemble standard deviation is more strongly impacted by these few high reflectivity values than the ensemble mean. Therefore, the ensemble standard deviation is larger than the deviation between ensemble mean and observation.

Two scores are defined to quantify the potential of the ensemble spread to correctly represent the uncertainty of the precipitation product in time and space. The aim is to prove that the spatial and temporal structure of the uncertainty estimate is not random, but actually provide valuable additional information. The first score indicates the amount of data points for which the absolute product error falls within the uncertainty range predicted by the ensemble spread—the percentage of ‘hits’. The ensemble spread is interpreted as the tolerated margin of error. The score

((2) )
represents the reliability of the product uncertainty estimate, where ϵi and σi are the absolute forecast error and the uncertainty estimate for each available verification data point i, respectively. The second score measures the deviation from the perfect spread-skill relation, i.e. from a one-to-one relationship. It is calculated following the definition of the root mean square deviation between ensemble spread and absolute model error at the verification grid points,
((3) )

To assess the areal uncertainty estimate for the combined precipitation product, above scores are computed with the spatially and temporally variable product uncertainty described by the ensemble spread after the update step at all available verification grid points. The reliability RELvar is computed with σi=σvar,i, where σvar,i is the product ensemble spread, and ϵi, the corresponding absolute difference between forecast ensemble mean and observations from X-band radar MOD at the same grid points and time. For the spread-skill deviation DEVvar, it is taken into account that the relation between ensemble spread and forecast skill is a statistical one. Therefore, the ensemble spread σvar,i is binned into j classes of 0.5 dB width with class centre σbin,j, as shown in Fig. 9. DEVvar is then computed using σbin,j and the median of the model error distributions in each class, ϵbin,j. This also reduces the effect of outliers on the DEV score.

The benchmark against which the spatially and temporally variable uncertainty estimate gained by the flow-dependent ensemble spread is evaluated is a constant uncertainty estimate valid for all grid points at all-time steps. Constant uncertainty information still is most common for precipitation data from radar or other sources (e.g. Tong and Xue, 2005; Bick et al., 2016). To facilitate the comparison, the constant uncertainty should conserve the total amount of uncertainty given by the variable uncertainty estimate among the considered data points. In this way, the focus of the analysis is on the ability of the presented method to correctly distribute the available amount of uncertainty information in space and time. Therefore, the chosen reference uncertainty estimate is a mean ensemble spread of the system, σ¯, calculated as an average over the ensemble standard deviation at verification grid points, including all analysis time steps. In order to compare the spatially and temporally variable uncertainty estimate of the combined precipitation product σvar,i to the benchmark, RELconst and DEVconst are computed from Eqs. 2 and 3 by replacing σi by the constant value for σ¯.

Both scores are summarised in Table 1 for the four precipitation events introduced in Section 4: case 1 with σ¯= 2.71 dB (convective), case 2 with σ¯= 3.10 dB (stratiform), case 3 with σ¯= 2.37 dB (convective) and case 4 with σ¯= 3.28 dB (frontal). The reliability REL of the uncertainty forecast improves for all four cases when using the areal uncertainty estimate provided by the data assimilation cycling method compared to using constant benchmark uncertainty estimates. RELvar improves by 36.80 (19.64, 8.73, 8.56) percentage points compared to RELconst for case 1 (case 2, case 3, case 4). This means that estimate σvar,i provides a more accurate distribution of uncertainty among the system than the reference σ¯. The resulting information on the reliability of the precipitation product is therefore improved in all analysed cases, i.e. the tolerated margin of error represented by σi is satisfied more often.

The results of the DEV score do not show such a clear signal. DEVvar is reduced by 1.92 dB (3.37 dB), i.e. by more than a factor of two (one and a half), compared to DEVconst for case 1 (case 3). For case 2 (case 4), the DEV score does not perform better when using σbin,j instead of σ¯. This is due to areas where the ensemble is overdispersive, i.e. where large uncertainty is predicted (ensemble spread) but the error of the precipitation product is moderate (exemplarily in the right half of Fig. 9b for case 2). These values occur in a minority of cases but have a large impact on the DEV score after the binning. The large spread among the ensemble members in those few cases may be due to the fact that both events presented in cases 2 and 4 cover a large part of the considered area and have moderate spatial gradients in the measured precipitation fields, which complicates the derivation of displacement vectors for the nowcast. An improvement of the results is therefore expected by using a more physically-based probabilistic nowcasting scheme that yields more adequate ensemble members.


Summary and conclusions

This article presents a method to estimate spatially and temporally variable uncertainty for a combined areal precipitation product. The method makes use of data assimilation to merge precipitation measurements from different sources and an ensemble nowcasting method to evolve the uncertainty estimate over time. The potential of the areal uncertainty estimate provided by the method is demonstrated in a proof of concept study.

Requirements for this improved uncertainty estimate are an accurate representation of the actual error of the initial product, an adjustment to additional observations merged into the product through data assimilation, and flow dependency. In this study, a research network of four X-band weather radar provides reflectivity data for four precipitation cases. Reflectivity can be converted to rain rate using empirical relations that are not part of the analysis.

An ensemble data assimilation method, the LETKF, is used to statistically merge precipitation observations and is coupled to an extrapolation-based nowcasting scheme. The implemented nowcasting scheme computes the cross-correlation between subsequent radar composite images and extrapolates the evolution of the precipitation field using the deduced displacement. The nowcast is started from composite data of the research radar network and a probabilistic forecast is generated using an ensemble technique. The ensemble is generated by stochastic perturbation of the computed precipitation displacement vectors.

Four data assimilation experiments are performed to test the presented method with emphasis on its potential to provide an improved spatial and temporal uncertainty estimation. We use an ensemble precipitation nowcast and continuously merge observations into the forecast using the LETKF, generating a combined precipitation product. The uncertainty of the precipitation product is estimated by the ensemble spread of the nowcast after each update step of the data assimilation cycle. Two scores are introduced for the assessment of the method. Both are based on the definition of a perfect uncertainty estimate, for which the actual observed error statistically corresponds to the predicted uncertainty. The first score describes the reliability, REL, of the uncertainty estimate. It indicates the percentage of cases in which the actual error of the precipitation product lies within the estimated uncertainty range. The second score measures the deviation of the uncertainty estimate from a perfect spread-skill relation, DEV, as a root mean square deviation between actual error and uncertainty estimate.

The scores are computed at verification grid points selected on a regular grid and for all available analysis time steps. The potential of the obtained areal uncertainty estimate σvar,i is assessed by comparison to a benchmark which must be outperformed. The benchmark for this study is defined as the mean spread of the system. In this way, the ability of the method to correctly distribute the available mean uncertainty of the system over time and space can be studied.

The reliability of the uncertainty forecast, REL, is improved by the presented method. The number of hits, i.e. how often the error of the precipitation product is within the tolerated margin of error dictated by the predicted uncertainty range, is increases in all analysed cases. The spread-skill deviation, DEV, is improved in two out of four cases. This score yields unsatisfying results when the spread of the forecast ensemble is large compared to the error of the obtained precipitation product. The large ensemble spread is most likely caused by difficulties in predicting a good sample of displacement fields and large gradients in the precipitation at the edges of the cells in some cases. We expect the results to improve for those cases with a more elaborated probabilistic nowcasting scheme. Despite mixed results for the DEV score, the clear improvement of the uncertainty forecast visible throughout the analysed precipitation events already show that the presented areal uncertainty estimate allows for a more accurate spatial and temporal distribution of the uncertainty information. Since we assume a large impact of the nowcasting scheme on the results for the DEV score, we therefore see potential for further improvement by implementing a more sophisticated scheme.

The evaluation of both considered scores demonstrates that the provided areal uncertainty estimate outperforms constant benchmark uncertainty values, but that the method is sensitive to the quality of the probabilistic nowcasting. In subsequent work, the next logical step would therefore be the use of a more sophisticated nowcasting tool. Nevertheless, the proof of concept shows the potential of the developed method and establishes the groundwork for further studies. The possible applications of the method are numerous in hydrology, nowcasting or data assimilation. It provides the combination of areal and point measurements for comprehensive precipitation products and additionally yields probabilistic information allowing for the estimation of the product reliability and ensemble generation.