Interactive comment on “ Multiplicative cascade models for fine spatial downscaling of rainfall : parameterization with rain gauge data ” by D . E

text, the approach based on the moment-scale diagrams is deemed to be no suitable, as the gauge density is low relative to the desirable grid cell density. Therefore, the Authors exploit the properties of the -stable distributions to derive a set of analytical expressions that allows estimating the model parameters, resorting to the ratio between the rainfall values at the two limit scales of the cascade without using the intermediate scales. The proposed method is applied to calibrate eight versions of a beta-logstable model, in which the scale parameter of the distribution of the weight generator is assumed to be constant or varying with the scale and/or the large scale rainfall, whereas the intermittency parameter depends on the large scale rainfall. These models and the estimation method are tested on a rainfall data set recorded by 25 gauges located around Warsaw, Poland (the data are aggregated at 15-min temporal resolution). The performance is assessed by comparing the cumulative distribution functions, and the semi-variograms of the observed and simulated rainfall values. The approach seems to be very interesting, hypotheses and possible weaknesses are clearly stated, and, in general, the paper is clear and well organized. In my opinion, the paper is suitable for publication in HESS after some minor changes. Some technical comments are provided in the next sections. Please, note that some of my remarks should be considered as a contribution to the open discussion rather than requests of changes.


Introduction
Urban catchments, due to their diminished damping properties relative to rural and natural catchments, are particularly responsive to bursts of local, high intensity rainfall.This makes characterization of the spatial distribution of rainfall at small time scales critical to evaluating the efficacy of urban stormwater drainage systems.Traditionally, design storms have been used to evaluate these systems in conjunction with rainfall-runoff and hydrodynamic models, but in recent years there has been a push towards stochastically downscaling long (e.g., multi-decadal) time series of coarse (e.g., daily) rainfall to higher resolution (e.g., minutes) with which to force models of stormwater drainage systems (e.g., Hingray and Ben Haha, 2005;Molnar and Burlando, 2005;Licznar et al., 2011a).Advantages of using long time series are that they allow for a statistical analysis of system performance and they eliminate the problem of defining the appropriate initial catchment water storage for a design storm (Hingray and Ben Haha, 2005).Furthermore, long time series of daily rainfall are already abundant and readily available, and time series of high-resolution rainfall with which to develop downscaling models are becoming more prevalent.
While using long time series station data provides advantages, their remains the issue that the rainfall field is continuously evolving through time.While one might simplify the problem by using a predefined and static dimensionless rainfall field, this takes away a key strength of the stochastic simulation approach.An alternative is to stochastically downscale the rainfall field as well as the time series.For this purpose, a number of models for stochastically downscaling rainfall fields have been developed.Following Ferraris et al. (2003), most can be grouped into three general types: autoregressive models, point-process models that randomly position rainfall "cells", and fractal and multifractal cascade models.Additionally, there are hybrid models that combine features of these different approaches.For an overview of these various types, see Ferraris et al. (2003) and references therein.We focus on multifractal cascade models because, as noted by Veneziano et al. (2006), multifractal models are simpler and have fewer parameters, and furthermore, though we do not consider these properties in this study, one can deduce the frequency distribution of rainfall intensities and rainfall extremes from their multi-fractal structure.
Parameter estimation for spatial downscaling models requires observations of the rainfall field.With multifractal cascade models, parameter estimation has mostly been done using radar-derived rainfall fields, though in a small number of cases rainfall fields were generated by interpolating rain gauge data (Svensson et al., 1996;Jothityangkoon et al., 2000;Sharma et al., 2007).However, when the gauge density is coarse relative to the final spatial resolution of interest, the interpolation methods will fail because they smooth out the fine-scale variability.
It is common for large metropolitan areas to have in excess of twenty rain gauges installed, whereas reliable fine scale radar-rainfall is less common (e.g., Thames Water, 2010).Even where radar imagery is available, there is value in estimating model parameters directly from rain gauge data, given that accurate rainfall estimation from radar is complex and continues to be a focus of research (Krajewski and Smith, 2002;Pepler et al., 2011).Therefore, we propose a simple method of calibrating a multifractal cascade model for generating rainfall fields of short-duration rainfall (e.g., 15 min) when information across the full field is not available, or specifically, when only data from a network of rain gauges is available.The number of gauges needs only be sufficient to adequately estimate the variance in the ratio of the rain rate at the rain gauges to the areal average rain rate across the entire spatial domain.We expect the particular number will depend on the degree of spatial variability across the domain of interest.We apply the calibration method to precipitation over Warsaw, Poland, and discuss error and bias in the estimation of the model parameters.
In this study, we do not consider temporal evolution of the rainfall fields, which is required for a complete spacetime downscaling model.Various cascade-based space-time models based on cascades have been proposed (e.g.,Over and Gupta, 1996;Venugopal et al., 1999;Deidda, 2000;Jothityangkoon et al., 2000;Kang and Ramirez, 2010).Parameterization of a space-time model will be a topic of a subsequent paper.

Data
Rainfall data were collected from a network of 25 rain gauges distributed throughout Warsaw, Poland.The gauges were installed by Warsaw Waterworks in the fall of 2008 to better characterize storm systems with the specific objective of modeling combined sewer-stormwater systems.Individual gauges were located to obtain best representative meteorological observations in urban settings (Oke, 2006) and to have approximately constant gauge density over the entire city (Fig. 1).The gauges were connected to a single data acquisition system by means of general packet radio service (GPRS) modems.The data used in this study were recorded with a temporal resolution of 1 min and cover the period from the 38th week of year 2008 up to the 49th week of year 2010.For our analysis, data were included only when 21 or more gauges were operating.
All gauges were weighing-type instruments suitable for both liquid and solid precipitation (MPS systém Ltd.,model TRwS 200E).The manufacturer's claimed accuracy was 0.1 % and the resolution was 0.001 mm.Field tests of the installed gauges were conducted prior to operational use.Good agreement between total depth of known and recorded precipitation was observed.However, at a 1-min resolution, the output signal was detectably more damped and broader than the input signal.As a consequence, rain was at times still being recorded for up to a few minutes after water was no longer being added to the gauge funnel.To reduce the relative error caused by this modulation of the signal, we aggregated the data to 15-min intervals.
Precipitation occurs in Warsaw as rain and snow and is generated during both frontal and convective storms.Scaling statistics may vary by precipitation type and storm type, so categorizing data by distinct meteorological processes can be revealing (Lovejoy and Schertzer, 1991;Harris et al., 1997).In Warsaw, snow is limited mainly to the months of November through April, and averages 61 % of all precipitation in February (De ¸bski, 1959).Though data on Warsaw storm types were not available to us, a recent analysis of precipitation and circulation patterns at Kraków, Poland, 268 km south of Warsaw, showed daily precipitation events in the summer were nearly evenly divided between frontal rainfall and non-frontal rainfall (Twardosz et al., 2011).In winter, non-frontal rainfall was half as frequent as frontal rainfall, while non-frontal snowfall was 50 % more frequent than frontal snowfall.The implication is that simply dividing the dataset by season would not be adequate, and we do not have sufficient information to categorize the Warsaw data by meteorological process nor by precipitation type.Given that the focus of this paper is not the precise characterization of Warsaw precipitation, grouping all the data does not detract from our primary purpose.From hereon, we make no distinction between rain and snow (as rain equivalent) and refer to all precipitation as "rainfall" to be consistent with the existing modeling literature.

Spatial downscaling model
Our downscaling model is based on a discrete multiplicative random cascade (MRC).In the discrete MRC model of rainfall fields, the small-scale rainfall rate per unit area in a square cell at the nth cascade level is given by where the area of n is given by L 2 0 b −n .Here the largescale rainfall rate R 0 is the rainfall amount over some interval of time per unit area over the host cell with area L 2 0 .The constant b is the branching number, or number of sub-cells (in our case, 4) into which rainfall from a cell is partitioned at the next level in the cascade (Fig. 2).For each level, the index pair (j , k) represents the cell along the path to the nth level cell.The cells at the n-th cascade level are indexed by n,k , k = 1, 2, . . ., 4 n (see Over and Gupta, 1996).The cascade weight W is a random variable with a prescribed distribution function, of which various types have been proposed in the context of rainfall (e.g., Schertzer and Lovejoy, 1987;Gupta  and Waymire, 1993;Over and Gupta, 1996;Deidda et al., 1999;Ahrens, 2003).
The weight W is generated as a random quantity with the following probability density: where P denotes probability, p is a parameter and W + are the non-zero (positive) weights (Over and Gupta, 1994).Equations ( 2) and (3) comprise the cascade generator: Eq. ( 2) generates the intermittency in the rainfall field (subareas of zero rainfall), while Eq.(3) generates the rainfall volumes greater than zero.The non-zero weights W + have a log-stable density, which is to say that X = ln(W + ) has a stable distribution with four parameters: the stability index 0 < α ≤ 2, the skewness parameter −1 ≤ β ≤ 1, the scale parameter σ > 0, and the shift parameter −∞ < µ < ∞.We denote the stable distribution by S(α,β,σ,µ).While the shift parameter can be defined in several ways, we follow the definition as given in  Samorodnitsky and Taqqu (1994).Properties of stable distributions in the context of multifractal rainfall fields have been discussed by Schertzer and Lovejoy (1987), Lovejoy and Schertzer (1990), and Gupta andWaymire (1990, 1993), for example.

Parameter estimation
Typically, estimation of spatial cascade model parameters relies on an analysis of the spatial scaling of the statistical moments of the observed rainfall quantities (e.g.,Over and Gupta, 1996;Deidda, 2000;Jothityangkoon et al., 2000;Pathirana and Herath, 2002;Sharma et al., 2007;Kang and Ramirez, 2010).The q-th moment M at each spatial scale λ is calculated as where the spatial scale λ n is given by L n /L 0 .Rainfall rates R n n,k at a particular scale λ n are determined by aggregating observed rainfall into grids with cells of area n .The relationship between the moments and scale is made through log-log plots of M(λ n ,q) versus λ n for various q.Linearity of the individual moments versus scale in log-log space implies either mono-or multifractality.The moment-scaling behavior of a fractal field has the form where τ (q) versus q is either a line (monofractal) or a curve (multifractal).Finally, the parameters of the cascade generator are estimated by fitting a distribution-dependent theoretical function to the empirical relationship τ (q).For examples on how the moment-scaling estimation method would be applied to the MRC model such as the one described in Sect.2.2, see Pathirana and Herath (2002) and Serinaldi (2010).The above estimation method requires observed quantities of R n across a range of scales λ n .Unfortunately, we are hindered by a low gauge density (∼0.25 gauges km −2 ) relative to a desirable grid cell density (on the order of 10 3 cells km −2 , or a resolution of 30 by 30 m).
If we imposed fine resolution grids over our gauge network, very few cells would contain enough gauges to adequately estimate the areal-average rain rate for those cells.
Though we could use deterministic spatial interpolation techniques (e.g., Thiessen polygons) to estimate the rainfall everywhere at every cell in the grid, this would likely result in a much too smooth rainfall surface.
Because of our inability to carry out a reliable analysis of moment scaling in space, we assumed a priori that there is power law scaling of the statistical moments.We based this assumption on previous observations of multifractality in rainfall fields for spatial scales under 30 km (e.g., Kumar and Foufoula-Georgiou, 1993a, b;Perica and Foufoula-Georgiou, 1996;Pathirana and Herath, 2002;Kang and Ramirez, 2010).Multifractality implies that realistic rainfall fields could be reasonably reproduced, in a statistical sense, by a family of parsimonious multiplicative random cascade models.
As an alternative to moment-scaling analysis, we parameterized our model using only the final product of the weights W 1 through W n , which we express by the variable Y as We consider the cases for Y = 0 and Y > 0 separately.From Eq. ( 3), p j is the probability that W j > 0 along a path in any j of n cascade levels.The probability that Y = 0 (which is to say that at least one W j equals zero along the path down all n levels) can be calculated as 1 minus the probability that W j is greater than zero in all n levels.From the binomial distribution function (Ross, 1998) we obtain the solution for the probability that Y = 0: To help us determine the distribution of Y when Y > 0, we defined the variable Y + : where Y + > 0. Noting from Eq. ( 3) that W = W + /p, Y can similarly be defined in terms of W + as Combining Eqs. ( 9) and (10) and using the substitution For constant stability index α, the log-stable distribution parameters for Y + can be easily determined from the logstable parameters of W + because the product of log-stable variables is also log-stable.Let given by W and for α = 1 (Samorodnitsky and Taqqu, 1994).
We estimated parameters for eight variations of varying complexity of the MRC model described in Sect.2.2.The simplest model used the log-normal distribution for W + with parameters that were scale-invariant and rainfallindependent, whereas the most complex used the stable distribution with parameters that depended on both scale and rainfall.Each model is summarized below: 1. SIσ /RIσ /LN: the σ parameter of the cascade generator was scale invariant (SIσ ) and W + was log-normally (LN) distributed and independent of rainfall intensity (RIσ ).
2. SIσ /RIσ /LS: the σ parameter of the cascade generator was scale invariant and W + had a log-stable (LS) distribution and was independent of rainfall intensity.
3. SIσ /RDσ /LN: the σ parameter of the cascade generator was scale invariant and W + was log-normally distributed and dependent on rainfall intensity (RDσ ).
4. SIσ /RDσ /LS: the σ parameter of the cascade generator was scale invariant and W + had a log-stable distribution that was dependent on rainfall intensity.
5. SDσ /RIσ /LN: the σ parameter of the cascade generator was scale dependent (SD) and W + was log-normally distributed and independent of rainfall intensity.
6. SDσ /RIσ /LS: the σ parameter of the cascade generator was scale dependent and W + had a log-stable distribution and was independent of rainfall intensity.
7. SDσ /RDσ /LN: the σ parameter of the cascade generator was scale dependent and W + was log-normally distributed and dependent on rainfall intensity.
8. SDσ /RDσ /LS: the scale parameter σ of the cascade generator was scale dependent and W + had a log-stable distribution and was dependent on rainfall intensity.
For those models where σ was scale invariant, σ was solved for uniquely in terms of σ Z by inverting Eq. ( 12): When σ was scale dependent, σ varied as the following function of the length scale λ = L/L 0 : where γ is a constant and σ 1 is the value of σ at λ = 1.For γ > 0, as the scale λ decreases the variance of W + decreases, which places it in the family of "bounded" cascade models (Marshak et al., 1994).Combining Eqs. ( 12) and ( 15) and using the substitution λ = (1/2) j −1 , σ α 1 can be solved for in terms of σ Z : for γ >0 In all eight models, the intermittency parameter was determined from P (Y > 0): Because it has been observed that spatial cascade parameters, and the intermittency parameter in particular, depend on large-scale rainfall (Over andGupta, 1994, 1996;Deidda, 2000;Jothityangkoon et al., 2000;Pathirana and Herath, 2002;Deidda et al., 2004Deidda et al., , 2006;;Sharma et al., 2007), we allowed some of the parameters of the cascade generator to vary with the large-scale rainfall depth R 0 .While it has been argued that for both space (Veneziano et al., 2006) and time (Veneziano et al., 2006;Rupp et al., 2009;Serinaldi, 2010) the parameters should vary with rainfall intensity at each scale (not just the largest scale), our dataset did not permit us to adequately examine rainfall dependency across scales, therefore we restricted the dependency to the largescale rainfall only.In all eight models, we allowed the intermittency parameter to depend on large-scale rainfall by varying P (Y > 0) in Eq. ( 17) with R 0 as where erf is the error function with parameters m and s (Rupp et al., 2009).In four of the models, the scale parameter σ was varied with rainfall by relating σ Z in Eqs. ( 14) and ( 16) to R 0 as where f () is an arbitrary function and c is a constant.We used cubic splines to determine f ().
Observations of rainfall at gauges were used to estimate values of Y , Y + , and P (Y > 0).Combining Eqs. ( 1) and ( 7) yields Y n = R n /R 0 , from which we see that an estimate of Y from observations of R at a given rain gauge can be calculated as where Ŷ is the estimate of Y , i = 1, 2, . . ., N obs indexes the i-th observation in time, and k = 1, 2, . . ., N gauges indexes the rain gauge.The areal average rainfall R 0,i at the reference length L 0 was approximated by taking the mean of the rainfall measured over all N gauges at time i.To estimate P (Y = 0), we used Pi (Y > 0) = (number of gauges with non − zero rain) i /N gauges ( 21) Finally, Y + was estimated with Because Y + is bounded by zero and positive infinity, whereas the upper limit to Ŷ + is N gauges , the distribution of Ŷ + is only approximately equal to the distribution of Y + .This limitation, plus instrument error at very low and very high rainfall intensities, introduces a bias into the estimation of σ Z .
For now, we simply accept this bias as a shortcoming of the estimation procedure, though we discuss it further in Sect.3. Free software packages for estimating the parameters of the stable distribution are rare, and we found none that suited our particular needs.For this reason, we used a simple procedure to estimate α Z and σ Z from the "observed" values ln Ŷ + .An optimization algorithm minimized the sum of squared differences between the following observed and theoretical quantiles: 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, and 0.95.For the normal distribution, we used the maximum likelihood method.
Note that except for at the largest spatial scale, we did not areally average the precipitation data.We also did not consider the particular location in space of the observed rainfall.Both of these characteristics distinguish our study from others.However, we did use the number of cascade levels n that brought us to the scale of the rain gauge itself.Given that the rain gauges have a diameter of approximately 0.15 m, and that the spatial extent of our rain gauge network corresponds to approximately L 0 = 20 000 m, approximately n = 17 cascade levels are needed.

Model evaluation
Direct comparison of the stochastically downscaled data to the observed rainfall requires disaggregation down to the capture area of the rain gauge through n = 17 cascades levels, which would result in a grid with over 17 × 10 9 cells.This would be impractical, especially given that very many such spatial fields would be generated for time increments of as little as 15 min.Instead of generating complete rainfall fields at the spatial resolution of a rain gauge, we followed the cascade process down a subset of the total number of possible paths down the cascade.Along a path at each subdivision, we randomly chose 1 of the 4 cells (with equal weight given to each cell), and tracked the position (x, y) of the cell for each of the n = 17 cascade levels.In total, we followed 240 paths  per 15-min increment.We used Eq. ( 1) to calculate the 15min rainfall amounts at each of these 240 locations, which served as our representative "gauges".
To further reduce the computational burden, downscaling was done on a sub-sample of all 18 723 15-min time steps where R 0 > 0. We selected 2000 time steps such that the cumulative distribution function (CDF) of the sub-sampled R n was similar to the CDF of the full record.
The spatial structures of the observed and simulated 15-min rainfall were compared using semivariograms of log(R n ) for R n ≥ 0.001 mm (15 min) −1 (the minimum rainfall amount recorded).The spatial structure of the intermittency was examined with semivariograms of presence (1) and absence (0) of rain.
In addition to the spatial structure, we assessed the ability of the models to reproduce the cumulative distribution frequency (CDF) of 15-min rainfall for R n > 0.
To examine the bias in the estimation method, we ran additional simulations as described above in which we recorded the rainfall at 24 locations and used all 18 723 15-min time steps.From both the 240-and 24-gauge datasets, we estimated the model parameters m, s, and α Z and σ Z using the methodology described in Sect.2.3.

Results and discussion
The proportion of gauges with zero rain in a 15-min period, P (Y = 0), was found to be strongly dependent on the largescale rainfall rate R 0 (Fig. 3a).Consistent with many other studies (e.g.,Over andGupta, 1994, 1996;Jothityangkoon et al., 2000;Pathirana and Herath, 2002;Sharma et al., 2007), the sparseness of the rain field was much greater when R 0 was low, while at high rainfall rates the tendency was for it to be raining everywhere.The sigmoidal shape of Eq. ( 18) appears suitable for simulating the rainfall intermittency, as it allows for P = 0) to go to 1 as R 0 goes to 0, and to go to 0 as R 0 goes to +∞.
The empirical histograms of lnY + were rightward skewed, thus more similar to a log-stable density with β = −1 than to a log-normal distribution (Fig. 4).At progressively lower values of R 0 , e.g., <∼0.01 mm (15 min) −1 , the empirical histograms were progressively more dominated by lnY + = 0, such that neither the log-stable nor the log-normal densities matched the observations.It is clear that by fitting theoretical distributions to lnY + at low values of R 0 , we are merely fitting to a data artifact and not to true rainfall behavior.
The value of the stable distribution parameter α Z showed a general increasing trend with increasing R 0 (Fig. 5).However, when α Z was fixed at a constant value of 1.47, the fits of the log-stable distributions were only marginally degraded (Fig. 4).This was fortunate because it allowed us to keep α Z as a constant parameter and only have to vary the scale parameter σ Z with large-scale rainfall.The particular value of α Z = 1.47 is the average α Z using all values of ln Ŷ + for R 0 ≥ 0.004 mm.The threshold of 0.004 mm was selected because below this value the empirical distributions appeared to be strongly influenced by data precision (Fig. 4).
The dependency of σ Z on R 0 was complex (Fig. 5), though cubic splines with no more than 6 knots reproduced the empirical relationship of σ Z with R 0 well.The relationship was similar in form for both the log-stable (α Z = 1.47) and lognormal (α Z = 2) distributions (Fig. 5).The range of σ R 0 are suspect, as we have already determined the empirical distributions to be unreliable.If we ignore σ Z for R 0 < 0.01 mm (15-min) −1 , the pattern in Fig. 5 (lower panel) implies a smoother field at intermediate rainfall rates of about 0.1 mm (15-min) −1 and a more variable field at lower and higher rates.This trend toward higher variability at the highest rainfall rates could be the result of localized, high-intensity rainfall generated from strong convective storm cells.This trend is not evident in the studies of Over and Gupta (1996), Jothityangkoon et al. (2000), Pathirana and Herath (2002), Veneziano et al. (2006), or Sharma et al. (2007), who only observed the scale parameter (i.e., variance) to decrease with increasing R 0 from intermediate to high R 0 .The difference between our and previous results may be due to differences in scales between studies: Jothityangkoon et al. (2000) and Sharma et al. (2007) analyzed daily rainfall, while Over and Gupta (1996), Veneziano et al. (2006) and Pathirana and Herath (2002) analyzed radar scans with resolutions ranging between 1 to 5 km.
The semivariogram of the observed rainfall rates shows covariability in rainfall intensity increasing strongly with increasing separation distance (Fig. 6).In contrast, all the scale-invariant models produced rainfall fields that showed little change in correlation with separation, though proximal rain was slightly more similar than distant rain.When the scale parameter σ as allowed to decrease with decreasing scale via Eq.( 15), however, the general variogram pattern of the observed rainfall could be reproduced for separation distances of less than about 10 km by using a value of γ ∼ 0.8 (Fig. 6).At separation distances above 10 km, the semivariances of the simulated rainfall become nearly constant, irrespective of the model or the value of γ .We believe this is an artifact of the discrete nature of the cascade procedure that was applied, which produces a blocky pattern.Note that at the first cascade level, the rainfall is first separated into four 10 km by 10 km cells.The rainfall simulated at a point in one of these first four cells will be equally correlated with the rainfall simulated at a point anywhere in one of the other three cells.Incorporating non-stationarity into the cascade process might remove this artifact; we return to this point briefly at the end of this section.
The semivariogram of the presence/absence observations shows that if it is raining (or not raining) at one location, it is more likely to be raining (or not raining) nearby than it is further away (Fig. 7).The simulated rainfall fields have this property as well (at least below the 10 km separation distance), but not to the degree of the observed rainfall field.Figure 7 gives semivariograms of the simulated presence/absence data using Models SIσ /RIσ /LN and SDσ /RIσ /LN only: semivariograms from all eight models were similar because the models are identical in how they simulate intermittency.
As mentioned above, the discrete cascade process produces a blocky pattern, which will have some influence of the semivariogram.To generate patterns that are more realistic, a filter may be applied to the discretely generated field (Schertzer and Lovejoy, 1987;Menabde et al., 1997;Watson, R. J. and Hodges, D. D., 2005), or one may opt for a continuous-in-scale cascade, such as the continuousin-scale universal multifractal (UM) model (Schertzer and Lovejoy, 1987;Lovejoy and Schertzer, 2010a, b).How a filter would affect the parameter estimation procedure presented here, and how the parameter estimation would be done in the framework of the UM model, are topics of future study.
Over most of the range of R, the log-normal models reproduced well the observed CDF of rainfall rates (Fig. 8).However, the simulated CDFs using the log-normal models diverged from the observed CDF below 0.02 mm (15-min) −1 .The rainfall-dependent (RD) models performed slightly better than the rainfall-independent (RI) models up to the very highest rainfall intensities.Above about 15 mm (15-min) −1 , the rainfall-dependent models overestimated the rainfall intensity at a given probability of occurrence by roughly a factor of two.The log-stable models did better than the log-normal models at reproducing the overall shape of the observed CDF (Fig. 9), including matching the curvature for R < 0.02 mm (15-min) −1 .The log-stable models also mimicked the upward curvature of the observed CDF for the highest values of R, though they did under-predict the probabilities of these extreme events.Making the scale parameter rainfalldependent resulted in an improved CDF, though at the highest intensities these models still underestimated the rainfall intensity at a given probability of occurrence by as much as a factor of two.
As mentioned previously, the estimation method is biased due to artifacts of the data.For one, the measurement instrument has a detection threshold, which results in sparser measured than true rainfall, and hence an underestimation of p, particularly when the larger-scale rainfall R 0 is low (Veneziano, et al., 2006).A second factor is the recording precision of the rainfall observations, which particularly affects parameter estimation at low rainfall intensities.The effect of precision can be seen in the preponderance of values of Y + = 1 at low values of R 0 (Fig. 4) which, in turn, results in an underestimation of σ .Similar observations were made by Rupp et al. (2009) and Licznar et al. (2011a, b) regarding the empirical weights W + when analyzing rainfall time series.A third factor is the underreporting of high rainfall intensities due to instrument error, which reduces the variance of Y + .Additional bias, as already introduced in Sect.2.3, arises from sampling the full rainfall field with a limited number of gauges.In the case of rain intermittency, sampling introduced error into the estimation of P (Y = 0), as seen from the deviations of the estimated from the assumed values of P (Y = 0) in Fig. 3b.However, bias in the estimation of the intermittency parameters m and s in Eq. ( 18) was very small (Fig. 3b).From the observations, we estimated the pair (m, s) to equal (−3.170, 1.804), while from the simulations using all the scale invariant models, (m, s) averaged (−3.177, 1.793) and (3.211, 1.809) with 240 and 24 gauges, respectively.
The bias effect of sample size was more prominent for the stable distribution parameters.The simulations using the scale-independent and rainfall-independent models provide a good illustration of this effect because the model parameters never varied.In general, when the sample consisted of 240 gauges, the estimation procedure accurately retrieved the assumed values of α Z and σ Z (Fig. 10).However, at progressively lower R 0 , the parameter values were increasingly underestimated, and at the lowest values of R 0 , there were simply to few observations to reliably fit the theoretical distributions to the data.Excluding low R 0 , when the sample consisted of 24 gauges, there was no notable bias in α Z when the rainfall came from a log-stable model, but there was a slight underestimation α Z when the rainfall came from a lognormal model.On consequence is that one might choose a log-stable model when in fact the simpler log-normal is more appropriate.Even so, given the high values of α Z (>1.9) estimated here, use of the more complicated log-stable model would hardly be justified.Again, when using 24 gauges, there was a slight underestimation in σ Z and it would appear that this bias would increase with decreasing sample size.We obtained similar results using the rainfall-dependent models (Fig. 11).
It is clear from Figs. 10 and 11 that estimation accuracy deteriorates at rainfall rates below about 0.01 mm (15min) −1 with 24 gauges and below about 0.001 mm (15min) −1 for 240 gauges.The decrease in σ Z with decreasing R 0 below 0.01 mm (15-min) −1 estimated from the rainfallindependent (RI) simulations means it is possible that the similar decrease in σ Z with decreasing R 0 below 0.01 mm (15-min) −1 from the observations is merely an artifact of the estimation procedure.A variety of procedures could be used to partly account for the bias.One is to iteratively adjust the model parameters until the estimated parameters from the simulated dataset are nearly the same as those from the observed dataset (Veneziano et al., 2006).Another procedure is to exclude some data while estimating parameters.For example, in our study we left out data where R 0 < 0.004 mm (15-min) −1 when estimating α Z and when estimating σ Z for the case where σ Z was assumed to be independent of R 0 .However, this excluded only a relatively small amount of data and thus did not greatly affect the values of α Z and of σ Z independent of R 0 .Furthermore, if the objective were to have rainfall-dependent parameters, excluding the low intensity data would provide no guidance as to which parameter values to actually apply at these low rainfall intensities.In another example of censuring data, Licznar et al. (2011a) simply eliminated what would be analogous in our study to all values of Y + = 1 from the empirical frequency distribution, under the assumption that most of these values were artifactual.A third procedure to deal specifically with recording precision is to add random noise to the rainfall observations, with the intent of replacing the information lost by round-off error and thus removing the discretization that leads to an excess of certain values of W + (or Y + ) (Licznar et al., 2011b).
Bias-correcting procedures such as those above should be explored, and we expect that they would improve the fits of frequency distributions.We know, for example, that both data precision and the finite number of gauges serve to decrease the estimated value of the scale parameter σ Z , in the former case by generating an overabundance of Ŷ + = 1 and in the latter case by imposing a maximum value to Ŷ + of N gauges .A bias-correcting procedure that led to an increase in the value of the log-stable parameter σ Z would produce more extreme events, resulting in a CDF more like the observed one in Fig. 9.It would also increase the semivariance overall, which was generally underpredicted by the log-stable models (Fig. 6, lower panel).
Lastly, we have assumed stationarity in the rainfall field, though there may be long-term spatial patterns across the Warsaw metropolitan area.With our short record length (less than 3 yr) it would be difficult detect any but very clear and strong large-scale patterns, which we did not see.Should continuing observations reveal deterministic patterns in the spatial distribution of rainfall, we could account for these within the MRC framework.Examples of how this might be done using a deterministic field of weights that are applied to the cascade generator are given by Jothityangkoon et al. (2000) and Pathirana and Herath (2002).

Conclusions
We have presented and evaluated a method for estimating the parameters of a multiplicative random cascade model for downscaling rainfall fields when observations of the full fields are not available either from radar imagery or from interpolation of very dense rain gauge network data.The estimation procedure still relies on rain gauge data, but the density of the network need only be such that (1) the rainfall rate over a given time interval averaged over the entire spatial domain can be reasonably approximated by averaging the rainfall rate from all the gauges, (2) the number and the spatial coverage of the gauges are adequate for generating a semivariogram of rainfall intensity.
When the cascade generator is independent and identically distributed (iid) throughout the cascade, the parameters can be estimated solely from the frequency distribution of the ratios of the rain rate at each gauge to the large-scale average Hydrol.Earth Syst.Sci., 16, 671-684, 2012 www.hydrol-earth-syst-sci.net/16/671/2012/ rain rate.We found, however, that an iid cascade generator failed to reproduce the spatial covariance structure of the rainfall over Warsaw, Poland: proximal rainfall was too dissimilar using an iid parameterization, and the simulated rainfall only showed a weak relationship between distance and covariability (or semivariance), whereas this relationship was strong in the observed data.
To better reproduce the spatial structure of actual rainfall fields (as summarized by the semivariogram), we added scale dependence to the cascade generator.The scale-dependent generator introduced an additional model parameter (γ ) that could not be estimated directly from the rain rate ratios.We therefore treated γ as a tuning parameter that was estimated by matching the observed and simulated semivariograms.To keep the model simple (i.e., to one tuning parameter) for this study, we considered only scale dependence in the generation of positive rainfall amounts, not in the generation of rainfall intermittency.A similar strategy, however, could be used for the intermittency parameter along with the semivariogram of rainfall presence/absence, though it would require the introduction of at least one additional parameter.
Overall, the scale-dependent MRC models generated the correct frequency distribution of short-duration rainfall intensities.We recommend, however, further research into bias in parameter estimation; we expect that through biascorrection procedures, improvements could be made at both the extreme lower and upper ends of the distribution.
We evaluated the model by using statistical properties of the rain gauge data as performance targets.This meant it was necessary to downscale to the approximate capture area of the rain gauge (15 × 15 cm).For most stormwater drainage system studies, generating fields at such a fine resolution would be impractical.However, an expedient property of the MRC model is that it lends itself nicely to downscaling to any spatial scale λ n , which can conveniently be used to generate gridded rainfall fields for use as input to hydrologic/hydrodynamic models.

Fig. 1 .
Fig. 1.Location of rain gauges used in study in Warsaw, Poland.

Fig. 3 .
Fig. 3. Probability of zero rain at a rain gauge, or equivalently, P (Y = 0), against large-scale areal-averaged rainfall rate R 0 for R 0 > 0. (a): P (Y = 0) as estimated from the observed rainfall (light gray symbols).(b): P (Y = 0) as estimated from rainfall simulated with Model SIσ /RIσ /LN.The dark gray symbols show the values estimated using 240 "gauges", represented by 240 grid cells in the finest resolution field (approximately 15 × 15 cm).The light gray symbols show the parameter values estimated using 24 such gauges.The colored lines show Eq.(18) fitted to the observed and simulated datasets.

Fig. 4 .
Fig.4.Empirical density (bars) of lnY + and densities of the fitted stable distribution with α Z as a fitting parameter (thin green line), with α Z fixed at 1.47 (heavy blue line) and with α Z fixed at 2 (thin red line).Distributions are shown for various large-scale areal-averaged rainfall rates R 0 (as mm per 15-min) for R 0 > 0. The rainfall value shown in each plot is the midpoint of the range of log-transformed rainfall used to bin that data in a given rainfall rate class.

Fig. 5 .
Fig.5.Estimated stable distribution parameters α Z and σ Z for Z = lnY + against large-scale areal-averaged rainfall rates R 0 for R 0 > 0. The open symbols indicate where estimation was clearly affected by artifacts arising from data precision.In the upper panel, the dashed line shows the estimate of α Z using all data for R 0 ≥ 0.004 mm (15-min) −1 and the dotted line is σ Z = 2 (the normal distribution).In the lower panel, values of σ Z assume constant α Z = 1.47 and α Z = 2 for the stable and normal distributions, respectively.The solid curves are fitted cubic spline functions.The dashed and dotted horizontal lines indicate the values of σ Z averaged over all R 0 ≥ 0.004 mm for the stable and normal distribution, respectively.