Automated convective and stratiform precipitation estimation in a small mountainous catchment using X-band radar data in Central Spain

High spatio-temporal resolution quantitative precipitation estimation (QPE) is an increasing requirement for many applications (e.g. weather forecasting, water resources assessment, flood risk management), and is especially important for hydrological modeling of floods during extreme events, allowing the uncertainty reduction of the inputs in the precipitation-runoff models and runoff forecasts (Zoccatelli et al. 2010; Gourley et al. 2011; Krajewski et al. 2011; Zhang et al. 2011; Lee et al. 2014). For most applications, the real-time estimation requirement adds a significant challenge to the interpolation procedure, as this should then be both automatic and efficient, providing robust interpolation of gauged data.

Traditionally, rain has been measured at ground level using rain gauges which record precipitation intensity (R); however, with them spatial representation is often low due to their sparse density (Villarini et al. 2008). Moreover, geostatistical interpolation methods (Creutin & Obled 1982; Grimes et al. 1999; Goovaerts 2000; Kyriakidis et al. 2001; Lloyd 2005; Carrera-Hernandez & Gaskin 2007; Chappell et al. 2013; Berndt et al. 2014; Jewell & Gaussiat 2015; Nanding et al. 2015; Rabiei & Haberlandt 2015) such as univariate (e.g. ordinary kriging) and multivariate (e.g. kriging with external drift (KED) or co-kriging) require a minimum of 50–100 data points to ensure an appropriate geostatistical analysis (Journel & Huijbregts 1978; Webster & Oliver 1993). This implies an unrealistic rain gauge density for very small catchments, impossible to attain even in developed countries.

Developed from the 1940s onwards, weather radars have been increasingly used for QPE in hydrological modeling due to their high spatial resolution, real-time availability (Cole & Moore 2008; Ehret et al. 2008) and also a higher coverage than ground observations (e.g. over the sea). In contrast to rain gauges, weather radars do not provide a direct measure of rain, but record the reflectivity (Z) of particles at a certain height above the ground (Rogers & Yau 1989). A weather radar has the capacity to sweep the entire study area at a high time resolution a few minutes for S- and C-band radars (Rogers & Yau 1989; Fabry 2015), or even less than a minute for some X-band radars. Still the data it yields are subject to various errors that require correction (Hong & Gourley 2015). Marshall & Palmer (1948) related drop-size distribution to R, giving a power law of the reflectivity-to-rainfall relationship (or Z-R equation).

Z-R coefficients may vary depending on the origin of rainfall, the region, temperature, radar features, etc. For a further discussion on the power law relationship between Z and R, and the corresponding values of a and b coefficients for different kind of precipitations, Chapter 9 of Fabry (2015) can be consulted. All these variations prevent the direct use of radar rainfall fields obtained from classical Z-R equations, even with adjusted coefficients (Borga 2002). Therefore, sophisticated merging methods combining rainfall, radar and other available short time step data are required. Among the high number of radar-gauge merging techniques, geostatistical estimators in particular have been subject to intense development (Krajewski 1987; Creutin et al. 1988; Seo et al. 1990; Seo 1998; Sinclair & Pegram 2005; Schuurmans et al. 2007; Ehret et al. 2008; Velasco-Forero et al. 2009; Heistermann & Kneis 2011; Verworn & Haberlandt 2011, among others). Even within the area of geostatistical methods, there is a wide range of algorithms with varying degrees of complexity and success (e.g. KED or co-kriging). The method named ordinary kriging of radar errors (OKRE) by Erdin (2013) is similar to the combination approach suggested by Sinclair & Pegram (2005), known as conditional merging (MERG). Several authors have found that these methods outperform other geostatistical algorithms such as KED (Erdin 2013; Legorgeu 2013; Berndt et al. 2014), or perform very slightly worse (Goudenhoofdt & Delobbe 2009).

Small-scale rainfall variability should be taken into account to improve rainfall data resolution. Both QPE uncertainty and variability will affect the confidence of hydrological model predictions. All the references mentioned above used a C-band radar resolution of ∼1 km, rendering it difficult to capture precipitation microscale variability, particularly strong in convective events (Sideris et al. 2014). Three of these studies were devoted to very small mountainous catchments, with quick rainfall-runoff response. Ehret et al. (2008) merged a C-band radar and gauge data in a 75 km² catchment with a resolution of 500 m, preserving both the relatively reliable rain gauge data and the accuracy concerning spatial variability of the radar image. Meanwhile, Germann et al. (2009) conducted a hydrological experiment analyzing C-band radar data with a resolution of 1 km in a small, steep catchment of 44 km². They found that radar estimations yielded overestimations compared with rain gauge measurements for convective situations with strong spatial rainfall variability. Heistermann & Kneis (2011) benchmarked six different QPE methods in two small mountainous basins of 15 and 49 km². They used one C-band radar with spatial resolution of 1 km and 16 rain gauges situated close by the study area, but only one of them is located inside one runoff gauge sub-basin. They found that the merging method used by Ehret et al. (2008) performs the best according to cross validation and hydrological verification based on runoff simulation.

X-band radar measurements can have a resolution of 100 to some meters (Rico-Ramirez et al. 2007), registering small-scale precipitation heterogeneities (Anagnostou et al. 2010). Thus, X-band radar is ideally suited for local rainfall estimation (Delrieu et al. 1997; Gires et al. 2012), as illustrated in the PhD dissertation of Legorgeu (2013), who applied geostatistical techniques in an area of 1,300 km². Furthermore, for stream flow simulation purposes Shrestha et al. (2013) found X-band weather radar based estimations of precipitation to be more accurate than estimates using C-band radar.

The objective of the present study is to obtain a reliable spatial and temporal distribution of precipitation intensity for a small mountainous catchment. Geostatistical interpolation is used to merge measurements in rain gauges and X-band radar data. Analysis of the results is performed distinguishing between stratiform and convective rainfall.

A maximum time scale of 1 hour is required to capture time variability of fast rainfall events (Sideris et al. 2014). Smaller accumulation time periods present degradation due to the fact that discrepancies between radar and rain gauge measurements are larger for such short aggregation periods (Sideris et al. 2014). The authors found that, in the context of real-time no-human-intervention applications, hourly input accumulations give a satisfactory balance between very high temporal resolution (non-robust), aggregations and very robust but low temporal resolution. Thus, as used in other geostatistical studies of QPE (Harberlandt 2007; Velasco-Forero et al. 2009; Verworn & Haberlandt 2011), 5 min rainfall values were aggregated to hourly values. Accumulated hourly radar images were computed assuming that precipitation fields move at a constant velocity between images and vary linearly in intensity with time between each interval (Bellon et al. 1991).

A map of false echoes was available, obtained on different clear days. To correct these false echoes, this map was automatically subtracted from each reflectivity image. It was later proven that for clear days without precipitation, the reflectivity image is almost white, with tiny echoes remaining which are usually present despite correction (Meischner 2004). Quantifying attenuation was difficult since this is not a polarimetric radar and there is no receiving antenna in the basin, but a comparison of different radar images for intense rain events confirmed that attenuation in the area was not significant. Furthermore, the radar was so close to the basin that attenuation would not be high.

According to the power law established by Marshall & Palmer (1948), the relation between dBZ (10 log₁₀Z) and log₁₀R must be linear. This renders it possible to obtain linear correlation maps (centered on a rain gauge) between log₁₀R and dBZ at nearby pixels. In the absence of anomalous propagation, the most correlated point will be the center of the map, making it possible to compare rainfall data with reflectivity at rain gauge coordinates. When there are strong temperature inversions, the radar beam suffers super refraction and it bends toward the ground. In these situations, the most correlated point would show a radial displacement with respect to the gauge coordinates, indicating a problem of anomalous propagation (Bean & Dutton 1966; Pamment & Conway 1998; Fabry 2015).

It is important to note that the same approach employed by Goudenhoofdt & Delobbe (2009) and Velasco-Forero et al. (2009) was followed, and rainfall data were used in their original form rather than transforming them into Gaussian space.

Another important and challenging issue affecting radar QPEs is related to the precipitation origin: either convective or stratiform events. Systematic discrimination between convective and stratiform events can be performed with different techniques (Anagnostou & Kummerow 1997; Awaka et al. 1997; Hong et al. 1999; Llasat 2001; Verworn & Haberlandt 2011; Nanding et al. 2015). Monjo (2016), using a dimensionless index to separate events by origin, defined according to the relative temporal distribution of maximum intensities. The main advantages of this method are that it does not require thresholds, so it can be applied for each rain gauge.

Definition of the n index is based on the fact that stratiform rainfall (advective or frontal) normally has a constant intensity during all events, whereas convective rainfall is characterized by variable intensities within the same event. Therefore, n values close to 0 are associated with constant intensities (stratiform event), while values close to 1 indicate variable intensities of convective events. In the present study, hourly n indexes were computed for each rain gauge and the median of the n distribution was used as the threshold for considering rainfall as predominantly stratiform. Comparison of convective available potential energy values and n indexes classification validated the stratiform–convective separation in the present work.

A theoretical model must be fitted in order to deduce the variogram values for any possible distance required by interpolation algorithms (Goovaerts 2000). The variogram model was calculated for each time step, as each calculation time was treated independently from the previous ones. Hence, an automatic approach was applied with the geostatistical R package, gstat (Pebesma 2004), and the model type that best fit experimental variograms was established. The sill, nugget and range were determined by the fitting procedure. The sill is the upper boundary of the variogram function, equal to the variance of the underlying studied population. The nugget represents any discontinuity at zero separation and the range is the distance over which the covariance is valid (i.e. the extent of the spatial structure in the data).

This was then followed by subtracting the interpolated radar error from the radar rainfall values estimated from Z-R relations to obtain the actual precipitation estimate. Differences in the variogram structure between the radar and radar errors explain the miscalculations of the present method. Ordinary kriging was built using the geostatistical R package gstat (Pebesma 2004).

To validate the interpolation results, the ‘leave-one-out’ or cross-validation method was applied. A successive estimation of all gauge locations was performed by using all other stations and always omitting the sample value at the location in question. Three performance measures or scores were used to compare estimations and observations, as described below. The corresponding analysis is presented in the Results and discussion section.

This was employed to distinguish between dry and wet areas. The term discriminant that is sometimes used refers to this score's characteristic of measuring discrimination between yes and no cases. The score has a range of −1 to +1, where 0 represents no skill or a random estimate and 1 represents a perfect estimate. HK is also the difference between the hit rate and false rain rate. Therefore, it makes it possible to examine whether forecasting an event, particularly a rare event, more often leads to a large increase in false rains or not.

This is similar to a root-mean-squared error, but is adapted to the skewed nature of precipitation. This behavior is reflected in the fact that intense events are less predictable than light precipitations or dry records. In addition, it is normalized by the mean root-squared deviation in order to mitigate dependency on the observed precipitation deviation. RMRTE = 0 indicates a perfect estimation, whereas RMRTE = 1 is equivalent to approximate precipitation with the mean rainfall value.

The study area, the Venero Claro Basin (Ávila, Spain), has suffered severe flash floods in the past, and these thus need to be predicted in order to improve awareness measures and mitigate future damage. Although the catchment (Figure 2) may appear to have a high density of instruments, six rain gauges are still not enough to directly derive spatial variability in rainfall. Venero Claro is the drainage basin of the La Cabrera Stream, which measures 5.5 km long, and it covers an area of 15.5 km² with a maximum altitude difference of 1,224 m between the highest point (1,959 m a.s.l.) and to the point where the stream flows into the Alberche River (735 m a.s.l.). The average slope of the stream is 21.6%. Geologically, the area is located in the western part of the Spanish Central System. The bedrock consists of biotite granites and granodiorites from intrusions during the Variscan orogeny. Superficial Quaternary formations are made up of gravels, sands and silts. These cover the slopes, the valley bottom and endorheic depressions. Its morphology and composition is common to other small catchments in the Spanish Central System.

The climate of the study area is Continental Mediterranean, which is typical of inland Spanish mountain ranges. The climate of these mountains is determined by the frequent arrival of Atlantic depressions, mainly from the SW, during fall, winter and spring, and by the predominant Azores anticyclone which causes very dry summers (accounting for only 10% of the annual precipitation). The mean annual precipitation is 554 mm in the lower areas and up to 2,000 mm at the highest point (Palacios et al. 2011). Six recording stations provided continuous rainfall data between 2006 and 2014 over 5 minute time intervals. Rain gauge values are considered as estimates of true precipitation rate at the station location. In this study area there is a sparse rain gauge density when compared to the spatial variability of precipitation (Ruiz-Villanueva et al. 2011, 2013). During the rainfall data period, 6,006 hours had at least one rain gauge measuring. Only five rain gauges were operating during rainfall events with radar data (317 hours with rain measured by at least one rain gauge and with reflectivity images).

An X-band weather radar was installed in 2012, located 2,330 m west of the mouth of the La Cabrera Stream (Figure 2). This is an appropriate location as it is close enough to avoid strong attenuations and the radar is able to cover the entire vertical profile of the basin, with a fixed elevation angle of 12°. The technical specifications of the radar are given in Table 2. Recording spatial resolution was defined as 80 × 80 m or 120 × 120 m, while time resolution was 5 min. A region of 50 × 50 km centered on the radar was defined as the study area.

Regarding radar data analysis, the existence of anomalous propagation was investigated by looking for displacement patterns in the radar data. Generally, the best correlated pixel was the one corresponding to the gauge. Hence, no anomalous propagation correction was required. Due to the low number of rain gauges, it was not possible to apply the MERG or KED method as no variogram model could be inferred for rainfall. In contrast to Erdin (2013), rainfall data were not transformed in the present application of OKRE.

The lack of knowledge about spatial rainfall variability is one of the major sources of uncertainty in the stream flow model-derived information (Willem & Berlamont 2002). Hence, even if there are not as many rain gauges as in the applications of OKRE performed by Erdin (2013) with 75 rain gauges or by Nanding et al. (2015) with up to 161 rain gauges, the study of the flood predictions in the area needs to improve as much as possible the spatial heterogeneities of QPE. In the present work it is possible to avoid the reduction of the significance of the results due to the small number of gauges by taking the radar correction variability directly from the radar data. Assigning radar spatial variability to rainfall spatial variability will end up in some uncertainties on the results. However, the reduction on the uncertainty linked to the experimental variogram justifies this assumption. It may be here recalled that the experimental variogram is statistically robust when it is computed from a minimum of 30 data pairs (Journel & Huijbregts 1978). Similarly, Germann & Joss (2001) conducted precipitation variogram estimation using radar data and reported that high resolution radar images provided good information about the spatial continuity of precipitation. The same approach was taken by Berndt et al. (2014). The novelty of the methodology proposed here is the application of OKRE using the radar variogram model to krige the radar errors.

Theoretically, variogram estimation must be carried out separately for each time step. If this is done manually, the procedure is very time-consuming and subjective. A solution proposed by Schuurmans et al. (2007) is to use pooled variogram models based on all the studied events. They recommended this practice when there is a lack of data or when an automatic prediction procedure is implemented without variogram estimation. Berndt et al. (2014) and Rabiei & Haberlandt (2015) used seasonal variogram models, averaging variograms for all time steps within each of the two split seasons. Other authors such as Velasco-Forero et al. (2009) have used non-parametric automatic evaluation of the variogram model for each time step. In the present study, fitting a parametric variogram for each hour of estimation was considered of great importance to capture temporal variation in rainfall. As with Schuurmans et al. (2007) and Erdin (2013), automatic experimental calculation and modeling of the variogram was successfully implemented. Previous research by Ehret (2003); Verworn & Haberlandt (2011) and Erdin (2013) showed that the choice of variogram model had little impact on rainfall estimation performance. In this work, after different trials, the variogram model selected as the best option was the Gaussian model.

Bias values differ considerably between rain gauge locations. The only rain gauge for which OKRE gave a slight underestimation was the one located at the top of the basin (Peña Parda in Figure 6). This can be explained by the fact that when it was raining at this point, it also rained in the rest of the basin, but probably with less intensity. In these cases, interpolation was always performed taking into account smaller rainfall values. The elevation of the radar, 12°, also affected the reflectivity measurements, being less accurate at higher elevations.

With regard to the RMRTE, an improvement was achieved with the geostatistical method in comparison with the Z-R relation. The best estimations were obtained at the Dehesa location (closest to the radar). This good performance of the method with X-band radar data was also observed with C-band radar information by Ehret et al. (2008) and Nanding et al. (2015) when using the similar MERG method. Ehret et al. (2008) reported that MERG gave better results than when using only radar data with a static or updated Z-R relation, when using radar but determining the Z-R relation with a disdrometer, or when using KED.

Valuable information was extracted from the stratiform and convective ensembles of variogram models. The Gaussian model parameters were analyzed to determine characteristics of spatial continuity that differentiated stratiform and convective events, and quartiles of parameters in both ensembles were computed (Table 3). Sill median values were not very different for the two kinds of event (0.13 for convective and 0.40 for stratiform). The higher sill values for stratiform events (Q2 0.4 and Q3 3.73) indicated higher hourly precipitations than convective ones (Q2 0.13 and Q3 0.39). Median range values indicated that stratiform events in the area had a spatial correlation up to ∼5 km, whereas this was up to ∼4 km for convective events. This higher correlation length is explained by the larger dimensions of the stratiform rain processes, corroborating the consistency of the methodology. Similar characteristics can be seen in the experimental variograms shown in Jewell & Gaussiat (2015) for convective and stratiform events.

Figure 7 shows the three error scores obtained for all the points and times, but distinguishing between stratiform and convective origin. For the estimation using Z-R relations and no geostatistical merging, the bias error measurement indicated that it was almost unbiased for stratiform rainfalls, but highly underestimated for convective ones. This can also be seen in Figure 4, where points of higher intensity are far from the adjustment. This underestimation of spatially distributed precipitation estimates is confirmed by other authors (Van de Beek et al. 2010; Thorndahl et al. 2014; Nanding et al. 2015; Rabiei & Haberlandt 2015), and is particularly dangerous when predicting flash floods during convective events. When using Z-R adjustments, the probability of successfully distinguishing wet and dry periods (HK) is much higher in stratiform events compared with convective ones. For both rain origins, HK was much improved in the OKRE results. Error values computed with the RMRTE for estimations applying the Z-R estimations were very high, reaching the value of 1, and hence estimated values were comparable with the average precipitation of all rain gauges. These RMRTE's were also much improved by means of OKRE.

In general, the geostatistical method described here yielded better predictions of stratiform events than of convective ones. In comparison with stratiform rainfall, convective events presented higher overestimations, worse discrimination between dry and wet periods and a larger RMRTE. This has been also found by Velasco-Forero et al. (2009), Erdin (2013), Legorgeu (2013), Jewell & Gaussiat (2015), among many others. It remains to be determined how significant this improvement is in terms of rainfall product resolution and accuracy in flood simulations for this basin. As was also observed by Van de Beek et al. (2010); Verworn & Haberlandt (2011) and Jewell & Gaussiat (2015), the results obtained highlight the importance of meteorological conditions on the quality of the merged product.

The main innovation of the present work is the application of a merging radar-rain method to a small catchment (15 km²) with very small concentration times, where only a dense rain gauge network (e.g. 10 times higher than the World Meteorological Organization recommendations for small mountainous islands with irregular precipitation) combined with one X-band radar with a high spatial resolution, allows an adequate characterization of precipitation. Cross-validation results throughout the basin and the data proved a good performance of the interpolation approach. In comparison with the classical Z-R relationships, obtained here from a linear fit, the method described here showed a greater capacity for distinguishing between wet and dry areas and only presented minor errors in the estimate with respect to the observed values.

Another novel aspect of this work was the use of the n index to improve the Z-R adjustment for convective and stratiform precipitation. The main advantages of this method are that it does not require thresholds, so it can be applied for each rain gauge. Moreover, this index allows automatic separate variogram modeling of stratiform and convective events, which opens up new possibilities for further investigating the consequences of this separation.

This classification of rainfall origin, which is also the contribution of the proposed method in comparison with other similar studies using OKRE, is implemented in the real-time QPE estimation process. The present study shows that QPE in convective events involved larger errors than for stratiform events. This finding may, in part, be explained by the greater spatial irregularity of convective precipitation, corroborated by the modeled variogram ranges. Moreover, spatial correlation for convective events can be characterized up to ∼4 km, while for stratiform events up to ∼5 km, both values taken from the median ranges of all fitted models.

Further research could involve using more sophisticated algorithms such as KED and parametric variogram modeling. Another future application of these variogram models would be the simulation or disaggregation of convective and stratiform events. This study has provided the basis for subsequent hydrological calibration and validation based on rainfall runoff modeling, while also offering a methodology that can be applied in basins with similar characteristics.

X-band radar installation was funded by the FEDER project UNCM08-1E-086. This study also forms part of the MARCoNI project (funded by the Spanish Ministry of Economy and Competitiveness, ref: CGL2013-42728-R), the MIDHATO Venero project (funded by the Spanish Geological Survey, ref: IGME 2313/2010) and projects CGL2013-48367-P and CGL2016-80609-R (funded by the Spanish Ministry of Economy and Competitiveness). We would also like to thank technicians and rangers from the Regional Government of Castile and León, Navaluenga City Council and the Fundación Caja de Ávila, for their essential support during the study.