Singularity-sensitive gauge-based radar rainfall adjustment methods for urban hydrological applications

Gauge-based radar rainfall adjustment techniques have been widely used to improve the applicability of radar rainfall estimates to large-scale hydrological modelling. However, their use for urban hydrological applications is limited as they were mostly developed based upon Gaussian approximations and therefore tend to smooth off so-called “singularities” (features of a non-Gaussian field) that can be observed in the fine-scale rainfall structure. Overlooking the singularities could be critical, given that their distribution is highly consistent with that of local extreme magnitudes. This deficiency may cause large errors in the subsequent urban hydrological modelling. To address this limitation and improve the applicability of adjustment techniques at urban scales, a method is proposed herein which incorporates a local singularity analysis into existing adjustment techniques and allows the preservation of the singularity structures throughout the adjustment process. In this paper the proposed singularity analysis is incorporated into the Bayesian merging technique and the performance of the resulting singularity-sensitive method is compared with that of the original Bayesian (non singularity-sensitive) technique and the commonly used mean field bias adjustment. This test is conducted using as case study four storm events observed in the Portobello catchment (53 km) (Edinburgh, UK) during 2011 and for which radar estimates, dense rain gauge and sewer flow records, as well as a recently calibrated urban drainage model were available. The results suggest that, in general, the proposed singularity-sensitive method can effectively preserve the non-normality in local rainfall structure, while retaining the ability of the original adjustment techniques to generate nearly unbiased estimates. Moreover, the ability of the singularity-sensitive technique to preserve the non-normality in rainfall estimates often leads to better reproduction of the urban drainage system’s dynamics, particularly of peak runoff flows.


Introduction
Traditionally, urban hydrological applications have relied mainly upon rain gauge data as input.While rain gauges generally provide accurate point rainfall estimates near the 40 ground surface, they cannot properly capture the spatial variability of rainfall, which has a significant impact on the urban hydrological system and thus on the modelling of urban runoff (Gires et al., 2012;Schellart et al., 2012).Thanks to the development of radar technology, weather radar data have 45 been playing an increasingly important role in urban hydrology (Krämer et al., 2007;Liguori et al., 2011).Radars can survey large areas and better capture the spatial variability of the rainfall, thus improving the short-term predictability of rainfall and flooding.However, the accuracy of radar mea-50 surements is in general insufficient, particularly in the case of extreme rainfall magnitudes (Einfalt et al., 2004(Einfalt et al., , 2005)).This is due to the fact that, instead of being a direct measurement, radar rainfall intensity is derived indirectly from measured radar reflectivity.As a result, both radar reflectivity 55 measurements and the reflectivity-intensity conversion process are subject to multiple sources of error.
First, errors in radar reflectivity measurements may arise from blockage of the radar beam, attenuation, ground clutter, anomalous propagation of the signal, among other sources 60 (Collier, 1996;Einfalt et al., 2004;Harrison et al., 2000).Although the radar reflectivity measurements undergo a number of corrections before they are converted into rainfall intensity, it is virtually impossible to have error-free reflectivity measurements.Second, the conversion between radar reflec-65 tivity (Z) and rainfall rate (R) uses the so-called Z − R relationship, Z = aR b (Marshall and Palmer, 1948), where a eral studies, which indicate that rain drop size distribution is a highly dynamic process and may significantly or suddenly change within a storm event (Smith et al., 2009;Ulbrich, 1983).Because of this, the Z − R derived rainfall intensity cannot effectively reflect the short-term dynamics of true rainfall intensities and may statistically compromise to intermediate rainfall intensities.
In order to overcome these drawbacks of radar rainfall estimates while preserving their spatial description of rainfall fields, it is possible to dynamically adjust them using rain 90 gauge measurements.Many studies on this subject have been carried out over the last few years, though most of them focus on hydrological applications at large scales (Anagnostou and Krajewski, 1999;Fulton et al., 1998;Germann et al., 2006;Goudenhoofdt and Delobbe, 2009;Harrison et al., 95 2000Harrison et al., 95 , 2009;;Seo and Smith, 1991;Thorndahl et al., 2014).Only few studies have examined the applicability of these adjustment techniques to urban-scale hydrological applications and concluded that they can effectively reduce rainfall bias, thus leading to improvements in the reproduction of hy-100 drological outputs (Smith et al., 2007;Vieux and Bedient, 2004;Villarini et al., 2010;Wang et al., 2013).However, despite the improvements achieved with the current adjustment techniques, underestimation of storm peaks can still be seen after adjustment and this is particularly evident in the case 105 of small drainage areas (such as those of urban catchments) and extreme rainfall magnitudes (Wang et al., 2012;Ochoa-Rodríguez et al., 2013).This may be due to the fact that the underlying adjustment techniques, mainly based upon Gaussian (1st and/or 2nd statistical moments) approximations (Goudenhoofdt and Delobbe, 2009;Krajewski, 1987;Todini, 2001), cannot properly cope with the so-called "singularities" (which implies non-normality and often corresponds to local extreme magnitudes) observed at small scales (Schertzer et al., 2013;Tchiguirinskaia et al., 2012).In fact, it is often the case that the radar image captures the spatial structure of striking local extremes (albeit the actual rainfall depth/intensity may be inaccurate), but these structures are lost or smoothened throughout the merging process.This deficiency may cause large errors in the subsequent urban hy-120 drological modelling.
To address this limitation and improve the applicability of adjustment techniques at urban scales, a method is proposed herein which incorporates a local singularity (identification) analysis (Cheng et al., 1994;Schertzer and Lovejoy, 1987; 125 Wang et al., 2012) into existing adjustment techniques and enables the preservation of the singularity structures throughout the adjustment process.The singularity-sensitive method is particularly intended to improve geostatistical-based merging techniques (e.g.Cokriging (Krajewski, 1987), Bayesian 130 merging (Todini, 2001), Kriging with External Drift (Wackernagel, 2003) and conditional merging (Sinclair and Pegram, 2005)), which seek to represent the spatial covariance structure of the rainfall field or its errors by making use of the semi-variogram (Goudenhoofdt and Delobbe, 2009).Be-135 ing a 2nd-order tool, the semi-variogram cannot adequately capture higher-order features of the rainfall field, thus causing these to be lost in the merging process.
The proposed singularity-sensitive method was initially developed and preliminarily tested in the reconstruction of 140 a storm event which led to reported flooding in the Maida Vale area, Central London, in June 2009 (Wang and Onof, 2013;Wang et al., 2014).The radar rainfall product for this event showed strong and localised singularity structures, but the accuracy of the actual estimates was poor.A dynamic 145 gauge-based adjustment was conducted using the Bayesian data merging method (Todini, 2001), which in previous studies had been shown to outperform other adjustment methods (Mazzetti and Todini, 2004;Wang et al., 2013).Nonetheless, for this particular event records from only a few rain 150 gauge sites were available and these were located away from the area of interest and at points where less intense radar rainfall was observed.Under these circumstances, the aforementioned shortcomings associated with existing adjustment techniques became evident.The Bayesian data 155 merging method proved inadequate as it smoothed out the singularity structures, which had the effect of considerably reducing the peak rainfall intensities.It was then that the singularity-sensitive method was devised and effectively incorporated into the Bayesian merging technique.The result-160 ing singularity-sensitive Bayesian merging method led to rainfall fields which better preserved the spatial structure as captured by the radar and better reproduced peak rainfall intensities.
In the present paper the formulation of the proposed 165 singularity-sensitive method is explained in detail and new numerical strategies aimed at improving the use of singularity information are introduced.Moreover, the method is further tested using as case study four storm events observed in the Portobello catchment (53 km 2 ) (Edinburgh, UK) during 170 2011 and for which radar estimates, a spatially dense network of rain and flow gauges, as well as a recently-calibrated urban drainage model were available.The paper is organised as follows.In Sect. 2 a detailed explanation is provided of the theoretical development of plemented numerical strategies.In Sect. 3 we present the case study, including a description of the study area and dataset, the performance criteria used to evaluate the proposed methodology, and the results of the testing.Lastly, in Sect. 4 the main conclusions are presented and future work is discussed.
2 Formulation of the singularity-sensitive bayesian data merging method Firstly, a description is provided of the two key techniques used in this paper: the Bayesian data merging method and the local singularity analysis.Afterwards, the proposed method for integrating these two techniques is explained.Intermediate results of each of the steps described in this section, which help illustrate the main features of the proposed methodology, can be found in the supplementary document "Formulation and demonstration of the singularity-sensitive Bayesian data merging method".

Bayesian radar-rain gauge data merging method
The Bayesian data merging method (BAY) is a dynamic adjustment method (applied independently at each time step) intended for real-time applications (Todini, 2001).The underlying idea of this method is to analyse and quantify the uncertainty of rainfall estimates (in terms of error co-variance) from multiple data sources -in this case, radar and rain gauge sensors -and then combine these estimates in such a way that the overall (estimation) uncertainty is minimised.The BAY merging method consists of the following steps (illustrated in Fig. 1a): a.For each time step t, the point rain gauge (RG) measurements are interpolated into a synthetic rainfall field using the block kriging (BK) technique.The result of this step is an interpolated rain gauge rainfall field, with areal estimates at each radar grid location (y RG t ) and which are accompanied by the associated estimation error co-variance function (V RG εt ), representing the uncertainty of rain gauge estimates.
b.The interpolated rain gauge rainfall field is compared against the radar field (y RD t ), based upon which a field of errors (estimated as the bias at each radar grid location: ) is obtained empirically.Assuming that areal rain gauge estimates are unbiased, the expectation value (µ ε RD t ) and the co-variance function (V RD εt ) of this radar-rain gauge error field at each time step is used to represent, respectively, the mean bias and the uncertainty of radar estimates.c.Using a Kalman filter (Kalman, 1960), the two rainfall fields are optimally combined such that the overall estimation uncertainty is minimised.In the Kalman filter the radar data and the interpolated rain gauge estimates act, respectively, as "a priori estimate" and "measurement".The degree of "uncertainty" of each type of estimate constitutes a gain value (the so-called Kalman gain, K t ) at each radar grid location, and determines the proportion of each type of estimate that is used to 230 compute the merged output.The use of this gain value ensures the minimisation of the overall estimation uncertainty and is expressed as: and the optimally merged output, BAY (i.e. the a poste-

235
riori estimates y t in the Kalman filter) can be obtained from where y t , is the "unbiased" radar rainfall estimate (i.e.It can be seen that the Kalman gain is a function of the co-variances of radar and rain gauge estimation errors.When e. radar estimates have significantly higher uncertainty than the rain gauge 245 ones), the radar estimates are less trustworthy and the output estimates will be very similar to the interpolated rain gauge field.In contrast, when , the output will be closer to the radar estimates.
It is in steps b and c where the problems associated with the 250 Bayesian merging technique, and geostatistical techniques in general, arise.The (2nd-order) co-variance function that these techniques employ to characterise radar-rain gauge errors cannot well capture local singularity structures.Instead, in 2nd-order models singularities may be mistakenly re-255 garded as errors in the radar data, thus leading to higher estimated radar uncertainty, V ε RD t .As a result, the radar data will be trusted less, leading to smoother merged outputs, which are closer to the interpolated rain gauge field.

Local singularity analysis 260
Various types of hazardous geo-processes, including precipitation, often result in anomalous amounts of energy release or mass accumulation confined to narrow intervals in time and/or space.The property of anomalous amounts of energy release or mass accumulation is termed a singularity and it 265 is often associated to structures depicting fractality or multifractality (Agterberg, 2007;Cheng, 1999;Lovejoy and Mandelbrot, 1985;Schertzer and Lovejoy, 1987).Several mathematical models and methodologies have been developed to respectively characterise and treat singularities.In this work, 270 the local singularity analysis proposed by Cheng et al. (1994) has been adopted to identify and extract singularities form rainfall fields.Cheng's method, which has been widely used for estimation of geo-chemical concentrations (Agterberg, 2007;Cheng and Zhao, 2011;Cheng et al., 1994), employs the definition of coarse Hölder exponent to characterise singularities.According to this model, singularities are defined by the fact that the areal average measure (in this work, areal rainfall) centred on point x (taken as the centre of a radar pixel) varies as a power function of the area (Evertsz and Mandelbrot, 1992).This power-law relationship can be formulated as an equation (Cheng et al., 1994): where ρ(x, ) represents the density of measure (e.g.concentration of geo-data; in the context of this paper, rainfall intensity) over a square area with side-length l and associated scale ( = l/L, where L is the side-length of the largest square area under consideration) centred at a specific location x; c(x) is a constant value (in the context of this paper, a constant intensity value) at x; α(x) is the singularity index (or the coarse Hölder exponent); and E = 2 is the Euclidean dimension of the plane.
A schematic of the estimation of the constant value c(x) and singularity index α(x) from gridded data is provided in Fig. 2. For a given pixel with centre x (centre of top plot in Fig. 2), the mean rainfall intensities at different spatial scales (centred in x) can be calculated (i.e.rainfall intensities ρ 1 , ρ 2 , . . ., ρ n , respectively at scales 1 , 2 , . . ., n ).Then, the logarithms of these mean values and the associated spatial scales are compared (bottom plot in Fig. 2).The constant value c(x) and singularity index α(x) of the dataset can be derived by applying a simple linear regression analysis, where the slope and the y-intercept of the regression line correspond, respectively, to the terms (α(x) − E) and log c(x).A detailed explanation of the computation of c(x) and α(x) can be found in previous studies (Agterberg, 2012b;Chen et al., 2007;Cheng et al., 1994).It is worth mentioning that the estimation of c(x) and α(x) can be trusted only if a good linear relation is observed (i.e. if the scaling behaviour is well followed).
Going back to the definition of singularity, Eg. 3 constitutes a useful tool to decompose an areal rainfall intensity at a given location x into two components (Wang et al., 2012): (1) the background (or non-singular, NS) magnitude c(x), which is invariant as measuring scale changes and is more approximately normal than the original field; and (2) a local "scaling" multiplier, the magnitude of which changes as the measuring scale changes, according to the local singularity index α(x).When α(x) < 2, the rainfall magnitude strikingly increases as the measuring scale decreases (namely local enrichment); this corresponds to a "peak" singularity.In contrast, when α(x) > 2, the rainfall magnitude decreases as decreases (i.e.local depletion), and it is therefore a "trough" singularity.When α(x) = 2, there is no singularity: the rainfall intensity ρ(x, ) within a × area remains the same as scale changes (i.e.ρ(x, ) = c(x)).
In practice however, there is a drawback to this local singularity analysis.Because it carries out a "local" analysis, the singularity exponents are usually obtained from a small number of data samples.This increases the uncertainty of 330 the estimation of α(x).The consequence of this drawback is that the singularity is incorrectly estimated or incompletely extracted; therefore, c(x) is an unreliable or incomplete nonsingular value.To circumvent this, two numerical strategies were employed in this study.The first one involves constrain-335 ing the value of the estimated singularity exponents within a certain range.This can avoid obtaining unreasonably large or small singularity exponents.A number of ranges, symmetric to the non-singular condition (i.e.α(x) = 2), were selected for testing.They are (from the widest to narrow- 75, 2.25].These "truncated" singularity ranges were empirically chosen according to the authors' experience and the fact that the distribution of α(x) is seldom largely skewed to a specific side of α(x) = 2.It 345 can be generally expected that, the wider the range is, the more singularity information (both local enrichment and depletion) from the radar images is taken into account in the merging process.The impact of different singularity ranges in rainfall estimation and hydraulic simulation is further dis-350 cussed in Sect.3.3.
The second numerical strategy is to decompose the rainfall field using an iterative procedure (Agterberg, 2012b;Chen et al., 2007):

355
where the iterative index k = 0, 1, 2, . .., n.As k = 0, c (−1) (x) = ρ(x, ) (i.e. the original value) and c (0) (x) is the "calculated non-singular" value from the first iteration, which is equal to c(x) from the non-iterative calculation above (Eq.3).This c (0) (x) is then used as the left-hand-side 360 value of Eq. ( 4) to calculate the "non-singular" value at the next iteration, and so on.Substituting Eq. (4) into Eq.(3), one can obtain an iterative local singularity analysis equation: where The criterion to terminate the iteration procedure is when That means the singularity components have been clearly removed from the data.
such, the selected spatial scale range was deemed to represent a good balance between estimation uncertainty (which depends upon the number of samples employed in the calculations) and local feature preservation.Secondly, a scaling break at approximately 8-16 km has been reported in studies in which 1 km radar rainfall data were analysed (Gires et al., 2012;Tchiguirinskaia et al., 2012).This means that the rainfall data at spatial-scale regimes ranging from 1 to 8-16 km comply with the same or similar statistical or physical behaviour.This scaling range has also been used in other applications to represent relatively local characteristics of rainfall fields (Bowler et al., 2006).Lastly, a 10-iteration singularity analysis was applied in order to ensure that most of the singularity exponents could be extracted.The downside of conducting many iterations is the longer computational time, which may be an issue for real-time applications.Nonetheless, in practice, approximately 4-6 iterations are sufficient for effectively removing 395 most of the singularity.

Incorporation of the local singularity analysis into the Bayesian merging method
The underlying idea of the proposed method is to use the local singularity analysis to decompose each radar image into 400 a non-singular image and a singularity map before applying the Bayesian merging (step (i) in Fig. 1b).The non-singular radar image (NS-RD), which has a distribution closer to normality (thus being more suitable for Gaussian-based treatments), is merged with the point rain gauge data following 405 the Bayesian procedure.This yields a non-singular Bayesian merged field (NS-BAY).Afterwards, the singularity map is applied back and proportionally to the NS-BAY merged field (step (ii) in Fig. 1b), thus yielding a singularity-sensitive merged field (SIN).This is done by multiplying each pixel value of the NS-BAY field by the following ratio: which corresponds to the ratio difference between the original radar field (RD) and the non-singular radar field (NS-RD).In this way, a singularity-sensitive merged field (SIN), which better retains the local singularity structures embedded in the original RD field is obtained.
It is worth noting that the proposed singularity-sensitive merging method does not always increase the reliability of RD estimates.Such increase only happens when the RD estimates exhibit high singularity and thus cannot be well handled using Gaussian approximations.
A particular phenomenon which may cause problems in the application of the proposed methodology and is therefore worth highlighting is the eventual presence of singularity structures in the interpolated rain gauge field (i.e.BK field) and in the resulting supposedly non-singular Bayesian (NS-BAY) merged field.While BK fields are generally highly smooth, singularity structures may appear in the special case in which a rain gauge is located within a convective cell or 430 a local depletion.Singularity structures in the BK field may be preserved in the NS-BAY field.When this is the case, the application of the singularity map back and proportionally to the NS-BAY field may result in double-counting of singularities.This can ultimately result in a merged (SIN) product 435 with more singularities than those originally observed in the radar image.In order avoid this, a "moving window" smoothing has been applied to the BK field before it is merged with the NS-RD field.That is, each pixel value of the BK field is replaced by the mean of the original value and neighbouring 440 pixel values within a 9 km diameter (which is equal to the coarsest scale considered in the local singularity analysis).In this way singularity structures potentially present in the BK field are smoothed-off.

Case study 445
The proposed SIN merging method is tested using as case study four storm events observed in the Portobello catchment (Edinburgh, UK) during 2011 and for which radar estimates, dense rain gauge and flow records, as well as a recentlycalibrated urban drainage model were available.Portobello 450 is a coastal town located 5 km to the east of the city centre of Edinburgh, along the coast of the Firth of Forth, in Scotland (Fig. 3a).The catchment is predominantly urban, of residential character.It stretches over an area of approximately 53 km 2 , of which 27 km 2 are drained by the sewer system.Of the drained area, 46 % corresponds to impervious surfaces.This includes a small western part of Edinburgh city centre and the surrounding south-western region.The storm water drainage system is predominantly combined with some separate sewers and drains from the south-west to the north-460 east (towards the sea).

Urban storm-water drainage model
A semi-distributed model of the storm-water drainage system of the Portobello catchment, including its sewer sys-465 tem (Fig. 3b), was setup by the water utility of the area in the commercial modelling package InfoWorks CS v13.0.In this model the whole catchment surface is split into subcatchment units through which rainfall is applied (within each sub-catchment rainfall is assumed to be uniform).Each 470 sub-catchment comprises a mix of pervious and impervious surfaces whose runoff drains to a common outlet point, which corresponds to an inlet node of the sewer system (i.e. a gully or a manhole).Each sub-catchment is characterised by a number of parameters, including total area, 475 length, slope, proportion of each land use, amongst others.Based upon these parameters, the runoff volume at each sub-catchment is estimated using the NewUK rainfall-runoff model (Osborne, 2001).The estimated runoff is then routed to the sub-catchment outlet using the Wallingford (double linear reservoir) model (HR Wallingford, 1983).Sewers are modelled as one-dimensional conduits and flows within them are simulated based on the full de Saint-Venant equations (i.e.fully-hydrodynamic model).
The Portobello model contains a total of 1116 subcatchments, with areas ranging between 0.02 and 24.42 ha and a mean area of 2.3 ha.Sub-catchment slopes range from 0.0 to 0.63 m m −1 , with a mean slope of 0.031 m m −1 .The model of the sewer system comprises 2917 nodes and 2907 conduits, in addition to 14 pumps.The total length of modelled sewers is 250 km.The sewer system ranges in height from 186.6 mAOD at Comiston to 3.8 mAOD along the Firth of Forth. 2 % of the modelled pipes have a gradient between 0.1 and 0.25 m m −1 , 55 % have a gradient between 0.01 and 0.1 m m −1 , and 43 % have a gradient < 0.01 m m −1 .
Following UK standards (WaPUG, 2002) and using solely rain gauge data as input, the model of the Portobello catchment was manually calibrated in 2011 based upon three storm events recorded during the medium term flow survey described below.For the three storm events used for model calibration, the following mean performance statistics were obtained: Nash-Sutcliffe Efficiency of 0.5782; Root Mean Square Error of 0.0373 m 3 /s; coefficient of determination (R 2 ) of 0.7756; and regression coefficient (β), which provides a measurement of conditional bias, of 0.8826.

Local monitoring data (medium term survey)
Local rainfall and flow data were collected in the Portobello catchment through a medium term flow survey carried out between April and June 2011.The survey comprised 12 tipping bucket rain gauges and 28 flow monitoring stations (each comprising a depth and a velocity sensor, based upon which flow rates were estimated).Both rain gauge and flow records were available at a temporal resolution of 2 min.However, rain gauge records were linearly interpolated to 5 min, in order to ensure agreement with the temporal resolution at which radar estimates were available (see Sect. 3.1.3).The location of the local flow and rain gauges is shown in Fig. 3b.

Radar rainfall data
The Portobello catchment is within the coverage of C-band radars operated by the UK Met Office (Fig. 3b).Radar rainfall estimates for the same period as the local flow survey (i.e.April-June 2011) were obtained through the British Atmospheric Data Centre (BADC).These estimates were available at spatial and temporal resolutions of 1 km and 5 min, respectively, and correspond to a quality controlled multiradar composite product generated with the UK Met Office Nimrod system (Golding, 1998), which incorporates corrections for the different errors inherent to radar rainfall measurements, including identification and removal of anoma-lous propagation (e.g.beam blockage and clutter interference), attenuation correction and vertical profile correction (for a full description of the Nimrod system, the reader may refer to Harrison et al. (2000Harrison et al. ( , 2009))).

Storm events selected for analysis 535
During the monitoring period (April-June 2011), four relevant storm events were captured which comply with UK standards for calibration and verification of urban drainage models (i.e.these events have instantaneous rainfall rates > 5 mm h −1 and accumulation > 5 mm) (Gooch, 2009).Three 540 of these events (referred to as Storms 1, 2 and 3) were used for calibration of the urban drainage model of the Portobello catchment (following UK standards, as mentioned above).In this study, all four storm events, including one not used in the calibration of the model (Storm 4), were used to test the pro-545 posed singularity-sensitive merging method.The dates and main statistics of the four selected events are summarised in Table 1.It is worth mentioning that, given the response time of the catchment, as well as inter-event time definition thresholds (IETD) recommended in the literature (Guo and 550 Adams, 1998), a minimum IETD of 6 h was used as criteria to differentiate rainfall events.
As can be seen in Table 1, Storms 1, 2 and 4 are of relatively short duration, whilst Storm 3 is of much longer duration.Moreover, well-structured storm cell clusters crossing 555 the catchment area can be found in Storms 1, 3 and 4, but not in Storm 2. This is reflected in the lower RG peak intensity observed in Storm 2.

Evaluation methodology
The performance of the proposed singularity-sensitive 560 Bayesian method (SIN hereafter) is assessed by intercomparison against radar (RD), rain gauge (RG) and blockkriged (BK) interpolated RG estimates, as well against adjusted estimates resulting from the original Bayesian (non singularity-sensitive) technique (BAY) and the commonly-565 used mean field bias (MFB) adjustment method.It is important to note that, in this work, the MFB was implemented in a relatively dynamic way by computing a sample cumulative bias (B t ) at each time step as B t = RG catchment / RD co-located , where RG catchment and 570 RD co-located represent the rain gauge and the co-located radar grid rainfall estimates over the experimental catchment during the last hour.The MFB adjusted estimates are obtained by multiplying the bias (B t ) at each particular time step by the original rainfall field; that is, MFB = B t • RD.

575
Two evaluation strategies were applied: 1. Through analysis of the different rainfall estimates, using as main reference local rain gauge records, while also inter-comparing the behaviour of other estimates.
feeding the different rainfall estimates as input to the hydraulic model of the Portobello catchment and comparison of these with available flow records.Note that the RG estimates were applied to the model using Thiessen polygons.
Both evaluation strategies have inherent limitations which are next described.However, they provide useful and complementary insights into the performance of the proposed merging method.
The first strategy is a natural and widespread way of as-590 sessing the performance of rainfall products.However, the fact that all precipitation estimates entail errors and that the true rainfall field is unknown, in addition to the differences in the spatial and temporal resolutions of RG and RD estimates (and the resulting merged rainfall products) renders 595 any direct comparison of rainfall estimates imperfect (Brandes et al., 2001).In the particular case of rain gauge records, which are used as main reference in the present investigation, errors can arise from a variety of sources.In order of general importance, systematic errors common to all rain gauges in-600 clude errors due to wind field deformation above the gauge orifice, errors due to wetting loss in the internal walls of the collector, errors due to evaporation from the container and errors due to in-and out-splashing of water.In addition, tipping bucket rain gauges, such as the ones used in this investiga-605 tion, are known to underestimate rainfall at higher intensities because of the rainwater amount that is lost during the tipping movement of the bucket (La Barbera et al., 2002).This is a systematic error unique to the tipping bucket rain gauge and is of the same order of importance as wind induced losses.
Besides these systematic errors, tipping bucket rain gauge records are also subject to local random errors, mainly related to their discrete sampling mechanism (Habib et al., 2008).The order of magnitude of the two main error sources associated to tipping bucket rain gauges is 2-10 % for wind induced losses (depending on wind speed, weight of precipitation and gauge construction parameters) (Sevruk and Nešpor , 1998) and up to 10 % for rainfall intensities of 100 mm/h and 20 % for intensities of 200 mm/h, for water losses during the tipping action (Luyckx and Berlamont, 2001;Molini et al., 620 2005).Given that the storm events under investigation were not extremely heavy and that a basic quality-control was applied to ensure the quality of rain gauge measurements (e.g.manual comparison between neighbouring rain gauge data), the rain gauge records used in this study can be deemed acceptable as reference for the evaluation of other rainfall products.However, they are not perfect and the reader must keep in mind that the true rainfall field remains unknown.
The second strategy (i.e.hydraulic evaluation) allows some of the limitations of the rainfall evaluation strategy to be overcome, and is particularly useful when dense flow records are available, as is the case in the Portobello catchment.However, it has two main deficiencies: the fact that flow records (obtained based upon depth and velocity measurements) used in the evaluation contain errors, and the fact 635 that the hydraulic modelling results encompass uncertainties from different sources in addition to rainfall input uncertainty (Deletic, 2012;Kavetski et al., 2006).In this regard, it is worth reminding the reader that the hydraulic model used in this study was calibrated using as input RG data.In fact, the 640 model was calibrated based upon the data of Storms 1 -3, which are used for testing in the present paper (Storm 4 is the only event not used in the calibration of the model).Since the model was "attuned" for RG inputs, it favours all RG-derived (and RG-emulating) rainfall estimates.Moreover, the rela-645 tively coarse spatial resolution of rainfall inputs (in this case ∼ 1 RG/4.4 km2 ) may have led to further biases in the model (Kavetski et al., 2006).It would be desirable to re-calibrate the hydraulic model using as input the different rainfall products analysed in this study.However, this would entail a sig-650 nificant amount of work which falls outside of the scope of the present study.In spite of these limitations, we believe that the hydraulic evaluation strategy using the currently available hydraulic model still provides useful insights into the performance of the proposed SIN merging method in relation to 655 other rainfall estimates, and complements the findings of the rainfall analysis.

Methodology for analysis of rainfall estimates
The performance of the SIN rainfall products in relation to other rainfall estimates (including RD, RG, BAY and MFB) 660 is evaluated in terms of accumulations and rainfall rates at the areal level (i.e. at a scale corresponding to the area over which the Portobello catchment stretches) and at individual point gauge locations.In addition, a qualitative assessment of the spatial structure of the different (gridded) rainfall prod-665 ucts is carried out based upon visual inspection of images of the rainfall fields at the time of areal average peak intensity.
In view of the high density and coverage of the RG network over the Portobello catchment, the areal average RG estimates in the areal level analysis, are assumed to be a good 670 approximation of the "true" areal (average) rainfall over the experimental catchment (i.e. the areal reduction effect is expected to be minor (Bell, 1976)) and are therefore used as reference to evaluate the performance of the different areal gridded rainfall estimates.The areal analysis includes com-675 parison of event areal average accumulations and peak intensities, as well as comparison of intensities throughout the duration of each event (through scatterplots using RG areal average intensities as reference).
In the analysis of rainfall estimates at rain gauge point lo-680 cations a cross-validation strategy was adopted and three performance statistics are used.The cross-validation strategy, also referred to as "leave-one-out", is an iterative method in which, at each iteration, data from one RG site is omitted from the calculations and the value at the "hidden" (i.e. formance statistics are then computed from the comparison between the estimated and the known (but not used) values (Velasco-Forero et al., 2009).The following three performance statistics are used in the present study.Firstly, a sample bias ratio (B) is used to quantify the cumulative bias between gridded rainfall estimates (i.e.RD, BK, MFB, BAY and SIN) and RG estimates at each RG location over the event period under consideration.B = 1 means no cumulative bias between the RG and the given gridded rainfall 695 estimates (i.e.equal rainfall accumulation recorded by RG and the gridded product at the given gauge location); B > 1 means that the accumulations of the gridded estimates at the point locations are greater than those recorded by RG, and B < 1 means the opposite.In addition to the comparison of rainfall accumulations, a simple linear regression analysis is applied to each pair of "instantaneous" (rain rate) point RG records and the co-located gridded estimates.The results of the regression analysis are presented in terms of R 2 (coefficient of determination) and β (regression coefficient, i.e., the slope or gradient of the linear regression).These two measures provide an indication of how well RG rates are replicated by the different rainfall estimates at each gauging location, both in terms of pattern and accuracy.The R 2 measure ranges from 0 to 1 and describes how much of the "observed" (RG) variability is explained by the "modelled" (RD/BK/adjusted) one.In practical terms, R 2 provides a measurement of the similarity between the patterns of the observed (i.e.RG) and "modelled" (i.e.gridded estimates) rainfall time series at a given gauging location.However, systematic bias (under-or over-estimation) of the modelled estimates cannot be detected from this measure (Krause et al., 2005).The regression coefficient, β, is therefore employed to provide this supplementary information to the R 2 .β ≈ 1 represents good agreement in the magnitude of the rainfall rates recorded by RG and those of the gridded estimates; β > 1 means that the rain rates associated to the gridded estimate are higher in the mean (by a factor of β) than those recorded by RG; and β < 1 means the opposite (i.e.rain rates of gridded estimates are lower in the mean than RG ones).

Methodology for analysis of hydraulic outputs
A qualitative analysis of the hydraulic outputs is carried out based upon visual inspection of recorded vs. simulated flow hydrographs (for the different rainfall inputs) at different points of the catchment.Furthermore, similar to the rainfall analysis, a simple linear regression analysis is applied to each pair of recorded and simulated flow time series (at each flow gauging location).The performance of the associated hydraulic simulations is evaluated using the R 2 and β statistics obtained from the linear regression analysis.In addition, the "weighted" coefficient of determination (R 2 w ) is employed to quantify the joint performance of hydrological efficiency.This measure is defined as (Krause et al., 2005): where higher R 2 w values correspond to better hydraulic per-740 formance.
In order to minimise the influence of the errors in the flow measurements, the available flow records were qualitycontrolled (QC) before carrying out the statistical analysis of hydraulic outputs.The QC was carried out following UK 745 guidelines (WaPUG, 2002).It included analysis of depth vs. flow scatterplots at each monitoring location (the shape and spread of the resulting scatterplots provides insights into the quality and consistency of depth and velocity records at each site), as well as visual inspection of the observed hy-750 drographs at each location.Whenever a flow monitor was deemed unreliable, it was manually removed from the analysis.Likewise, with the purpose of preventing systematic hydraulic modelling errors from affecting results, whenever the model was found to be unable to replicate the recorded flows 755 at a given location, the given flow monitor was also manually removed from the analysis.This left us with a total of 16 flow monitors for analysis.Moreover, all records associated with depth measurements below 0.1 m were left out when estimating performance statistics; this is due to the fact that at low 760 depths, both velocity and depth records become unreliable (the 0.1 m threshold was adopted based upon UK guidelines and recommendations from studies focusing on the performance of flow gauges (Marshall and McIntyre, 2008)).

Areal rainfall estimates
Table 2 shows the areal average (AVG) accumulations and peak intensities for the different rainfall products for the four storm events under consideration.As can be seen, the differ-770 ence between RD and RG event areal average accumulations is generally small, yet it is event-varying.For Storms 1 and 2, the areal average RG totals are slightly overestimated by the RD estimates, while for Storms 3 and 4 they are underestimated, with the underestimation being largest in Storm 3. In 775 terms of areal average peak intensities, RD estimates appear to consistently underestimate the areal average peak intensities recorded by RG.The relative difference between RG and RD areal average peak intensities is approximately 20-30 % for Storms 1 and 2, and it is as high as 60 % for Storms 3 and 780 4.This indicates that the RD estimates could not satisfactorily capture instantaneous rainfall rates, particularly high rainfall rate values, and corroborates the need for dynamic adjustment of RD estimates using local RG measurements.
As would be expected, the BK estimates exhibit areal av-785 erage accumulations and peak intensities similar to those of the RG.Small differences are observed (in general BK values are slightly lower than RG ones) which can be generally attributed to the area-point rainfall differences (Anagnostou and Krajewski, 1999)).These differences become more evident when analysing results at individual gauge locations (in next section).When looking at the adjusted rainfall products (i.e.MFB, BAY and SIN), it can be seen that all of them can improve the original RD estimates, but the degree of improvement is different for each method.As expected, the MFB successfully reduces the difference in event areal average accumulations (i.e.bias), leading to areal average accumulations close to those recorded by RG.In terms of peak intensities, the MFB method leads to some improvement, but the resulting rainfall estimates for urban hydrological applications.
The BAY estimates show the least improvement in terms of event bias, with a general tendency to underestimate RG areal accumulations, which is even more marked than for BK estimates.This is particularly the case in Storms 3 and 4, 810 in which strong singularity structures, as represented by the high frequency of α values different from 2 (Fig. 4), were observed.This tendency to underestimate can be attributed to a combined effect of the BAY method "over-trusting" the BK estimates (which show a slight underestimation tendency at the areal level), in addition to smoothing off the singularity structures (often associated to strong intensities) originally present in the RD image.With regard to peak intensities, the BAY estimates display a larger improvement than the MFB ones, which demonstrates the benefits of more dynamic and spatially-varying adjustment methods.However, the areal average peak intensities of the BAY estimates still underestimate the "true" (RG) areal peak intensities.Lastly, the SIN estimates, particularly SIN ranges 1-4, exhibit very good performance: both areal average accumulations and peak intensities of SIN estimates are close to those of RG, and no systematic over-or underestimation is observed in SIN estimates.The better performance of SIN estimates suggests that the singularity analysis can in fact improve the original BAY merging method.With regard to the impact of the singularity range, it can be seen that, as the range becomes narrower (from SIN1 to SIN5), the areal accumulations and peak intensities tend to be higher.A particular steep increment is observed in the areal accumulations and peak intensities of SIN5, in relation to the other singularity ranges, which results in a slight overestimation of accumulations and peak intensities, as compared to the RG estimates.This is especially apparent in Storms 3 and 4. The increment in rainfall accumulation and peak intensities as the singularity range becomes narrower can be explained by the fact that at narrower ranges, only part of the singularity structures are re-moved before merging, and a big proportion of them remains in the radar image.This is evident from Fig. 4, where it can be seen that for Storms 3 and 4 a significant number of singularity exponents spread beyond the narrowest singularity 845 range (i.e.SIN5: α ∈ [1.75, 2.25]).The problem arises because some singularity features remain in the radar image before the BAY merging is applied, and these may be partially preserved throughout the merging process.Afterwards, when the extracted singularity component is applied back 850 and proportionally to the merged rainfall field, it may interact (in a non-linear fashion) with the singularity structures preserved throughout the BAY merging, thus leading to an overestimation of extremes (whether these are enrichments or depletions).These results suggest that a better approach 855 is to be sure to remove most of the singularity structures before carrying out the merging.Therefore, very narrow ranges such as SIN5 should be avoided.In fact, it can be seen that the intermediate range of SIN3 (i.e.α ∈ [1, 3]) covers most singularity exponents (see Fig. 4).Indeed, using wider sin-860 gularity ranges leads to very similar results as those obtained when using the SIN3 range (notice the similarity between SIN1, SIN2 and SIN3 estimates in Table 2).This suggests that the singularity range of [1,3] and 4).Unlike BK estimates, the difference between areal RD and RG peak intensities is too large to be entirely explained by areal-point differences.The fact that the RD esti-880 mates display relatively good performance in (long-duration) accumulations but not in (short-duration) instantaneous intensities corroborates the claim that RD estimates fail to capture the short-term dynamics of small-scale rainfall.As mentioned in Sect. 1, the main reason for this lies in the use of 885 a static Z−R conversion function, which represents a statistical compromise for the range of rainfall rates that frequently occur (whereas the occurrence of very small and large intensities is relatively rare).Concerning the MFB method, from Fig. 5 it can be seen that it fails to satisfactorily improve RD 890 instantaneous rainfall rates; this is particularly evident at high intensities, at which, similarly to RD estimates, MFB estimates perform poorly.An accurate representation of peak intensities is of outmost importance in the modelling and forecasting of urban pluvial flooding.This confirms that, being 895 a 1st-order technique, the MFB adjustment method may be insufficient for urban-scale applications.In contrast, it can be seen that over-and underestimation errors in RD estimates can be improved by 2nd-or higher-order adjustment techniques, such as BAY and SIN.In fact, in terms of instantaneous rainfall rates the BAY estimates display a significantly better performance than the MFB ones.However, the BAY estimates still fail to properly reproduce the highest intensities.These shortcomings of the BAY method seem to be overcome by the incorporation of the singularity analy-905 sis: indeed, the SIN estimates exhibit the best overall performance, particularly at peak intensities.In agreement with the results displayed in Table 2, in Fig. 5 it can be seen that as the singularity range becomes narrower, rainfall estimates with slightly higher intensities are generated.Moreover, it can be noticed that wider singularity ranges lead to more conservative results and appear to be a good choice.

Rainfall estimates at gauging locations
The aforementioned features of the different rainfall estimates are further highlighted through analysis at each rain 915 gauge location; the associated statistics, including sample bias (B), regression coefficient (β) and coefficient of determination (R 2 ), are summarised in Fig. 6.As expected, the RD estimates (before adjustment) display the largest differences from point RG estimates: in general, 920 they possess the largest cumulative bias (B) and the lowest R 2 , and their statistics show great variability.Moreover, the distribution of the β values indicates that RD estimates tend to largely underestimate RG instantaneous rainfall rates.This is the case for all storm events, except for Storm 2, for which 925 rainfall intensities were low on average.
Similarly to the results of the areal (average) analysis, the individual-site BK estimates display the closest behaviour to the RG ones.This is of course expected given that the BK estimates are obtained by simple interpolation of point RG 930 data.It can be seen that the BK estimates are nearly unbiased (B is very close to 1) and possess the highest R 2 medians (i.e.closest to 1), as well as the narrowest R 2 boxes and whiskers.However, when looking at the distribution of β values of the BK estimates, it can be seen that most of the time the whole 935 boxes and whiskers are below 1.This reflects a systematic underestimation of rainfall rates at point gauging locations, which is discussed below.
With regards to the adjusted rainfall estimates, the MFB method is found to bring original radar estimates slightly 940 closer to RG ones, but the improvement seems insufficient.As expected, the main improvement of MFB estimates is found in the bias (B), which is significantly reduced (thus becoming closer to 1).In terms of instantaneous rainfall rates, the improvement provided by MFB is very limited.This is reflected in the low R 2 values and in the poor β scores, which remain remarkably close to those of the original RD estimates.Similar to the results of the areal analysis, these results suggest that the MFB adjustment method is insufficient for urban-scale applications, in which small scale rainfall dy-950 namics are critical.
When looking at the statistics of the BAY estimates, it can be noticed that these behave similarly to the BK ones: their bias is also small, the R 2 is generally high and the β values are systematically below 1.The similarity in the behaviour of 955 BK and BAY estimates suggests that for the selected events, the BAY method tends to trust the (smooth interpolated) BK estimates more than the RD estimates.As explained in Sects. 1 and 2, this is the main shortcoming of the BAY method.The systematic underestimation of rainfall rates ob-960 served in BK and BAY estimates (reflected by β values systematically below 1) can be partially explained by the areal reduction effect.However, considering that the individualsite comparison was conducted using instantaneous rainfall estimates for very short time intervals (i.e. 5 min), one would 965 expect the tendency of the areal-point differences to be of higher randomness, rather than of a systematic nature.This suggests that the systematic underestimation may be a joint consequence of the areal reduction effect and of the underlying 2nd-order approximation (which smooths off some local 970 extreme magnitudes).
With regards to the SIN estimates, it can be seen that their bias is small (close to 1) and that the distribution of their R 2 values is somewhere between that of the BAY and RD estimates.This indicates that, as compared to the original 975 BAY estimates, the SIN method incorporates more spatial features from the RD estimates throughout the merging process, while retaining the accuracy of the RG estimates.In terms of β, it can be seen that although the median values are usually below one, they are generally much closer to one 980 than other rainfall estimates.Moreover, their distribution is more variable than that associated to BK and BAY estimates (which display a systematic underestimation).In line with the results of the areal analysis, this serves to further highlight two important features of the proposed SIN method: 985 it generally respects the commonly-observed areal reduction effect, and it integrates more small scale randomness from RD data in the data merging process.Regarding the singularity ranges, a similar behaviour can be observed for all of them, although a slight tendency can be observed for the bias 990 (B) and β values to increase as the singularity range becomes narrower.This is in agreement with the results of the areal analysis and suggests that working with intermediate ranges, as opposed to very narrow ones which can truncate the singular structures is advisable.It can be seen that the spatial structure of the BK rainfall 1005 field (fully based upon rain gauge data) is highly symmetric and smooth, and is rather unrealistic.With regards to the adjusted rainfall products (MFB, BAY and SIN), it can be noticed that the proportion of radar (RD) and BK interpolated rain gauge features that are preserved varies according to the method.The MFB fields fully inherit the spatial structure of the RD fields; the only change is that the actual intensity values are scaled up or down by an areal ratio derived from the sample bias between mean rain gauge and radar rainfall estimates.In agreement with the quantitative 1015 results presented above, it can be seen that the structure of the BAY peak rainfall fields is often similar to that of the BK ones and is smoother than the original RD image.Singularity structures are often present in rainfall fields during peak intensity periods (such as the ones shown in Fig. 7).As 1020 explained in Sect.2, the presence of these structures causes the RD fields to be considered highly uncertain and therefore these are less trusted in the BAY merging process.This results in BAY peak intensity merged fields closer to the BK ones, instead of to RD ones.Some spatial features from RD 1025 can be still observed in the BAY fields, for example, in the lower-right area of the BAY image of Storm 1 and in the middle-left area of the BAY image of Storm 3 (see Fig. 7, top).However, these features appear to be much smoother and spreading over a larger area, as compared to their struc-1030 ture in the original RD image.As compared to the BAY peak rainfall fields, the SIN fields display less smooth and more realistic structures which preserve more features of the original RD fields.In general, the inspection of the snapshot images of the different rainfall products confirms the findings 1035 of the areal and point gauge analyses regarding the ability of the SIN method to better preserve the singularity structures present in rainfall fields (and captured by RD) throughout the merging process, as compared to the original BAY merging method. 1040

Hydraulic modelling results
Fig. 8 and Fig. 9 shows example observed vs. simulated flow hydrographs for the different rainfall inputs at four gauging locations (see locations in Fig. 3) during Storm 3. Note that Storm 3 is an event with high rainfall accumulations, rainfall 1045 rates and strong singularity structures.Figure 10 summarises the performance statistics resulting from the simple linear regression analysis (i.e.β, R 2 and R 2 w ) conducted at each flow gauging station for each storm event.
From the hydrographs in Figures 8 and 9 it can be seen that, in terms of pattern and timing, all simulated flows (resulting from the different rainfall inputs) are generally in good agreement with observations.This indicates that all rainfall products, including the RD before adjustment, can well capture the general dynamics of rainfall fields.The main difference between the simulated flows lies in their ability to reproduce flow peaks, which in turn is a function of the ability of the different rainfall estimates to reproduce peak rainfall rates in terms of magnitude, timing and spatial distribution.In line with the results of the rainfall analysis, the flow 1060 hydrographs associated with BK estimates are close to the ones associated with RG records; however, the former display a smoother behaviour and generally lead to flow peaks that are lower than the recorded ones and the ones resulting from RG inputs.The RD outputs can well match some of the 1065 observed flow peaks, but significantly underestimate others (e.g.flow peak at around 23 h in FM10); this can be attributed to the underestimation of peak intensities observed in RD estimates.The MFB associated flows show little improvement over the RD ones and in many cases they even lead to a worse 1070 performance (e.g.see overestimation of peak flows at around 02:00 UTC in FM14 and FM19).This is further confirmed by the statistics in Fig. 10 and corroborates the claim that the application of a "blanket" MFB correction over the area of interest is insufficient for urban applications.The BAY outputs, 1075 on the other hand, show a consistent improvement over the RD outputs.Nonetheless, in agreement with the results of the rainfall analysis, the BAY estimates behave similarly to the BK ones and lead to smooth flow peaks which often underestimate observations.The SIN outputs also show a consistent 1080 improvement over RD estimates, but, unlike BAY outputs, the SIN outputs do not smooth off flow peaks and instead show a better ability to reproduce these, sometimes leading to a better match of observed flow peaks than RG associated outputs (which were used as input for the calibration of 1085 the model).With regards to the singularity ranges, it can be seen that the differences observed in the SIN1-SIN5 rainfall estimates are mostly filtered out when converting rainfall to runoff.As a result, the flow outputs of all SIN estimates are very similar and their hydrographs can barely be differenti-1090 ated.
The preliminary conclusions drawn from the visual inspection of the selected hydrographs are corroborated by the statistics in Fig. 10.As would be expected, the RG outputs generally show the best performance in all statistics (except 1095 for Storm 2).It is nonetheless noteworthy that the β values for RG estimates (for which the model was calibrated) are generally below 1.This reveals a slight bias of the model to underestimate flows and partially explains the fact that β values associated with all rainfall inputs are mostly below 1100 1.With regards to the BK (i.e.interpolated RG) associated outputs, the tendency to underestimate, which is observed in the rainfall analysis, becomes even more evident in the hydraulic outputs: BK's β values are significantly lower than 1 and lower than the β values associated with RG estimates.

1105
Different from the results of the rainfall analysis, in which those products closest to RG estimates (including BK) displayed the best performance in terms of R 2 , in the hydraulic analysis BK outputs generally lead to a deterioration in R 2 values (see statistics of Storms 2, 3 and 4).This suggests 12 L.-P.Wang et al.: Singularity-sensitive data merging that the smoothing caused in the BK interpolation affects the small scale dynamics of rainfall fields, leading to poor representation of associated flow dynamics.Contrary to RG associated outputs, RD flow estimates generally display the worst performance.In agreement with the results of the rainfall analysis, RD outputs show a tendency to largely underestimate flows (as indicated by β values well below one and much lower than those obtained for RG outputs).Moreover, they display relatively low R 2 and associated R 2 w values.Nonetheless, a special case is observed in Storm 2, when RD estimates, which in the rainfall analysis displayed the poorest performance, yielded the best flow simulations, thus emphasising the added value that RD estimates can provide, as well as the complementary information provided by the hydraulic evaluation strategy.In the cases in which RD outputs perform poorly (i.e. in Storms 1, 3 and 4), all adjusted rainfall estimates lead to improvements over RD hydraulic results, with the degree of improvement varying according to the adjustment technique.In Storm 2, when RD outputs displayed the best performance, the different merging methods showed to retain different degrees of RD features.Overall, it can be seen that the MFB estimates provide little improvement over the original RD estimates.In the cases in which RD led to systematic underestimation of flows, the MFB estimates managed to slightly reduce this underestimation by bringing β values closer to 1, as compared to those of the original RD estimates.However, in Storm 2, in which RD outputs performed best, MFB caused a large deterioration of β values, whereas the SIN estimates managed to keep β scores closer to 1.In terms of R 2 , the MFB estimates do not provide much improvement and can actually lead to a deterioration of this statistic (e.g.Storms 2 and 4), suggesting that the application of the MFB adjustment can alter the spatial-temporal structure of the original RD fields.The BAY outputs show a greater improvement than MFB, particularly in terms of R 2 .Nonetheless, in agreement with the rainfall analysis and with the visual inspection of hydrographs, the BAY outputs behave remarkably similarly to BK ones.One of the main features of the BAY outputs is that they lead to systematically lower flows than RG estimates (note β values consistently lower than those of RG outputs).This confirms the smoothing of rainfall peaks that occurs when 2nd-order approximations are applied.Lastly, it can be seen that the SIN outputs display the greatest improvement over original RD outputs, both in terms of correcting systematic bias, as well as in terms of well reproducing rainfall and associated flow patterns.As compared to MFB and BAY outputs, the SIN outputs generally display β values closer to 1, higher R 2 values (sometimes even higher than those of RG outputs) and consequently higher R 2 w values.The better performance of SIN hydraulic outputs over BAY ones, particularly in terms of β, provides an a posteriori confirmation that it was right, in the SIN method, to view the singularities in the radar field as actual features of the real rainfall.Same as was observed in the hydrographs in Fig. 8 and 9, the performance of the dif-ferent SIN ranges is very similar.In line with the results of the rainfall analysis, SIN5 shows a slight tendency towards higher flows, which sometimes resulted in better hydraulic statistics.Nonetheless, the differences are small and based upon the findings of the rainfall analysis, the adoption of 1170 a wider SIN range and removal of most singularities appears to be a more conservative option.However, this aspect must be further investigated using a wider range of storm events and pilot catchments.
Regarding the difference in hydraulic performance be-1175 tween the events used for model calibration (i.e.Storms 1-3) and the independent event (i.e.Storm 4), it can be seen that the statistics of the hydraulic outputs during Storm 4 are generally lower than for the other three storm events (Fig. 10).This can be expected, given that the model was "tuned" to 1180 give a good fit for the calibration events.However, as discussed above, the general features and relative performance of the hydraulic outputs associated to the different rainfall inputs was generally consistent in all storm events under consideration.The particular differences that were observed 1185 were due to the nature of a given storm and not to the fact that a given event was used for model calibration or not.

Conclusions and future work
In this paper, a new gauge-based radar rainfall adjustment method was proposed, which aims at better merging rainfall 1190 estimates obtained from rain gauges and radars, at the small spatial and temporal scales characteristic of urban catchments.The proposed method incorporates a local singularity analysis into the Bayesian merging technique (Cheng et al., 1994;Todini, 2001).Through this incorporation, the merging 1195 process preserves the fine-scale singularity (non-Gaussian) structures present in rainfall fields and captured by radar, which are often associated with local extremes and are generally smoothed off by currently available radar-rain gauge merging techniques, mainly based upon Gaussian approxi-1200 mations.
Using as case study four storm events observed in the Portobello catchment (53 km 2 ) (Edinburgh, UK) in 2011, the performance of the proposed singularity sensitive Bayesian data merging (SIN) method, in terms of adjusted rainfall esti-1205 mates and the subsequent runoff estimates, was evaluated and compared against that of the original Bayesian data merging (BAY) technique and the widely-used mean field bias correction (MFB) method.This analysis clearly brought out the benefits of introducing the singularity-sensitive method.The 1210 results suggest that the proposed SIN method can effectively identify, extract and preserve the singular structures present in radar images while retaining the accuracy of rain gauge (RG) estimates.This is reflected in the better ability of the SIN method to reproduce instantaneous rainfall rates, rain-1215 fall accumulations and associated runoff flows.This method clearly outperforms the commonly used MFB adjustment, which simply fails to reproduce the dynamics of rainfall in urban areas, and the original BAY method, which shows an overall good performance but smooths off peak rainfall magnitudes, thus leading to underestimation of runoff extremes.
In this study the sensitivity of the SIN results to the "degree" of singularity that is removed from the radar image and preserved throughout the merging process was also tested.While the impact of it was found to be generally small, the results suggest that partially removing singularities could have a negative impact on the results.Therefore, removing most singularity exponents from the original radar image is advisable.
While the proposed singularity method has shown great potential to improve the merging of radar and rain gauge data for urban hydrological applications, further testing including more storm events and pilot catchments is still required in order to ensure that the results are not case specific and to draw more robust conclusions about the applicability of the proposed method.Other aspect on which further work is recommended are the following: • The current version of the singularity-sensitive method shows a slight tendency to overestimate rainfall rates and accumulations.This is likely to be due to one of two aspects, or a combination of them: -In the eventual case in which a rain gauge is located within the core of a convective cell, the resulting interpolated (block-kriged; BK) field may end up having singularity structures and, as explained in Sect.2.3, this may ultimately lead to "double-counting" of singularities.For the time being, a moving-window smoothing has been applied on the BK field before it is merged with the nonsingular radar field, so as to remove singularity structures potentially present in the BK field.While this has proven to be an acceptable solution, we believe it can be further refined.Other methods for dealing with this particular problem have been tested, including extraction of singularities from the point rain gauge records using the singularity exponents derived from the co-located radar pixels, and extraction of singularities from the BK field using local singularity analysis.However, these have proven unstable and highly uncertain.Further work to better deal with this issue is required.
-The asymmetric distribution of singularity exponents and the numerical stability of singularity extraction from a small set of data samples.This drawback could be improved by forcing the mean of non-singular components to remain equal to the original radar estimates (Agterberg, 2012a).Alternatively, other techniques for singularity identification and extraction could be used.For example, the wavelet transformation (Kumar and Foufoula-Georgiou, 1993;Mallat and Hwang, 1992;Robertson et al., 2003;Struzik, 1999), Principal Component Analysis (Gonzalez-Audicana et al., 2004;Zheng et al., 2007) and Empirical Mode Decomposition (Nunes et al., 2003(Nunes et al., , 2005) ) techniques are 1275 widely recommended in the literature.
• Given that the proposed singularity-sensitive merging method is particularly intended to improve rainfall estimates for (small-scale) urban areas, it would be interesting to test it using higher spatial-temporal resolution 1280 data (e.g. from X-band radars).
Lastly, a suggestion often made to us and therefore worth briefly discussing is to use a transformation in order to bring the distribution of the radar field closer to normality before the merging (be it with the Bayesian or other geo-statistical 1285 method) is conducted.However, doing this would somehow miss the point of the proposed method.The key point here is that for a non-Gaussian structure, moments beyond the second order are important, as each brings new information worth preserving.To create a more 'normal' field is not the 1290 purpose of the singularity extraction; instead, it is the consequence after removing singularities from the rainfall fields, which can be physically associated with abnormal energy concentration, such as 'convective' cells, and which in the proposed method are set aside to ensure their preservation 1295 throughout the merging process.
model, and for their constant support with the hydraulic simulations.Thanks are also due to the UK Met Office and the BADC (British Atmospheric Data Centre) for providing Nimrod (radar) data, to Innovyze for providing the InfoWorks CS software, and to Dr. Cinzia Mazzetti and Prof. Ezio Todini for making freely available to us the 1310 RAINMUSIC software package for meteorological data processing.Lastly, the authors would like to thank the reviewers, Dr. Scott Sinclair and Dr. Cinzia Mazzetti, and the Editor, Dr. Uwe Ehret, for their insightful and constructive comments which helped improved the manuscript significantly.
used as the a priori estimate in the 240 Kalman filter. 765

995
Spatial structure of rainfall fieldsSnapshot images of the different gridded rainfall products at the time of peak areal intensity for the 4 storm events under consideration are shown in Fig.7.Due to space constraints, images of only one of the SIN ranges (the intermediate one, 1000 SIN3: α ∈ [1, 3]) are shown in this figure.As mentioned in the previous sections, the SIN3 range covers most of the sin-gularity indices and consistently led to good results for the storms under consideration.

Fig. 3 .
Fig. 3. Portobello catchment (a) general location; (b) sensor location, sewer network and radar grid over the catchment.Flow gauges FM 3, 10, 14 and 19 (the circled round markers in (b), respectively located at the up-, mid-and downstream parts of the catchment) are particularly selected for visual inspection in Sect.3.3.2.

Fig. 4 .Fig. 5 .
Fig. 4. Histogram of singularity exponents (α) at the time of areal peak intensity for the four selected storm events.

Fig. 6 .
Fig. 6.Boxplots displaying the distribution of sample bias ratio (B) (left column), regression coefficient (β) (middle column) and coefficient of determination (R 2 ) (right column) estimated between the different gridded rainfall estimates and RG records at individual rain gauge locations following a cross-validation approach.

Fig. 7 .
Fig. 7. Snapshot images of the different spatial rainfall products at the time of peak areal intensity for Storms 1 (top) to 4 (bottom) over the Portobello catchment.From left to right: RD, BK, MFB, BAY and SIN3 (with singularity range [1, 3]) estimates.The black polygon indicates the boundary of the Portobello catchment, and the black and white markers respectively represent the location of flow and rain gauges.

Fig. 8 .
Fig. 8. Observed flows vs. simulated flows with RG, RD, BK, BAY and SIN3 rainfall inputs at selected flow gauging sites of Portobello catchment during Storm 3. Selected gauging sites: (top) FM3: upstream end of the catchment; (bottom) FM10: mid-stream area.The location of the selected monitoring sites is shown in Fig. 3.

Fig. 9 .
Fig. 9. Observed flows vs. simulated flows with SIN1-SIN5 rainfall inputs at selected flow gauging sites of Portobello catchment during Storm 3. Selected gauging sites: (top) FM14: upstream end of a small branch of the sewer system; (bottom) FM19: downstream end of the catchment.The location of the selected monitoring sites is shown in Fig. 3.

Fig. 10 .
Fig. 10.Boxplots displaying the distribution of regression coefficient (β) (left column), coefficient of determination (R 2 ) (middle column) and weighted coefficient of determination (R 2 w ) (right column) statistics derived from the linear regression analysis conducted for each pair of recorded and simulated flow time series at each gauging location.
is appropriate.Fig.5shows a further comparison of instantaneous areal 865average RG intensities vs. areal average BK, RD, BAY and SIN intensities throughout each of the storm events under consideration.As expected and in line with the analysis above, the areal BK estimates are generally in good agreement with areal RG estimates.Some underestimations can be 870 observed at high RG rainfall intensities; nonetheless, most of them are fairly minor and are still within a reasonable range which can be attributed to areal-point differences.With regard to the RD estimates, it can be seen that they tend to overestimate small rainfall intensities and underestimate the 875 peak intensities, with underestimation being more evident in the events with relatively high intensities (i.e.Storms 1, 3

Table 1 .
Selected rainfall events over the Portobello catchment.:Theaccumulation and peak intensity values shown in this table correspond to areal mean values for the entire domain under consideration (as shown in Fig.3b). Note

Table 2 .
Areal average rainfall accumulations and peak intensities for the different rainfall products.