Regionalizing non-parametric precipitation amount models on different temporal scales

Parametric distribution functions are commonly used to model precipitation amounts corresponding to different durations. Precipitation amounts themselves are crucial for stochastic rainfall generators or weather generators. Non-parametric kernel density estimates (KDEs) offer a more flexible way to model precipitation amounts. As it is already stated in their name, these models do not exhibit parameters which can easily be regionalized to run rainfall generators at ungauged locations as well as at gauged locations. To overcome this deficiency we present a new interpolation scheme for non-parametric models and 5 evaluate it for different temporal resolutions ranging from hourly to monthly. During the evaluation the non-parametric methods are compared to commonly used parametric models like the two-parameter gamma and the mixed exponential distribution. As water volume is considered to be an essential parameter for applications like flood modeling, a Lorenz-curve based criterion is also introduced. To add value to the estimation of data at sub-daily resolutions, we incorporated the plentiful daily measurements in the interpolation scheme and this idea was evaluated. The study region is the federal state of Baden-Württemberg in 10 the southwest of Germany with more than 500 rain gauges. The validation results show that the newly proposed non-parametric interpolation scheme works, and additionally seems to be more robust compared to parametric interpolation schemes, and that the incorporation of daily values in the regionalization of sub-daily models is very beneficial.


Introduction
Rainfall time series of differing temporal resolutions are needed for various applications like water engineering design, flood modeling, risk assessments or ecosystem and hydrological impact studies (Wilks and Wilby, 1999;Burton et al., 2008).As many precipitation records are too short and contain erroneous measurements, stochastic precipitation models can be used to generate synthetic time series instead.Starting from single-site models (summarized in Wilks and Wilby, 1999), multi-site models for simultaneous time series at various sites (e.g., Wilks, 1998;Buishand and Brandsma, 2001;Bárdossy and Plate, 1992) and finally models which allow for gridded simulations are developed (e.g., Wilks, 2009;Burton et al., 2008).
For modeling precipitation one crucial variable is the precipitation amount, which follows a certain distribution.Distributions of daily precipitation amounts are strongly right skewed, with many small values and few large values (Wilks and Wilby, 1999;Li et al., 2012;Chen and Brissette, 2014).This also holds true for different temporal resolutions with increasing skewness for higher temporal resolutions and vice versa.This means that rainfall intensity distributions depend on the temporal scale of the observed values.Applying single-site or multi-site precipitation models at ungauged locations requires regionalization of precipitation amount distributions.This can be done in two different ways: 1. Interpolate precipitation amounts from observation points for every time step to the target location(s) and set up a distribution with the interpolated values.
2. Fit a distribution function to the precipitation amounts separately for each gauge and interpolate the distribution functions to the target location(s).
The first approach seems more straightforward, but exhibits several deficiencies such as overestimation of the rainfall probability, underestimation of the variance and underestimation of the maximum rainfall value.In supplement section S1 an example demonstrates these problems.Due to the relative inefficiency of the first interpolation approach, the second is preferred.
In most stochastic rainfall models theoretical parametric distribution functions are fitted to the empirical values using, e.g., the exponential distribution or the two parametric gamma distribution (Wilks and Wilby, 1999;Papalexiou and Koutsoyiannis, 2012).It is possible to either interpolate the parameters of the theoretical distribution or to interpolate the moments (e.g.mean and standard deviation) of the rainfall intensities (Wilks, 2008;Haberlandt, 1998).Lall et al. (1996) introduced a more flexible non-parametric single-site rainfall model, where they used non-parametric KDEs with a prior logarithmic transformation to model daily rainfall intensities.They mentioned the problem of regionalization by using non-parametric estimates of distribution functions.However, a different interpolation scheme is required for non-parametric estimations, as they do not use any parameter that can be simply interpolated.
In the present work we introduce a regionalization strategy for non-parametric distributions and compare it to the traditional regionalization of parametric distributions for varying temporal resolutions from hourly to monthly scale.The common procedure to interpolate parametric distribution functions is: 1. Fit a parametric distribution (e.g., a gamma or exponential distribution) at each sampling site to the empirical distribution function (EDF).
2. Interpolate the moment(s) or parameter(s) of the fitted parametric distribution.
3. Set up the theoretical cumulative distribution function (CDF) at every interpolation target with the interpolated moment(s) or parameter(s).
The newly proposed procedure for non-parametric distribution functions is: 1. Fit a non-parametric distribution to log-transformed rainfall values using a Gaussian kernel.
2. Estimate the interpolation (kriging) weights with the precipitation values of a certain quantile.
3. Apply these weights to the values of certain discrete quantiles.4. Linearly interpolate the remaining quantile values to receive a continuous CDF for all target locations.
In Arns et al. (2013) a similar approach is used to interpolate quantile value differences of water levels for a bias correction between empirical distributions of observed and modeled values at the German North Sea Coast.In contrast to their work, entire theoretical distribution functions through interpolation are estimated in our work.Goulard and Voltz (1993) introduced a curve kriging procedure to regionalize fitted functions, which was further developed by Giraldo et al. (2011).Based on their work Menafoglio et al. (2013) developed a universal kriging approach for the non-stationary interpolation of functional data, which was applied in Menafoglio et al. (2016) for the simulation of soil particle distribution functions.As CDF curves are special functions, which are monotonically non-decreasing between 0 and 1, the curve kriging procedure additionally needs to be constrained to these conditions.Our approach can deal with these conditions directly.
After describing the study region Baden-Württemberg in section 2, the concept of precipitation amount models is introduced in section 3. The data selection in section 4 is followed by an investigation of the spatial dependence of the precipitation amount models in section 5.The theory of precipitation amount models is addressed in section 6 and the basis of the proposed interpolation procedure for non-parametric models is established in section 7. The application of different regionalization procedures for precipitation amount models is explained in section 8.The implementation of daily rainfall observations within the interpolation of sub-daily distribution functions is outlined in section 9.The resulting performance of different precipitation amount models at point locations and their regionalization is depicted in section 10.

Study region and data
The study region is the federal state of Baden-Württemberg, which is located in the southwest of Germany.The mountain range Black Forest in the western part and the mountain range Swabian Alps extending from southwest to northeast exhibit the highest elevations in Baden-Württemberg.The rising of large scale moist air masses across mountainous regions causes higher rainfall amounts on the windward side and lower amounts on the leeward side.In the summer months, slopes with differing inclinations lead to a warming of the air triggering convection currents, leading to a greater number of showers and thunderstorms over mountainous regions.This shows a dependence of rainfall on elevation with seasonal differences.The rain-bearing westerly winds lead to high rainfall amounts in the Black Forest.The relatively lower altitude of the Swabian Alps results in lower rainfall amounts as they lie in the shadow of the Black Forest (Landesanstalt für Umwelt, Messungen und Naturschutz Baden-Württemberg (LUBW), 2006).
The years from 1997 to 2011 are chosen as investigation period, as the German Meteorological Service (DWD) set up many new rain gauges in 1997.A relatively homogeneous dataset is obtained by only choosing gauges with observation periods greater than or equal to five years providing rainfall measurements for at least 80 % of the time steps within their observation period.It turned out that we had access to (i) 242 hourly and 5 min resolution and (ii) 347 daily gauges available in the study region, with 80 sites having both high and daily resolution instruments.The observations are provided by the DWD and the Environmental Agency of Baden-Württemberg (LUBW).The high resolution rain gauges are mostly equipped with tipping buckets and gravimetric measurement devices (Beck, 2013).Fig. 1 shows the study region with the locations of the two sets of rain gauges.
Modeling precipitation amounts in our context means estimating distribution functions.The usage of these distribution functions includes the implicit assumption of temporally independent and identically distributed (i.i.d.) variables.This assumption is generally accepted for daily rainfall as the autocorrelation of consecutive nonzero daily precipitation is relatively small and usually of less importance.For higher temporal resolutions, such as hourly, autocorrelation needs to be incorporated in the model (Wilks and Wilby, 1999).In practice different methods exist to take such a correlation into account.One approach is to include autocorrelation prior to the sampling procedure by using conditional distributions.Conditions may be event statistics like the duration of a rainfall event (e.g., Acreman, 1990) or varying statistical moments depending on the hour of the day (e.g., Katz and Parlange, 1995).Another approach is introducing autocorrelation after the sampling procedure.Bárdossy (1998) uses empirical distributions of hourly rainfall intensities to sample values whose random order is subsequently changed within a Simulated Annealing scheme to consider autocorrelation.In Bárdossy et al. (2000) theoretical representations (CDFs) of the empirical distributions are used to allow for regionalization of the distributions and enable simulations at ungauged locations.
The non-exceedance probabilities of a CDF are referred to as quantiles in this work and their corresponding rainfall values are called quantile values.

Data selection
For applications of rainfall estimates, like hydrological or hydraulic modeling, the correct representation of small rainfall values is not necessary as their contribution to decisive high discharge rates is rather small.Furthermore, tipping bucket gauges lead to wrong estimates especially for low rainfall values (Habib et al., 2001).Relative estimation errors are increasing for decreasing rainfall rates (Nystuen et al., 1996;Ciach, 2003) and they only represent a small part of the total water volume, but the number of smaller rainfall values is rather high.To avoid the negative effect of this high number of inaccurate values and due to their minor importance for further applications, this study focuses on medium and high rainfall values.
Therefore, the quantile threshold (Q th ) for hourly (1H) values is set to 0.95.This means, that values smaller than the quantile value at Q th =0.95 are excluded.To investigate the total water volume represented by rainfall values above this quantile at point locations, the Lorenz-curve (Lorenz, 1905) is used.We considered a water volume analysis for varying quantiles as important, to show that high quantiles not only represent the decisive higher rainfall intensities, but also a large proportion of the total water volume.So focusing on these quantiles during the model setup is likely to lead to a better model, as lower quantiles would disturb the model estimation due to measurement errors and the higher quantiles already represent a great percentage of the total water volume.The volume of the lower quantiles can then be modeled by simple and robust methods as they do not require a very precise estimation due to their high inaccuracy and minor importance.
After arranging the n observations x i in non-decreasing order, the Lorenz-curve L i can be calculated from a population (in our case rainfall values at a single gauge) with the following formula: The hourly threshold quantile values (QV th ) range between 0.2 and 1.6 mm for Q th = 0.95 depending on the location of the gauge (see Table 1).The Lorenz curve in Fig. 2 (a) shows that hourly values above Q th = 0.95 represent between 70 and 95 % of the total water volume (1 -cumulative share of water volume).
Based on hourly values (1H) of the high resolution data set, aggregated rainfall values of different temporal resolutions are obtained: 2 hourly (2H), 3 hourly (3H), 6 hourly (6H) and 12 hourly (12H).Through aggregation of daily values (1D) of the daily data set, 5-daily (5D) and monthly (M) values are obtained.In order to exclude small values and still consider the values producing a high percentage of the water volume, the Q th for sub-daily resolutions are defined with the mean Lorenz-curves in Fig. 2 (b).The mean hourly Lorenz-curve yields 0.15 as cumulative share of water volume for Q th = 0.95 (85 % of the total water volume is represented by larger values), which is also defined as target share for the remaining sub-daily resolutions.
This target share of 0.15 results in the following values of Q th for sub-daily resolutions: 0.93, 0.92, 0.9, 0.86 (see Table 1).For aggregations greater or equal to one day, the number of values is rather small and there estimation errors are lower due to an increasing accumulation time (Ciach, 2003;McMillan et al., 2012).Nevertheless, only values above the highest quantile of 1 mm in the study region are used for the daily (1D) and 5-daily (5D) resolution (see Table 1), as smaller values may still exhibit measurement errors.
For the estimation of basic statistics in Table 1 and for following calculations, rain values of the investigated aggregations smaller than 0.1 mm are set to 0 mm.The reason is to achieve homogenization of the data sets of different years and gauges, as the discretization ranges from 0.01 mm to 0.1 mm depending on the gauge.

Probability distributions of precipitation amounts in a spatial context
This section focuses on the spatial dependence of precipitation amount distributions, as the applied interpolation technique of ordinary kriging (OK) is based on the assumption that the variable of interest (the CDF) is more likely to be dissimilar with increasing distances.For the purpose of describing the development of the distribution functions in space, the test statistic T of the two-sample Cramér-von Mises criterion is used (Anderson, 1962).It evaluates the similarity of two CDFs, in our case the similarity of CDFs from observations of two different point locations.The test statistic T is defined according to Anderson (1962) as: where with N as number of observations of the first sample and M as number of observations of the second sample.Both observations are joined together in one pooled dataset and the ranks are determined in ascending order of all observations in the pooled dataset.r i are the ranks of the N observations of the first sample in the pooled dataset and s j are the sorted ranks of the M observations of the second sample in the pooled dataset.T can be interpreted as the mean difference of CDF values (quantiles) of observed rainfall intensities between both data sets.So, if T increases for increasing distances, the CDFs are less similar for increasing distances.
For the calculations of T , only rainfall values above the different Q th (see Table 1) are used.The graphs in Fig. 3 show a decreasing similarity of the distribution functions with increasing distances over all temporal resolutions, as the values of T are increasing with increasing distances.Note that the average T-values of the hourly (1H) data in Fig. 3 (a) are shown as the highest dashed line in 3 (b).So the continuity of the whole distribution changes in space, not only the continuity of values of a single quantile.This shows the applicability of interpolation techniques like OK.

Precipitation amount models
In the following subsections non-parametric and parametric models for precipitation amounts at single sites are introduced.
Before estimating the non-parametric or parametric distributions at each observation gauge, observations smaller than QV th are censored from the sample of each gauge and QV th is subtracted from the values above them to fit to the support of the theoretical distribution functions [0, ∞).QV th varies from gauge to gauge for different temporal resolutions (see Table 1).
After estimating the theoretical CDFs, the quantiles F are scaled with Q th and QV th is added to the quantile values.Only the monthly resolution is excluded from the whole scaling procedure, as all monthly rainfall values are used.

Non-parametric models
Non-parametric KDEs for precipitation amount distributions were previously used and are described for daily precipitation amounts in Rajagopalan et al. (1997) and Peel and Wilson (2008).By using this non-parametric method no theoretical distribution needs to be preassigned, only a kernel and its bandwidth needs to be chosen.That is why they are assumed to be more flexible.A kernel in this context is a function which is centered over each observation value and is itself a probability density function whose variance is controlled by its bandwidth (Bowman and Azzalini, 1997).The probability density function (PDF) or KDE f (x) of every data set is then constructed through a linear superposition of these kernels (Peel and Wilson, 2008), where n is the number of observed values, K is the kernel function, h is the bandwidth of the kernel, x are discrete kernel supporting points, and x i are observed rainfall values: The estimation of f (x) is performed with an R (R Core Team, 2015) implementation of Wand (2015).However, since our non-parametric interpolation scheme is based on CDFs and not on PDFs, the CDF is needed.In order to obtain a CDF from the KDEs an integration is required, which is done numerically with the composite trapezoidal rule (e.g., Atkinson, 1989).For numerical reasons quantiles slightly greater than 1 are sometimes obtained, which are simply set to 1 to remain in the correct range.
To model right skewed precipitation amounts with their bounded support on [0, ∞), either an asymmetric kernel like the Gamma kernel (Chen, 2000) or a symmetric kernel with a prior logarithmic transformation of the values (Rajagopalan et al., 1997) can be used to avoid boundary bias.A boundary bias occurs when kernels with infinite support are used for data with bounded support, as this would lead to a leakage of probability mass (Rajagopalan et al., 1997).
In this work the symmetric Gaussian kernel with a prior transformation of data to logarithms is chosen, as this is an implicit adaptive kernel method with increasing bandwidths for increasing values and therefore alleviates the need to choose variable bandwidths with skewed data (Lall et al., 1996;Charpentier and Flachaire, 2014).The Gaussian kernel is chosen as it is straightforward and its application is facilitated through several software implementations (Sheather, 2004).The Gaussian kernel K(t) is described in Eq. 6: If the density of the logarithmically transformed observed values y = log(x) is f Y and a Gaussian kernel is used for this density estimation, the density estimation f X of the original values x according to Charpentier and Flachaire ( 2014) is: Finally, the bandwidth h needs to be chosen, which is commonly indicated as the key step for KDEs (e.g., Bowman, 1984;Harrold et al., 2003;Sheather, 2004;Charpentier and Flachaire, 2014) as a poor bandwidth selection may result in a peakedness or an over smoothing of the density estimation.Due to this great importance of the bandwidth selection, the performances of different selection methods are investigated.
1.The simplest and most widely used selection method is Silverman's rule of thumb (Silverman, 1986), which is defined as to obtain the optimal kernel bandwidth h opt,SRT with n sample values, where s is the standard deviation and q 3 − q 1 is the interquartile range.Silverman's rule of thumb (SRT) is deduced from minimizing an approximation of the mean integrated squared error between the estimated and the true densities, where the Gaussian distribution is referred to as the true distribution (Charpentier and Flachaire, 2014).
2. The second method is a plug-in approach developed by Sheather and Jones (1991), which is widely recommended due to its good performance (Jones et al., 1996;Rajagopalan et al., 1997;Sheather, 2004).Instead of using a Gaussian reference distribution it uses a prior non-parametric estimate in the approximation of the mean integrated square error and, therefore, requires numerical calculation (Charpentier and Flachaire, 2014) to find the optimal bandwidth h opt,SJ , which is performed with the R implementation of Wand (2015) within this work.
Instead of minimizing the mean integrated squared error, Bowman (1984) recommended minimizing the integrated squared error through a least squares cross validation (LSCV), which is applied using the R package of Duong (2015).Another common cross validation method is the maximum likelihood cross validation (MLCV).Cross validation methods tend to produce small bandwidths and therefore tend to produce peakedness of the density (Rajagopalan et al., 1997;Sheather, 2004;Peel and Wilson, 2008), which we also observed in our applications.Due to this deficiency both cross validation methods are not considered in what follows.

Parametric Models
Within the parametric procedure four different parametric distributions are used to model the precipitation amounts of all aggregations in this study.Commonly used models are the exponential distribution and the two-parameter gamma distribution (Wilks and Wilby, 1999).The mixed exponential distribution was recommended in Wilks and Wilby (1999) and was firstly used for daily precipitation amounts by Woolhiser and Pegram (1979).In addition to these models the Weibull distribution, which showed good performance for modeling monthly precipitation amounts in Baden-Württemberg (Beck, 2013), is used.
The CDF F (x) and the PDF f (x) of each used parametric distribution are listed in the following.
1.For the exponential distribution with the parameter λ these functions are: 2. For the two-parameter gamma distribution they are: where Γ is the gamma function and γ is the incomplete gamma function.
3. For the two-parameter Weibull distribution F (x) and f (x) are: 4. The mixed exponential distribution exhibits the following functions: Parametric distributions with more than two parameters are not considered, as this would complicate the regionalization of the distributions due to dependencies among the parameters.For the three parameter mixed exponential distribution the parameter α is fixed for the whole study region (Wilks, 2008) transforming it into a two-parameter distribution.
In order to estimate the optimal parameter sets of the presented parametric distributions for each rainfall gauge and temporal resolution the maximum likelihood method (MLM) using numerical maximization via a Simplex algorithm and the method of moments (MOM) are applied.The MLM is applied to all mentioned parametric distributions.In the special case of the mixed exponential distribution the parameter α is varied between 0.01 and 0.5 within the parameter estimation.For each value of α the sum of the log-transformed likelihoods is calculated over all gauges with varying values of the remaining parameters, while the maximum sum defines the parameter set.To apply MOM, the mean x and standard deviation s x of the sample values need to be calculated for the gamma distribution.In order to use MOM for the Weibull distribution, the method described in Cohen (1965) is applied.For the estimation of the mixed exponential distribution parameters MOM is not applied due to its shortcomings described in Rider (1961).MOM is neither applied to the one parameter exponential distribution, as it would yield the same results as those from MLM.
7 Non-parametric distributions in a spatial context In order to establish the basis of the proposed regionalization procedure for non-parametric models and to get a more detailed idea of the spatial relationship of distribution functions, the EDFs of hourly and monthly rainfall intensities of the gauge Stuttgart / Schnarrenberg and its five closest gauges are plotted in Fig. 4. It is therefore not of importance which EDF belongs to which gauge, but rather the relationship that the EDFs have with each other.These two graphs show that the order of the EDFs stays quite persistent over different quantiles for both aggregations, as the EDFs do not cross each other very often.In other words, if one gauge exhibits the highest rainfall values for a certain quantile it also exhibits the highest rainfall values for other quantiles and vice versa.The red and purple EDFs on the left graph illustrate this quite nicely.
A more global look at the spatial relation between different EDFs can be obtained with Spearman's rank correlation ρ xy of quantile values of all gauges for different quantile pairs.As we want to investigate the persistence of EDFs for the whole study region, we are only interested in the ranks or rather the order of different quantile values for differing quantiles, which can be done by calculating ρ xy .In our application the two input datasets for calculating ρ xy represent quantile values of two different pairs of quantiles over all gauges in the study region.These pairwise rank correlations of quantile values of all gauge pairs are calculated starting from Q th until 1 in 0.001 steps for sub-daily aggregations and in 0.005 steps for aggregations greater than or equal to one day.This procedure is repeated until the rank correlation of every quantile with every other quantile is obtained.
Finally the mean values of the rank correlation belonging to each quantile are calculated (see the dotted gray lines in Fig. 5).
The greatest mean rank correlation is indicated with a red cross in this figure, which also defines the control quantile (Q c ) with the highest mean rank correlation.Rank correlations of Q c with the remaining quantiles lead to the dashed lines in Fig. 5.
Fig. 5 demonstrates that most of the rank correlations are greater than 0.85, indicating a persistence of quantile values over a great interval of quantiles as well as over the whole study region for hourly through monthly data.Lower correlations can be observed for the highest and lowest quantiles, which indicates a non persistent behavior for these quantiles.This behavior is similar for all temporal resolutions.Therefore, quantile values of Q c can be used to set up the interpolation weights.Applying these weights to the remaining quantiles from Q th until 1 should lead to good regionalization results for non-parametric CDFs.
In Table 2 the control quantiles Q c with the highest mean correlations are summarized for all temporal resolutions.As the precipitation mechanisms are different in summer and winter in Baden-Württemberg, the rainfall data sets are also analyzed separately for summer (from May to August) and winter (from September to April).Q c is mostly close to the center of the considered quantile ranges, which are also shown in Table 2. Nevertheless, it is worth noting the strong similarity of winter and summer control quantiles Q c .The proposed procedure to interpolate non-parametric distribution functions using the same interpolation weights for different quantiles seems feasible as the persistence of the order for quantile values of spatially distributed rain gauges is evident.Only values of very high and low quantiles show a non-persistent behavior.Therefore, quality measures, which focus on the difference of these values, will be introduced.

Regionalizing of precipitation amount models
In the following, the regionalization of point models in order to obtain precipitation amount models at ungauged locations is described.The used regionalization method OK is introduced first.Then the approaches to regionalize parametric and nonparametric distributions are explained.
As only a short overview of OK will be given, the interested reader is referred to the common geostatistical literature, like Kitanidis (1997), for further information.The empirical variogram γ e (h) is calculated using Eq. 17 where n(h) is the number of gauge pairs for distance h, x i represents the position of gauge i and z(x i ) is the variable value at gauge i.As the distances between rainfall gauges never provide a continuous set of distances, the h in Eq. 17 represents different distance intervals.For following applications the width of the interval of h is 10 km and the maximum distance is 100 km.For the theoretical variogram γ t (h) one single model out of the following four is chosen based on the least squares criterion.The s parameters represent the sills, the r parameters the ranges of the variograms.
1. Gauss model: 2. Spherical model: 3. Exponential model: 4. Matern model (Pardo-Iguzquiza and Chica-Olmo ( 2008), K v is the modified bessel function of second kind): The next step within OK is solving the corresponding equation system to estimate an interpolated value at an unobserved location x 0 : where n is the number of gauges included in the interpolation (10 within this work) and µ is the Lagrange multiplier.
As already outlined in the introduction, either the parameters (Kleiber et al., 2012) or the moments (Haberlandt, 1998;Wilks, 2008) of parametric distributions can be interpolated to regionalize parametric models.Within this work the moments are interpolated, when MOM is used for fitting the parametric distributions.If MLM is used, the parameters are interpolated.
Since only rainfall values above QV th (see Table 1) are used, QV th also needs to be interpolated within the parametric approach.
Kernel smoothed distribution functions do not provide a parameter that can be interpolated, thus a procedure other than that for parametric distributions needs to be applied.By analyzing the spatial relation of rainfall EDFs in section 7, a persistent order of quantile values over a wide range of quantiles is observed.Therefore, the interpolation weights of quantile values for the control quantile Q c (see Table 2) can be applied to the remaining quantiles.
For all gauges the quantile values QV c of the control quantile Q c are estimated with the inverse of the gauge-wise numerically integrated non-parametric CDF F np : With these QV c at the observation points, the interpolation weights φ j for the target locations are estimated with OK (see Eq. 22).Then, these weights are applied to the quantile values of quantiles between Q th and 1 in 0.0001 steps and, finally, the remaining quantile values are linearly interpolated to receive a continuous CDF for all target locations.In order to assure a monotonically increasing CDF only positive interpolation weights are allowed.This makes the use of OK problematic.It can only be used if the equation system (see Eq. 22) is solved with positive weights, which leads to additional constraints: Considering these additional constraints the OK equation system is solved with a SCIPY implementation (Jones et al., 2001) of a FORTRAN algorithm by Lawson and Hanson (1987), which solves the Karush-Kuhn-Tucker conditions for the non-negative least squares problem.In the following, this kriging procedure will be called positive kriging (PK).Another way to solve this extended optimization problem with an application of the Lagrange method is presented in Szidarovszky et al. (1987).The persistence of quantile values described in section 7 also implies the persistence of quantiles.The interpolation of quantiles for discrete rainfall values would therefore also be an option.However, this would complicate the regionalization as not only monotonicity needs to be preserved, but also the value range of quantiles from 0 to 1.

Dependence of sub-daily on daily values
As the high resolution rain gauge monitoring network in the study area is quite sparse and the corresponding time series are often incomplete, it would be useful to include more dense and complete secondary information in the interpolation of the sub-daily distributions.Therefore, the applicability of daily values to improve their interpolation is investigated, as the daily monitoring network has a higher density.A simple disaggregation strategy (rescaled nearest neighbor) of Bárdossy and Pegram (2016) is applied to all days to obtain distributions of sub-daily resolutions at the locations of the daily gauges, allocating subdaily values from the closest high resolution gauge to the daily target gauge.The procedure to incorporate daily values in the interpolation of sub-daily values should be the following: 1. Choose a daily target gauge and allocate sub-daily rainfall values of the closest (concerning horizontal distance) high resolution gauge to it.
2. Aggregate sub-daily values of the high resolution gauge to daily values p sub−daily (t) and calculate a scaling factor for every day t by additionally using the values of the daily target gauge p daily (t) : 3. Multiply all sub-daily values of the nearest gauge with this scaling factor.The scaling factor changes from day to day and simply assures that daily sums of disaggregated sub-daily values at the target gauge equal the daily values measured at the target.
5. Calculate the sub-daily statistic of interest from these scaled values at every daily gauge and incorporate them in the interpolation procedure.
The applicability of this procedure is tested with a cross validation, which is described in section S3 of the supplement.For the incorporation of daily values within the regionalization of parametric and non-parametric sub-daily distributions a special regionalization technique is not needed.The rescaling method (NNS) is applied to all available daily gauges.If for a certain day no hourly values are available for the closest gauge, the next closest gauge is used for the rescaling of that certain day in order to increase the sub-daily sample size at the daily gauge.After obtaining the sub-daily values at the daily gauges, they are simply treated as additional control points for the regionalization.

Performance
This section is divided into three parts.In 10.1 the quality measures are introduced, in 10.2 the performance of the precipitation amount models for point wise estimations are compared for all temporal resolutions.The regionalization of the precipitation amount models is addressed in 10.3.The precipitation amount models are fitted and regionalized separately for winter (from September to April) and summer (from May to August) months, as the rain-producing weather processes are different in these two seasons.

Quality measures
The validation of the precipitation amount models at point locations and their regionalization is evaluated with two different quality measures.These quality measures need to be measures considering the CDF and not the PDF, as the interpolation of the non-parametric distributions only provides CDFs for ungauged locations.
The most common goodness of fit test to estimate the quality of fitted distributions is the Kolmogorov-Smirnov test.As distributions of precipitation amounts are positively skewed, the greatest part of the values are small or medium values, which leads to the highest gradient of the CDF for these values.Therefore, a greater difference of the corresponding CDF quantiles would be more likely and would govern the Kolmogorov-Smirnov test.However, these medium values are less important than the greater precipitation amounts for most of the precipitation model applications.
For this reason the Cramér-von Mises criterion as a more integral measure and a Lorenz-curve based measure -which allows for conclusions about the representation of the water volume -are used.The Cramér-von Mises criterion W 2 for single samples is (Stephens, 1974): where F (x i ) represents the theoretical distribution (non-parametric or parametric) of the observed values x i in ascending order.
For sub monthly resolutions the Cramér-von Mises criterion is slightly modified, as only quantiles above Q th (see Table 1) are used: As already mentioned in section 7, a quality measure, which describes the representation of high quantiles, is needed.For Lorenz-curves, high vertical differences are supposed to appear more frequently for high quantiles as the slope increases with increasing quantiles.Therefore, a measure respecting the vertical differences of the Lorenz-curves is suitable.In section 4 the estimation of the Lorenz-curve with observed rainfall values was described.However, the Lorenz-curve L(F (x)) can also be estimated from the theoretical CDF F (x), which is a preferable approach, as random rainfall values do not need to be generated from the CDF previous to the Lorenz-curve estimation: where x(F ) is the gauge wise quantile function (the inverse of the CDF).The integrals of the quantile functions are estimated numerically, because the non-parametrically estimated distribution functions are not invertible analytically.The Lorenz-curve criterion L d used here is the squared difference of the observed L(F n (x)) and modeled Lorenz-curve L(F (x)): The differences of the Lorenz-curves are only estimated for values greater than QV th (see Table 1).Within the validation of the regionalization only values above the highest QV th among the observed and regionalized values for each gauge are evaluated, as they may differ for the different techniques.

Point models
To determine an overall performance ranking for the remaining models, at first the arithmetic mean and the median over the number of gauges of both measures of quality -the Cramér-von Mises criterion W 2 and the Lorenz-curve criterion L d -are calculated for each precipitation amount model.This leads to four different measures, which are shown for hourly values of the winter season in Table 3.Note that the mean values reflect the robustness and the median values represent a good average performance of one precipitation model for the whole study region.
To combine the four statistics (mean and median of W 2 and L d respectively) in one single performance measure, every value in Table 3 is then divided by the smallest (best) value (bold numbers) of its corresponding quality measure, indicating the relative performance with respect to the best model.This leads to one number for each statistic and precipitation model starting from 1 for the best performing model of each statistic.The bigger this number, the worse its relative performance.These four numbers are then combined by adding them together, which results in a single number for each precipitation amount model to define the performance ranking for each temporal resolution.A ranking number of 4 is the lowest possible number and implies that the related model shows the best performance for all four quality measures.In Table 4 the ranking numbers for all temporal resolutions and both seasons are shown.
With the ranking numbers the best performing precipitation amount model is estimated for each season and temporal resolution.Among the non-parametric methods (NP) Silverman's rule of thumb (SRT) and the plug-in approach of Sheather and of Beck (2013) for the same study region.
The performance ranking of the different methods is quite similar in winter and summer.The non-parametric methods always lead to better performances concerning the Cramér-von Mises criterion W 2 .The parametric estimations mostly lead to better results regarding the Lorenz-curve criterion L d (for details see Table S2 and S3 in the supplement).Fig. 6 may provide an explanation for the differences in performance regarding these two quality measures.The graphs show the CDFs and Lorenzcurves for the hourly (1H) and 12 hourly (12H) resolution for a chosen gauge.For the hourly resolution the non-parametric SRT method leads to better results for both measures.An equally good performance regarding the W 2 for the parametric and non-parametric method can be observed for the 12 hourly resolution.However, the non-parametric method performs worse regarding the L d measure, as it overestimates the water volume represented by the higher quantiles.The reason can already be observed in the CDF, where the non-parametric method systematically overestimates the values of high quantiles.The parametric method can lead to over-and underestimations.This influences the W 2 criterion in the same way as a constant overestimation (see squared differences in Eq. 26), but it seems to lead to better results regarding the L d criterion.
Parameter estimation through MOM in combination with the Weibull distribution performs better for higher aggregations, which exhibit more symmetric distributions.For daily and sub-daily aggregations MLM parameter estimation in combination with the mixed exponential distribution leads to better results.
The overall performance is best with the mixed exponential distribution for temporal resolutions between two hours (2H) and one day (1D) in both seasons.For five daily (5D) resolutions the Weibull distribution exhibits the best overall performance in both seasons.For the hourly distribution (1H) the non-parametric models show the best overall performance in both seasons.
Only for the monthly distribution (M) the best performing methods differ between the two seasons.In the summer season the Weibull distribution shows the best results and in the winter season the non-parametric models perform the best.

Regionalization
In order to estimate the quality of the regionalized precipitation amount models, a 2-fold cross validation (split sampling) is used.Two equally sized samples of observation points are randomly generated (Fig. 7).The most simple regionalization method is using the estimates of the nearest neighbor (NN) of the calibration set, which are therefore used as benchmarks for the quality of the regionalization procedure.Additionally, the daily rescaled nearest neighbors (NNS) are used as a benchmark.
In this case all daily gauges are used for the rescaling except for the daily observations at the locations of the respective validation sample.
Following the results of the point-wise estimation in the previous section only the Weibull-MOM and the Mixed-Exp-MLM models among the parametric models are investigated for the regionalization, as they show good performance for differing aggregations.They are both investigated for all aggregations to test the difference of interpolating moments or parameters, except for the monthly aggregation, for which only the Weibull distribution is investigated.In order to regionalize the Weibull-MOM model, the mean and standard deviation are spatially interpolated, and for the regionalization of the Mixed-Exp-MLM model, the parameters λ 1 and λ 2 are interpolated while its parameter α is kept constant for the whole study region.
As the two non-parametric approaches SRT and SJ show very similar results during the point wise estimation, only the SRT approach is interpolated.For the regionalization of the non-parametric models QV c (see Table 2 and Eq.23) values are used to estimate the interpolation weights, which are further applied to the remaining quantiles.Following the conclusions in section 9, daily gauges can be used to set up distribution functions for sub-daily values with a scaled nearest neighbor approach (NNS).

Variogram estimation
The first step during the regionalization procedure is the estimation of the theoretical variograms.The interpolation variables of the three precipitation amount models for which theoretical variograms need to be estimated for the two seasons and eight temporal resolutions are: 1. P-Mixed-Exp-MLM: λ 1 , λ 1 2. P-Weibull-MOM: mean, standard deviation 3. NP-SRT: QV c values (see Table 2 and Eq.23) During the estimation of the parameters of the Weibull distribution with MOM, QV th is subtracted from the rainfall values prior to the estimation of the mean and the standard deviation.As the mean of these values show lower spatial dependencies than the mean of the censored values without subtraction, QV th is added to the mean values of the parameter estimation before the regionalization.After the regionalization, they are subtracted again to determine the parameters of the Weibull distribution.
Variogram models are also fitted to QV th , as the corresponding values serve as starting points for the parametric models at the ungauged locations.Fig. S4 to S7 in the supplement show exemplary theoretical variograms of different parameters for temporal resolutions of 1H and 12H for the winter and summer season of calibration sample 2.
It is difficult to compare the spatial persistence of T (see Fig. 3) with the spatial persistence of the different distribution parameters, as T considers the whole distribution function and the distribution parameters only describe properties of the distribution.However, the range of T was about 35 km, which can also be observed for some of the parameters, especially the mean of P-Weibull-MOM, QV c of NP-SRT and QV th .

Precipitation amount models
The regionalization of the precipitation amount models is evaluated with the same quality measures as the point wise estimation, the Cramér-von Mises criterion W 2 and the Lorenz-curve criterion L d .The investigated interpolation approaches for the parametric distributions are: 1. OK -MOM: OK of the Weibull distribution, fitted with MOM.
2. OK -MLM: OK of the mixed exponential distribution, fitted with MLM.The interpolation approaches for the non-parametric models are:

OK -MOM
1. PK -NP: PK of the non-parametric models, which are estimated using SRT.
2. PK -NP Daily: PK of the non-parametric models including scaled NNS values of daily gauges (only for sub-daily aggregations).
In Fig. 8, parts of the interpolation procedure for PK -NP are shown for the daily aggregation, where the non-parametric QV c at the calibration gauges and three interpolation fields are shown.
In Table 5 and Table 6 the performance ranking numbers of the regionalized precipitation amount models are summarized for the winter season and for the summer season respectively.The differences between the two cross validation samples are quite small, so the performances are not just resulting from the positioning of the gauges in the samples but from the interpolation approaches.Among the parametric methods the MOM approaches mostly perform better than the MLM approaches for aggregations greater than or equal to 2H during the winter season.In the summer season the MOM approaches perform mostly worse than the MLM approaches for aggregations smaller than 6H and vice versa for higher aggregations.Interpolating moments, therefore, seems to be more robust than interpolating parameters of distributions as the performance ranking changed in favor of the MOM approaches compared to the point wise results (see Table 4).Only for stronger skewed distributions in the summer and smaller aggregations, the MLM approach still outperforms the MOM approach.
Comparing the non-parametric interpolation approaches with the parametric interpolation approaches shows that the nonparametric approach performs best for hourly (1H) and two hourly values (2H) for both calibration samples and for calibration sample 1 with three hourly values (3H) in the winter season.This seems to indicate a more robust non-parametric interpolation method for the winter season as the performance ranking changed in favor of it in comparison with the point wise estimation.
In the summer season the non-parametric methods only perform best for the hourly resolution (1H), which is similar to the results of the point wise estimation.
It is obvious that using scaled values of the daily gauges is very beneficial as approaches incorporating these values almost always include the best performing method, except for the 12H aggregation in the summer season.
As a benchmark, the interpolation results are also shown for parametric and non-parametric estimates of nearest neighbors (NN) and additionally using scaled daily gauges for sub-daily aggregations (NNS).Among the benchmark methods the NNS approaches perform better than the simpler NN approaches for the sub-daily aggregations, except for the twelve hourly (12H) resolution in summer.Since the best interpolation approach almost always -with only two exceptions -performs better than the best nearest neighbor approach, the regionalization of distributions seems to be worthwhile.

Conclusions
Comparing different modeling schemes for precipitation amounts at point locations (see Table 4) over different temporal resolutions has revealed several findings.The non-parametric estimates only perform better for the hourly resolution in both seasons and for monthly distributions in the winter season.They have problems especially in reproducing the volume correctly, as they seem to have difficulties with high quantiles.Causes for this deficiency could be the numeric interpolation or the small number of rainfall values at high quantiles.For temporal resolutions between two hours and a month the parametric distributions outperform the non-parametric distributions.Among the parametric methods the mixed exponential distribution performs better for sub-daily and daily aggregations, whereas the Weibull distribution has the advantage for higher aggregations.
The regionalization of the precipitation amount models showed (see Table 5 and 6) that the proposed interpolation scheme for non-parametric distributions is useful as it does not worsen its performance ranking compared to the estimation at point locations.Rather, it appears to be a robust interpolation scheme as it more often outperforms the parametric schemes comparing point wise estimation and regionalization.Among the parametric methods the interpolation of moments turned out to be more robust than the interpolation of parameters.
As auxiliary variables the use of daily gauges for sub-daily resolutions is very beneficial, as was suggested by our data analysis in section S3 in the supplement and is also proven by the evaluation of the regionalization.
In general, the regionalization of distributions seems to be worthwhile as it nearly always performs better than the nearest neighbor (horizontal distance) approaches, which would be the most simple estimate.As lower rainfall values were excluded in this study due to their minor importance and measurement errors, the results are not directly comparable to those of most of the other publications within this research field.
The difficulty of non-parametric distributions in representing water volumes may be reduced by using the Epanechnikov kernel with finite support as proposed by Rajagopalan et al. (1997).Additionally, ways of incorporating elevation within the regionalization of non-parametric distributions need to be tested.Regarding the parametric distributions, Chen and Brissette (2014) and Li et al. (2012) recommended Pareto type distributions instead of exponential type distributions, which could also be tested in the future.Finally, the non-parametric interpolation approach could also be applied to parametric or empirical distributions and should be tested for various study regions.
12 Competing interests The authors declare that they have no conflict of interest.

Data availability
The used sub-daily precipitation data sets were obtained from the LUBW during various research projects and are not available to the public as far as the authors know.Therefore, they can not be provided by the authors.The daily data set was downloaded from the Webwerdis homepage (http://www.dwd.de/DE/leistungen/webwerdis/webwerdis.html) of the DWD, where a personal account is required.However, to the knowledge of the authors any academic researcher can apply for a personal account and some of the used daily and sub-daily values also seem to be available on the homepage of the DWD Climate Data Center (http://www.dwd.de/DE/klimaumwelt/cdc/cdc_node.html).
Jones (1991) (SJ) show very similar results, especially in the winter season.The mixed exponential distribution with a MLM parameter estimation (Mixed-Exp-MLM) leads to the best results among the parametric methods (P) for daily and sub-daily resolutions.For temporal resolutions greater than 1D the Weibull distribution with a MOM parameter estimation (Weibull-MOM) leads to the best results.The best performance of the Weibull distribution for monthly values coincides with the results Daily: OK of the Weibull distribution including scaled NNS values of daily gauges (only for sub-daily aggregations).4. OK -MLM Daily: OK of the mixed exponential distribution including scaled NNS values of daily gauges (only for sub-daily aggregations).

Figure 1 .
Figure 1.Locations of high resolution (hourly and 5 min, left) and daily rain gauges (right) in Baden-Württemberg.

Figure 2 .
Figure 2. In (a) the range of the Lorenz-curves and the mean Lorenz-curve for hourly rainfall values of all rainfall gauges inside the study region are shown, in (b) the mean Lorenz-curves are shown for different temporal resolutions.

Figure 3 .
Figure 3. T statistic over distance: (a) shows the results for hourly distribution functions of all gauge pairs (grey crosses) and their mean calculated for 5 km classes.(b) shows mean values of the T statistic for different temporal resolutions (for more detail on the temporal resolutions of 1D, 5D and M see Fig. S2 in the supplement).

Figure 4 .
Figure 4. EDFs of hourly (a) and monthly (b) precipitation amounts for the gauge Stuttgart / Schnarrenberg and its five closest gauges for a quantile interval.It shows that the order of the EDFs is quite persistent over a wide quantile range for low and high resolutions.Note: As the daily and hourly data set are not the same, the colors in the two graphs do not correspond to the same gauges.

Figure 6 .
Figure 6.Exemplary empirical (data), non-parametric (SRT) and parametric (Mixed-Exp) CDF (left) and Lorenz-curve (right) for hourly (1H) and 12 hourly (12H) resolution of a chosen gauge.Also the values of the two quality measures L d and W 2 are indicated.

Figure 7 .
Figure 7. Locations of the two 2-fold cross validation samples for sub-daily (left) and daily gauges (right).

Figure 8 .
Figure 8. Illustrations for the kriging procedure of non-parametric distributions with daily values (1D) of the summer season using calibration sample 1 (see Fig. 7): In (a) the non-parametric QVc of Qc = 0.865 at the gauges are shown, which then lead to the interpolated values in (b) using interpolation weights φj resulting from PK.The same interpolation weights φj are used for the remaining quantiles, for which exemplary results are shown in (c) for the quantile = 0.72 and in (d) for the quantile = 0.98.An exponential variogram with a range of 41 km and a sill of 2.2 mm 2 is used.

Table 1 .
Basic rainfall information of the study region for different aggregations (agg): P0 is the probability of 0 mm rainfall, Q th stands for the defined quantile thresholds or threshold ranges and QV th represents the corresponding quantile values (rainfall) for the defined Q th .

Table 2 .
Control quantiles (Qc) which exhibit the highest mean pair wise rank correlations with other quantiles.They are shown for different

Table 3 .
Mean and median of the two quality measures W 2 and L d for the eight precipitation amount models over the study region for hourly values (1H) in the winter season.The bold numbers indicate the lowest (best) value of the corresponding measure.

Table 4 .
Performance ranking numbers of the precipitation amount models for point wise estimations.The underlined numbers indicate the best parametric (P) and non-parametric (NP) models.The bold numbers indicate the overall best model.

Table 5 .
Performance ranking numbers for the 2-fold cross validation of regionalized precipitation amount models in the winter season.The underlined numbers indicate the best parametric (P) and non-parametric (NP) models.The bold numbers indicate the best overall model for each validation sample and temporal resolution.

Table 6 .
Performance ranking numbers for the 2-fold cross validation of regionalized precipitation amount models in the summer season.The underlined numbers indicate the best parametric (P) and non-parametric (NP) models.The bold numbers indicate the best overall model for each validation sample and temporal resolution.