Spatial interpolation of hourly rainfall – effect of additional information , variogram inference and storm properties

Abstract. Hydrological modelling of floods relies on precipitation data with a high resolution in space and time. A reliable spatial representation of short time step rainfall is often difficult to achieve due to a low network density. In this study hourly precipitation was spatially interpolated with the multivariate geostatistical method kriging with external drift (KED) using additional information from topography, rainfall data from the denser daily networks and weather radar data. Investigations were carried out for several flood events in the time period between 2000 and 2005 caused by different meteorological conditions. The 125 km radius around the radar station Ummendorf in northern Germany covered the overall study region. One objective was to assess the effect of different approaches for estimation of semivariograms on the interpolation performance of short time step rainfall. Another objective was the refined application of the method kriging with external drift. Special attention was not only given to find the most relevant additional information, but also to combine the additional information in the best possible way. A multi-step interpolation procedure was applied to better consider sub-regions without rainfall. The impact of different semivariogram types on the interpolation performance was low. While it varied over the events, an averaged semivariogram was sufficient overall. Weather radar data were the most valuable additional information for KED for convective summer events. For interpolation of stratiform winter events using daily rainfall as additional information was sufficient. The application of the multi-step procedure significantly helped to improve the representation of fractional precipitation coverage.


Introduction
Precipitation data with a high resolution in space and time are the driving forces for hydrological modelling of floods.While the temporal resolution of the recording stations is suitable, the network density is often too sparse for a reliable spatial interpolation.So, the more variable convective rainfall in spring and summer is captured worse than the less variable stratiform winter precipitation (Kalinga et al., 2003).Precipitation data from non-recording stations are usually provided in a dense network, but have only a daily temporal resolution.Meanwhile radar data have been used more frequently as input for hydrological modelling due to its advantage of the high spatial resolution.A considerable number of studies have shown though, that uncorrected radar data are insufficient due to the frequently occurring large space-time variable bias (Cole and Moore, 2008;Ehret et al., 2008).To obtain high space-time resolution fields of precipitation for flood studies it is therefore necessary to apply sophisticated interpolation methods on the short time step rainfall data in combination with radar information and other available additional information.
Several mapping techniques for rainfall fields like ordinary kriging or spline-surface fitting have been in use over some time (Creutin and Obled, 1982;Dubois et al., 2003).They provide the basis for the application and further development of multivariate geostatistical methods like kriging with external drift or collocated co-kriging using various co-variables (Sarangi et al., 2005;Grimes et al., 1999;Lloyd, 2005;Carrera-Hernandez and Gaskin, 2007).Goovaerts (2000) uses annual and monthly rainfall observations to demonstrate, that the three multivariate geostatistical methods, all incorporating a digital elevation model, outperform three other univariate methods.Kyriakidis et al. (2001) considers seasonal average daily precipitation showing that the Published by Copernicus Publications on behalf of the European Geosciences Union.

570
A. Verworn and U. Haberlandt: Spatial interpolation of hourly rainfall -effect of additional information integration of atmospheric and terrain characteristics in a geostatistical framework leads to more accurate representations of the spatial distribution of rainfall.
Due to the better availability and improved accuracy of radar data a new sector of radar-raingauge merging methods has been evolved.Velasco-Forero et al. (2009) compare three geostatistical interpolation methods, all incorporating radar data as secondary information in combination with a non-parametric technique to automatically compute correlograms.Kriging with external drift shows the most accurate results.Ehret at al. (2008) have developed a merging method combining a mean precipitation field interpolated from rain gauge observations with information about the spatial variability from radar data.García-Pintado et al. (2009) combine a multiplicative -additive decomposition and an objective analysis scheme to estimate multisensor rainfall.Goudenhoofdt and Delobbe (2009) have evaluated several radargauge merging methods with various degrees of complexity and favour the geostatistical merging methods.Still several questions need further adressing for an improvement of interpolation like sufficient consideration of fractional coverage of rainfall occurrence (Seo, 1998), importance of preprocessing of radar data, value of topography for short time step interpolations, relevance for discriminating seasons and storm types when chosing interpolation methods etc.
Prior to the interpolation task, the spatial persistence structure of precipitation has to be analysed usually based on semivariogram analysis (Holawe and Dutter, 1999;Skoien and Blöschl, 2003;van de Beek et al., 2010).Especially in mountainous regions it is critical to determine semivariograms based on a sparse raingauge network only.Germann and Joss (2001) report it to be less critical if using evenly spaced radar data.As an alternative the application of a 3-D estimation of the variogram with rainfall duration as third coordinate is suggested (Bargaoui and Chebbi, 2008).This leads to significantly lower prediction error than the classical 2-D kriging.Kravchenko (2002) has shown that the interpolation performance of soil properties with kriging and known variogram parameters usually outperforms the case when variogram parameters are unknown.Furthermore, optimal sampling schemes for a minimal kriging variance are considerably influenced by variogram parameters ( van Groenigen, 2000).Yet the influence of the semivariogram estimation on the interpolation performance of short time step rainfall is relatively unknown.
This research extends on a case study of Haberlandt (2007) about hourly rainfall interpolation using rain gages and radar for one extreme rainfall event.Here, based on a set of 15 rain storms, investigations focus first on a detailed analysis of various approaches for variogram inference with differentiations between seasons and storm types.Then, several interpolation methods are compared using additional information from weather radar, topography and from the daily raingauge network.Special attention is given to consider the fractional rainfall coverage, to analyse effects of radar data preprocessing and relating the results to the different storm types and seasons.This study will provide the basis for a subsequent hydrological validation based on rainfall runoff modelling using the differently interpolated precipitation input fields.
This paper is organised as follows.After the "Introduction" the section "Methods" follows with a description of the semivariogram estimation techniques and a description of the interpolation method kriging with external drift (KED) considering the special cases which were applied here.The study area and the available data will be introduced in section three.Also, the processing of the data, especially the radar observations, is described here.Section four discusses the results and is divided into two subsections dealing with variogram inference results and performance evaluations of the applied interpolation methods using different additional information.In the final section the main findings are concluded and an outlook is presented.

Semivariogram estimation
The spatial structure of rainfall mainly depends on weather conditions, type of the precipitation, topography of the study area, spatial and temporal scale, and can be quite dynamic in space and time.Therefore, one particular problem for geostatistical interpolation of complete time series is the effective and reliable estimation of the semivariograms for each time step.
A semivariogram measures the spatial variability of a regionalized variable Z assuming the variable being stationary or intrinsic (Armstrong, 1998).The traditional experimental semivariogram is defined as follows: where N (h) is the number of data pairs, which are located a distance vector h apart.The fitting of a theoretical model is necessary in order to deduce semivariogram values for any possible lag h required by interpolation algorithms (Goovaerts, 1997).The spherical model with a nugget component is chosen for the investigations, showing a linear behaviour near the origin: where a is the range, c the sill and c 0 the nugget.The exact matching of the experimental and theoretical variogram is considered less important than the choice of variogram model itself (Wackernagel, 1995).Theoretically, the semivariogram estimation needs to be carried out separately for each time step.If this is done manually, the procedure will be very time-consuming.Therefore, Hydrol.Earth Syst.Sci., 15, 569-584, 2011 www.hydrol-earth-syst-sci.net/15/569/2011/ A. Verworn and U. Haberlandt: Spatial interpolation of hourly rainfall -effect of additional information 571 an automatic approach was applied in addition to some averaging techniques, which are explained in the following: -Event-specific semivariograms: For each event a specific experimental semivariogram was obtained by averaging over all time steps, which exceed a specific threshold of precipitation intensity.
Before averaging, the experimental semivariogram was standardized by the variance for each time step: where n is the number of time steps, γ (h,i) is the semivariogram value for the distance class h of time step i and var(i) is the variance of time step i.Afterwards a spherical variogram model was fitted manually.Additionally, two cases were distinguished: One with assumed isotropy ("event isotropic") and the other one where spatial anisotropy was taken into consideration ("event anisotropic").The zonal anisotropy was modelled with nested spherical structures and an anisotropy factor, respectively (Deutsch and Journel, 1992).The term "zonal" means, that not only geometric anisotropies with different ranges but same sills were taken into consideration, but also more complex structures with sills not being the same in all directions.
-Average semivariograms: Based on the "event isotropic" type a total average semivariogram over all events ("average isotropic") was derived by weighting the parameters of the theoretical semivariograms by the number of time steps of each event which were used for averaging: where θ k stands for the parameter k of the theoretical semivariogram of event i, m is the number of events and n(i) the number of time steps of event i.Furthermore, an averaged semivariogram was calculated for each season (summer and winter) separately ("seasontype isotropic").
-Automatic semivariograms: An automatic fitting method ("automatic isotropic") was applied for each time step, such that a weighted sum of squared differences between the experimental and the applied spherical semivariogram model was minimized (Cressie, 1985): where k is the number of data pairs for each lag l and time step i.The minimisation of Eq. ( 5) was achieved with the Nelder and Mead optimisation method (see e.g.Press et al., 1989) assuming that the variance is equal to the sum of nugget c 0 and sill c.This has been done with the objective to provide sufficiently robust results.This way only the two parameters range and ratio of nugget to sill are optimised, instead of three parameters range, sill and nugget.The nugget was forced to zero only, if the calculated value was below zero.The sill was then set equal to the variance.For each time step the obtained objective function value for the automatic variogram was compared with a value considering the averaged variogram instead and that variogram with the lower objective function value was selected.Also, if no convergence could be reached, the averaged semivariogram was chosen and scaled with the variance.
-Assumed-linear semivariograms: As the simplest version without using any data, a linear isotropic semivariogram with γ (h) = h was applied.

Kriging with external drift
For interpolation, kriging with external drift (KED), ordinary kriging (OK) and inverse distance weighing (IDW) were applied.Only KED will be described in the following.OK and IDW were used as reference for comparisons.Kriging with external drift (KED) is a simple and efficient algorithm allowing the incorporation of one or more secondary variables, which are assumed to be linearly related to the expected value of the primary variable.As opposed to the method "kriging with a trend model", the additional information is not modelled as a function of the coordinates, but is treated as a smoothly varying linear function according to the denotation as "external variable" (Goovaerts, 1997).
The KED estimator for the unknown point consists of a weighted sum of n observed points in the neighbourho "ordinary kriging" (Webster and Oliver, 2001): where λ i are the kriging weights.For KED it is assumed that the expected value of Z is linearly related to a number of m additional variables Y k (u) as The unknown parameters b k are estimated implicitly and do not appear in the kriging system: where n is the number of neighbours, m is the number of additional variables Y k included, while µ k are Lagrange multipliers.The variogram values γ are inferred simply from the original variable Z and not from the trend residuals.Using the semivariogram of original rainfall data instead of residuals for KED is certainly a simplification.The difficulty lies in simultaneously estimating the unknown trend and the residuals.To overcome this problem, the trend component could be computed with a slightly modified KED system (Deutsch and Journel, 1992), while calculating the residuals and the residual semivariogram afterwards.However, this iterative process would be very demanding.Another simpler approach could be to infer the semivariogram only from data pairs, that are unaffected by the trend.Investigations of Haberlandt (2007) have shown though, that there were no significant differences in interpolation performance, which was the reason, why the simplified approach is applied here as well.A more detailed description of external drift models can be found in Chilès and Delfiner (1999).
One main focus of this study was the refined application of kriging with external drift using various additional variables.In this context the following versions were implemented: -Due to the assumption of a linear relation between primary and secondary variable, an appropriate transformation of the latter could be useful in the form of logarithm or square-root in case of a non-linear relations.
-In order to avoid instabilities in the KED system, a socalled "conditional" version of KED was introduced.
Instabilities mostly occur at time steps, where rain is detected only at very few stations, and, therefore, the relation between primary and secondary variable is weak.
To overcome this problem, KED was only executed, if the correlation between both variables exceeded a certain threshold.Otherwise, OK was used for those time steps.
-In addition, a "multi-step" interpolation procedure was set up to consider sub-regions without rainfall in a more appropriate way.First binary indicator kriging was applied as follows: The precipitation time series of all recording stations were transformed into zero's for no rain and one's for rain.Afterwards OK was carried out on the binary time series to determine rainy and no-rainy cells.Estimated cell values below a threshold of 0.5 were counted as dry and were set to zero.Cell values equal or above of 0.5 were counted as wet and set to one.As a second step KED or OK was applied on the original time series.For obtaining the final estimate, the interpolated fields from both steps were multiplied, which resulted basically in interpolating just over the wet areas.However, it has to be mentioned that this procedure does not fully preserve the precipitation volume and therefore possibly produces a bias.
Finally, different combinations of various additional variables were investigated, given that KED allows the incorporation of more than one drift variable.This procedure was carried out under the assumption, that some additional information might only be valuable in concurrent use with other variables.

Performance assessment
The impact of the semivariogram inference on the interpolation performance as well as the interpolation performance itself was evaluated by precipitation cross-validation.The principle of the so called "leave-one-out-method" is to estimate rainfall successively for each sampled location using the known neighbours but always discarding the sample value at the particular location.An estimate was therefore obtained at every sample point, which was compared with the observed values using the following performance measures.
In addition to the simple bias criterion the root-mean-square-error (RMSE), standardized with the observed average Z, was applied where Z * (u i ) is the estimated value and Z(u i ) the observed value, both at location u i .The RVar coefficient measures the preserved variance of the interpolated values in relation to the variance of the observed values According to the known but unwanted smoothing character of precipitation interpolation, it was aimed for an RVar value close to 1, preserving the observed variance as much as possible.3 Study area and data

Study regions
The 125 km radius around the radar station "Ummendorf" in the northern part of Germany comprises the total study area including 21 recording rainfall gauges.Within the operating radius of the radar also 676 non-recording rainfall gauges are located providing daily precipitation values.The terrain is characterized mainly by flat land in the north-eastern part, while the south-eastern part is strongly influenced by the Harz Mountains and its foothills (Fig. 1).The elevation ranges from almost sea level to 1043 m a.s.l.(above sea level).The total study area was used for areal rainfall estimation and precipitation cross-validation.According to Fuchs et al. (2003), the mean annual precipitation for the study area varies between around 400 mm yr −1 on the leeward side of the Harz Mountains and about 1300 mm yr −1 at the mountain tops.This variability is caused by the windward/rain shadow effect with moist air masses usually coming from the west, leading to increased precipitation on the western side of the mountains, while on the eastern side there is a large drop in precipitation.

Data and pre-processing
The general idea was to carry out investigations on an event basis.The rainfall events were selected based on the inspection of flood hydrographs observed at streamflow gauges in the study area.Our intention was to validate the interpolation methods using hydrological modelling in a subsequent investigation.The recording stations provide data with 10min time steps, which were aggregated to hourly values.Events occurring in different seasons and caused by different weather conditions were considered.Both convective and stratiform precipitation structures might have contributed to one event.The selected events as well as some statistics are listed in Table 1.For the calculation of standard deviation and coefficient of variation (COV) only those time steps were used which exceed an average precipitation of 0.1 mm h −1 .This low threshold was used for the event statistics in order to exclude time steps with non-significant rainfall.A classification of the events into convective and stratiform storms was made based the COV and season.A COV value above 1.70 indicated rather convective precipitation, equal or below 1.70 a more stratiform event.Rather convective precipitation usually occurs during summer.The threshold of 1.70 was established based on visual inspections of specific radar images for all events.Three different types of additional information for interpolation using kriging with external drift (KED) were utilised here.The derivation of those data is explained in the following.
The first type consists of the two time-invariant external drift variables elevation and luv/lee index, both derived from a digital elevation model (DEM) with a resolution of 1 km × 1 km and a fixed window size of 230 km × 230 km.The determination of the luv/lee index was based on the direction of slope and wind.The slope direction was directly inferred from the digital elevation model while the wind direction was taken from the climate station "Brocken" (Fig. 1) providing daily observed wind data.The wind direction was  assumed to be spatially constant over the study area.The spatial resolution of the luv/lee index was of special interest though.The precipitation occurrence is usually affected by larger-scale climatologic characteristics.A small DEM resolution was not appropriate in this case, given that it would indicate luv cells in an overall lee area and vice versa through its dependence on the slope direction.To find the optimal resolution of the digital elevation model regarding derivation of topographic indices, the cellsize was changed, but with the condition of keeping the same extent.This procedure resulted in a cell size of 5.75 km × 5.75 km for the luv/lee index showing the highest correlation to observed precipitation.The subsequent downscaling of the index to the 1 km × 1 km resolution for the precipitation interpolation task was carried out with a bilinear resampling technique.
The second type of additional information concerns timevariant daily rainfall data from the non-recording stations.Three different versions were distinguished here: (a) rainfall data were used on the daily basis as measured (P daily ), (b) daily rainfall values were accumulated up to the current interpolation time step (P cum ), and (c) rainfall data were aggregated to event sums (P event ).The data were spatially interpolated with the inverse distance method (IDW).Another alternative would have been the application of ordinary kriging (OK).However, due to the dense network of daily stations the differences were expected to be low and the simpler procedure has been applied here.
As third type of additional information for KED the highly dynamic time-variant radar data were employed.The high spatial resolution of radar data is most valuable in combination with KED.For processing of radar data it should be noted that those data were used here only as background field and not independently as primary rainfall information.A good correlation between hourly rainfall from the recording stations and radar information was therefore more important than an optimal adjustment or unbiased estimation of radar rainfall.
Radar observations from the C-band instrument at "Ummendorf" were provided as raw reflectivities with a spatial polar resolution of 1 km × 1 km azimuth and a time discretisation of 5 min (dx-product of the German Weather Service, DWD).A statistical clutter filter was applied by the DWD on the raw data already.Prior to a transformation into rainfall intensities missing radar values were filled by interpolation and a circle with clutter close to the location of the radar instrument was cut out from the study area.Afterwards, the reflectivities were transformed into rainfall intensities applying the well-known Z − R relationship: where Z is the reflectivity in mm 6 m −3 and R is the rainfall intensity mm/h.The parameters were set to a = 200 and b = 1.60 constantly for each time step and for all events according to the standard Marshall-Palmer Z − R relationship.
Both reflectivities and rainfall intensities were used here as additional variables for KED.The variables were interpolated on a regular raster with a resolution of 1 km × 1 km as follows.If more than one radar observation in polar coordinates belongs to a raster cell, the radar data were averaged.In case no radar point falls into a raster cell, which occurs at further distances from the radar origin, the value of the nearest neighbour was allocated to the particular raster cell.
Table 2 shows estimated average correlations between rainfall at all recording stations and radar rainfall at the corresponding cells for each event.Correlations range from 0.28 to 0.93 depending on the specific event with an overall average of 0.70.This indicates that radar data could be expected to be a useful additional variable for KED.
Using radar data directly as rainfall information an attenuation correction is strongly recommended (Krämer, 2008;Illingworth, 2004).When using radar data as background field for interpolation as applied in this case, the importance of attenuation correction is not so clear.In order to evaluate the influence of attenuation correction, a simple test was carried out as follows.In addition, the transformations into radar rainfall intensities were carried out with a uniform non-linear attenuation correction for all time steps and for all events according to Krämer (2008).The correlation between precipitation at all recording stations and attenuation corrected radar rainfall intensities was calculated and is shown in Table 2.The effect of a uniform non-linear attenuation correction on the correlations was negligible.So, it could be omitted when using KED for interpolation.

Variogram inference and impact on interpolation
Semivariograms were inferred from radar data rainfall due to the higher resolution in space compared to the recording stations.One thousand randomly selected radar cells based on the 1 km × 1 km grid were used for the estimation of experimental semivariograms.Comparative simulations have shown that this number represents the spatial precipitation structure sufficiently.In order to put more weight on time steps with significant precipitation, the estimation of experimental semivariograms was applied only for time steps exceeding an average radar rainfall threshold of 0.5 mm h −1 .For some smaller events, the threshold was lowered to 0.25 mm h −1 .With higher thresholds not enough time steps would have been available for some events due to the systematic rainfall underestimation by radar.
Precipitation cross-validation was carried out for each event evaluating the influence of the different variogram types as described in Sect.2.1.The methods ordinary kriging (OK) and kriging with external drift (KED) in their original versions were applied for this interpolation task.For the latter only the radar rainfall intensities were used as additional variable.The root mean square error standardized with the observed average (RMSE) served as a performance indicator.The cross-validation exercise was applied for all time steps with an average precipitation intensity of the recording stations exceeding a rainfall threshold of 1.0 mm h −1 .In order to still focus on time steps with heavy precipitation, the threshold had to be raised here in comparison to the case when using radar data alone.Event-specific isotropic experimental semivariograms were estimated for each event.The manually estimated parameters nugget, range and sill for a spherical model are listed in Table 3.The detected low nugget effects are probably related to the used dense network of radar data considering one thousand randomly selected radar cells.Furthermore, averaged variogram parameters for summer and winter as well as in total were calculated.Although the averaging of theoretical variogram parameters is rather unusual it provided a straightforward way to estimate mean variograms.Comparisons have shown, that a mean experimental variogram calculated over all events matches well with a mean theoretical variogram based on averaged parameters.In Fig. 2 event-specific isotropic experimental and theoretical semivariograms for two selected precipitation events of different types are presented.The summer events and convective storms showed generally shorter range and higher sill compared to winter events and frontal storms.For instance, the range of the frontal rainfall event no.14 was with a = 115 km more than twice as large as the range of the convective summer storm event no.16 with a = 45 km.The low nugget effect is clearly visible for both events.Comparing semivariograms for different directions, zonal anisotropy was visible for all events.Figure 3 shows an example, where the higher sill has an azimuth angle of 90 • , while the greater range was found in the north-south direction.This behaviour corresponded to the general weather pattern with storms usually moving from west to east over the region.A comparison of interpolation performance for the different variogram types is shown in Figs. 4 and 5 using crossvalidation based on the methods OK and KED.Generally, the choice of method to calculate the semivariogram had only a small impact on prediction performance.It could even be noticed that for some events interpolation performs better when using an "assumed linear" variogram compared to specifically estimated variograms.Nevertheless, for the major part of events the specific variograms were the better choice.Being more precise regarding the differences between the selected variogram types, it was difficult to identify which variogram type shows the best performance as it varied from event to event.This was the case for both interpolation methods OK and KED.In addition, the relative performance of a specific variogram often changed with the applied interpolation method.A comparison of averaged RMSE values over all events is shown in  4. RMSE difference between the selected variogram type and an "assumed linear" variogram, averaged over time steps with mean precipitation intensity exceeding a threshold of 1.0 mm h −1 .OK is used for interpolation.specific variogram type was the best overall, since the absolute RMSE differences were very small (Table 4).The use of an averaged isotropic semivariogram seemed therefore sufficient, while the proposed automatic estimation method did not lead to further improvements.Figure 6 confirms the results from the cross-validation showing the spatial distribution of precipitation for one hour on the 17 July 2002 interpolated with OK based on the diverse semivariograms.Using specific semivariograms the spatial pattern looked similar, while an 'assumed linear' type produced a different and smoother map.

Interpolation using different additional information
Spatial interpolation of hourly precipitation was carried out with the average isotropic semivariogram (see Table 3).Here, the focus was on the evaluation of different versions of the multivariate geostatistical method kriging with external drift (KED) as described in Sect.2.2.The univariate interpolation methods inverse distance weighting (IDW) and ordinary kriging (OK) were used as reference.The interpolation performance was assessed on the basis of cross-validations (see Sect. 2.3) including all 15 events, but only time steps exceeding a threshold 1.0 mm h −1 of average precipitation intensity as for the variogram cross-validations.The performance indicators RMSE and RVar were averaged over summer and winter events weighted by the number of time steps for each event.The interpolation variants were separated into four cases: (a) reference methods IDW and OK, (b) KED interpolation without using weather radar as additional information, (c) KED interpolation using radar reflectivities as central additional information and (d) KED interpolation using radar rainfall as central secondary variable.Considering the conditional version of KED, tests showed that a correlation threshold of 0.5 for interpolation without radar data  and 0.3 for interpolation with radar data produced the lowest RMSE, which were therefore applied here.Higher thresholds reduced the number of cases when KED is applied, while at lower thresholds instabilities in the KED kriging system prevented better results.Theoretically KED should be equivalent to OK if there is no correlation between primary and secondary variables.However, experimental evidence showed that for some time steps with low correlation KED resulted in significant worse performance than OK.Applying the conditional versions this could be avoided.Results of precipitation cross-validation are shown in Table 5 for summer events and in Table 6 for winter events.Selected results are visualised in Figs.7 and 8, which are discussed in detail as follows.
Using elevation as additional information for KED generally did not show a reduction of the RMSE in comparison with OK.Neither a transformation in the form of logarithm or square-root of the additional variable nor interpolation in the conditional version improved the interpolation performance compared to OK for summer events.Only in the conditional version for winter events a marginal improvement was detected (Table 6, B1).A conclusion which types of rainfall benefit from KED interpolation with elevation, could not be made in this study.One reason for the weak correlation between precipitation and elevation is certainly the short hourly time step.Similar results were obtained for summer using the luv/lee index as additional information, while for winter storms a significant reduction of the RMSE was achieved (Table 6, B2).
The most valuable additional information for interpolation without radar came from the daily precipitation network favouring daily rainfall amounts, followed by cumulated multi daily rainfall and event rainfall.Either version was clearly preferred to the other additional variables no matter for which season.The RMSE could further be reduced by the logarithmic transformation of P daily in combination with conditional KED, but mainly for winter events.However, the variance was preserved significantly less (Tables 5 and 6, B7).The best performance was achieved with the application of the multi-step procedure.Due to the better discrimination of regions with and without rainfall, the RMSE was smaller no matter what additional variable was used.In addition, the variance of the observed rainfall was preserved better.The disadvantage of producing a negative bias using this method (see Sect. 2.2) was not a problem here as can be seen from Tables 5 and 6.Only for the cases D using radar rainfall as additional variable a small negative bias occurred.Even if time steps with rainfall smaller than 1.0 mm h −1 were included in cross-validations, the negative bias did not increase.However, regarding spatial interpolation, a loss of total precipitation volume needs to be considered.The use of more than one additional variable for cases without radar showed no improvement compared to the application of each additional variable alone.This was true for both seasons.
The fact that the elevation did not contribute to a better interpolation performance was related possibly to the short time step and to the topography of the study area.Over three quarters of the recording stations were located in flat terrain and the average correlation between rainfall and elevation or luv/lee index was only 0.32 and 0.16, respectively.The daily stations were available from a dense network showing a stronger average correlation of 0.61 to the hourly rainfall data.However, the daily rainfall data are not available for real-time applications like operational flood forecasts.
Comparing summer and winter events up to this step the relative reduction of the RMSE with the methods above turned out to be similar.However, the absolute values of the RMSE were almost twice as high for summer events.The RVar indicated a higher preserved variance of the observed values during winter.
Finally, radar data were employed as additional information for KED in form of reflectivities as well as rainfall intensities (see Sect. 3.2).Applying raw reflectivities as only secondary information for summer events, a significant improvement was made if rainfall as primary and reflectivities as secondary variable were both log transformed (Table 5, Hydrol. Earth Syst. Sci., 15, 569-584, 2011 www.hydrol-earth-syst-sci.net/15/569/2011/  C1).This leads to a linearization of the Z − R relationship (Eq.12) and satisfies the basic assumption of KED (Eq.7).
In case of interpolation on log-transformed values y, the estimated values were back-transformed taking exp(y) before calculating performance measures.With interpolation on log-transformed values, problems with biased results may occur.The exponentiation in the back-transformation process exponentiates any errors, which is very sensitive to interpolation outliers.However, in this case the results seemed not to be affected by that as the cross-validation showed values of the bias in the same magnitude as for interpolation methods without log-transformed data.In case of winter events (Table 6) the use of reflectivities as additional information did not lead to an increase in performance compared to using daily rainfall.A further reduction of the RMSE could be achieved if P daily was added as another additional variable (Tables 5 and 6, C4).Elevation led to a higher preserved variance, but only in winter (Table 6, C3).The conditional version was without effect if all variables were log transformed (Tables 5 and 6, C5).
Using rainfall intensities as additional information for KED (Table 5, D1), the RMSE was a bit larger than for the case with log-log transformed reflectivities (Table 5, C5), while the indicator for preserved variance rose.In case of winter events (Table 6) using rainfall intensities likewise an increase in performance was not detectable in comparison to the application of daily rainfall as external drift.Adding another secondary variable besides radar, only daily rainfall led to a further reduction of the RMSE and a significant higher preserved variance of the estimated values, which was the case for both seasons (Tables 5 and 6, D6).Neither elevation nor the luv/lee index was valuable in combination with radar data in this case.Solely the variance of the estimated rainfall was higher with elevation or the luv/lee index as additional  stratiform precipitation structures radar data were not really necessary.Almost the same performance could be achieved using KED with the daily rainfall as drift variable.
In addition, the potential value of using radar data directly as rainfall without calibration was evaluated.For that no cross-validation is carried out, but a direct comparison of radar rainfall with observed gauge rainfall (case E1 in Tables 5 and 6 as well as in Figs.7 and 8).The direct utilisation of radar rainfall showed the worst performance, especially for stratiform winter events.The negative bias revealed a strong underestimation of observed rainfall.
Comparing the RMSE values from Tables 5 and 6 it can be seen, that for summer events the RMSE is significantly larger than for winter events.The large RMSE values, especially in summer, show the overall high uncertainty in hourly precipitation interpolation.
In Fig. 9 maps are presented showing the spatial distribution of hourly precipitation interpolated using selected methods exemplarily for one hour on the 17 July 2002.While the application of OK (A2) shows quite a smooth map as expected, the addition of the daily rainfall (B3) did not significantly improve the representation of the spatial precipitation pattern.This was related to the low correlation between primary and additional variable (0.24) for this time step.Only when radar data were used as additional information for KED, the complex spatial structures of precipitation appeared in the maps.For the log-log transformed case with reflectivities (C1) the local extremes were most pronounced which might be a negative artefact of the transformation.5).
comparison to the case without elevation (D1).If radar rainfall intensities and daily rainfall were employed in combination with the multi-step procedure (D8), the no rainfall fraction in the region was largest, which is in most cases desirable.Geostatistical interpolation methods usually lead to a spatial smoothing of the rainfall distribution with underestimation of high values und overestimation of low values.Especially, no rain areas become rainy areas with low precipitation intensities after interpolation.These smoothing effects are not wanted for flood simulation.

Summary and conclusions
In this study first the effect of using different approaches for estimating semivariograms on the spatial interpolation performance of hourly precipitation was assessed.Furthermore, different versions of kriging with external drift (KED) using additional information from physiographic factors, from the daily rainfall measurement network and especially from weather radar were investigated.The main results and conclusions can be summarised as follows: -The impact of different semivariogram types on the precipitation interpolation performance was low.Nevertheless, the performance varied over the events, whereas the event-specific types showed on average no improvements in comparison with average types.However, the worst performance was found for the simple "assumed linear" type.Using an average semivariogram seemed therefore sufficient for this kind of interpolation.
-Weather radar data proved to be the most valuable additional information for KED for convective summer events.When using reflectivities both primary and secondary variable should be log transformed.A uniform non-linear attenuation correction applied on the raw radar data did not improve the interpolation performance.
-For stratiform winter events daily rainfall as additional information was sufficient, while radar data did not significantly contributed to a better interpolation performance.Generally, in absence of radar data the use of rainfall from the daily network as secondary variable was recommended.
-Using KED with additional variables like elevation or luv/lee index could hardly improve the interpolation performance according to the cross-validation results compared with the univariate method ordinary kriging (OK).However, when including topography, spatial distributions of precipitation showed sharper differentiations of spatial structures, which seemed more plausible.So, the use of additional information from topography can generally be recommended.
-Applying log-transformations only on the additional variables except for radar reflectivities the interpolation performance could be further improved, although the preserved variance decreased.
-The application of the multi-step procedure significantly contributed to a better representation of fractional precipitation coverage.
This study was conducted on the basis of 15 flood events caused by precipitation of different characteristics with hourly discretisation in time.It is expected, that the main findings regarding precipitation interpolation performance will generally hold for other events and also for similar regions.For mainly mountainous regions additional information from topography and climatology might have a greater influence on the precipitation interpolation performance as in this case.In subsequent work investigations should extend to other regions with different characteristics as well as to analyses with different time resolutions.The hydrological validation of the interpolation performance using rainfall-runoff modelling is in progress for selected catchments within this study area and will be reported elsewhere.

Fig. 5 .
Fig.5.RMSE difference between the selected variogram type and an "assumed linear" variogram, averaged over time steps with mean precipitation intensity exceeding a threshold of 1.0 mm h −1 .KED with radar rainfall intensities is used for interpolation.

Fig. 7 .
Fig. 7. Precipitation cross-validation results using all summer events averaged over time steps with average precipitation exceeding a threshold of 1 mm h −1 for selected interpolation methods (red bars: univariate methods; yellow bars: KED without radar; green bars: KED with radar reflectivities; blue bars: KED with radar rainfall intensities; grey bars: RVar coefficient, black bar: radar rainfall intensities only).

Fig. 8 .
Fig. 8. Precipitation cross-validation results using all winter events averaged over time steps with average precipitation exceeding a threshold of 1 mm h −1 for selected interpolation methods (red bars: univariate methods; yellow bars: KED without radar; green bars: KED with radar reflectivities; blue bars: KED with radar rainfall intensities; grey bars: RVar coefficient; black bar: radar rainfall intensities only).

Fig. 9 .
Fig. 9. Spatial distribution of precipitation in mm/h for one hour between 06:00 p.m. and 07:00 p.m. on the 17 July 2002 interpolated with different methods (variant number see Table5).

Table 1 .
Precipitation statistics of the selected events with hourly time step data; standard deviation, coefficient of variation and no-rain fraction are spatially calculated and averaged over time.

Table 2 .
Correlation between precipitation at all recording-stations and radar rainfall intensities for time steps exceeding a threshold of 1.0 mm/h of average precipitation.

Table 3 .
Estimated semivariogram parameters considering a spherical model for each event, averaged over seasons and averaged over all events (sill is standardised with variance; only time steps with radar rainfall intensity >1 mm h −1 are used for variogram estimation).

Table 4 .
Interpolation performance from cross-validation of the various variogram types with respect to RMSE criterion, averaged over all events and time steps exceeding a threshold of 1.0 mm h −1 of average precipitation.

Table 5 .
Precipitation cross-validation results using all summer events averaged over time steps with an average precipitation exceeding a threshold of 1 mm h −1 ("log": all Y k are log-transformed, "log-log": both Z and all Y k are log-transformed; highlighted no. are presented in Fig.7).