Radar Estimation of Intense Rainfall Rates through Adaptive Calibration of the Z-R Relation

Rainfall intensity estimation from weather radar is still significantly uncertain, due to local anomalies, radar beam attenuation, inappropriate calibration of the radar reflectivity factor (Z) to rainfall rate (R) relationship, and sampling errors. The aim of this work is to revise the use of the power-law equation commonly adopted to relate radar reflectivity and rainfall rate to increase the estimation quality in the presence of intense rainfall rates. We introduce a quasi real-time procedure for an adaptive in space and time estimation of the Z-R relation. The procedure is applied in a comprehensive case study, which includes 16 severe rainfall events in the north-west of Italy. The study demonstrates that the technique outperforms the classical estimation methods for most of the analysed events. The determination coefficient improves by up to 30% and the bias values for stratiform events decreases by up to 80% of the values obtained with the classical, non-adaptive, Z-R relations. The proposed procedure therefore shows significant potential for operational uses.


Introduction
The use of weather radar to estimate rainfall rates and their spatio-temporal distribution is a well-established practice. Although the radar does not provide direct measures of rainfall fields, it is possible to estimate rainfall rates from the radar reflectivity measures. Weather radar provides spatial data with higher spatio-temporal resolution and lower maintenance costs than traditional rain gauging networks.
Rainfall rate and radar reflectivity are linked by the drop size distribution (DSD), as recognized by Marshall and Palmer [1]. The use of standard diameter distributions leads to a power law relation between the radar reflectivity factor Z (mm 6 /m 3 ) and the rainfall rate R (mm/h): where a and b are coefficients related to the DSD. In what follows, Equation (1) will be referred as the Z-R relationship and the "radar reflectivity factor" just as "reflectivity". As the Z-R relationship is inherently indeterminate, the coefficients a and b are usually estimated through an empirical approach, based on the comparison of radar and rain gauge data (e.g., [1,2]). The estimation procedure is characterized by a high level of uncertainty, as can be recognized by the wide variety of coefficient values collected in [3][4][5]. A variety of papers highlights the variability of a and b according to the type of event ( [6,7]) and, during the event itself, in time and space ( [8][9][10]). Even with measured drop size distributions, recorded over a physically homogeneous period, the instrumental noise related to small sampling effect and the observational noise, due to the drop sorting effect, can lead to Z-R relationships which are only true in a statistical sense [9].
Moreover, the expected estimation quality is considerably influenced by radar characteristics, sampling methodologies, and by the characteristics of the area in which the estimation is performed ( [11][12][13]). The validity of an approach based on the use of a constant Z-R relationship should therefore be strictly limited to the calibration domain of the relation itself, making the relation unsuitable for any systematic operational use. The importance of distinguishing the search for a physical dependence between reflectivity and rainfall intensity in a specific context from the goal of obtaining the best radar estimation of the rainfall rate is well discussed in Ciach and Krajewski [14].
Optimal radar estimation of the rainfall rate, the main objective of this study, has been addressed by several papers in the literature. Some authors [15][16][17] propose to reduce the uncertainty of estimation by introducing additional information, e.g., type of precipitation, distance from the radar, etc.
Generally, a significant improvement in the estimation quality is obtained by calibrating Equation (1) at the local scale, on the rain gauge nodes. This kind of approach reduces the problems related to orographic and climatological spatial variability. Two main categories of methodologies can be adopted: (i) A long-term (static) approach that involves the calculation of a single corrective term from an historical dataset of radar and rainfall measurements, to be subsequently applied to all the time frames in an identical way; and (ii) a short-term (dynamic) approach which requires the definition of a different calibration factor for each time step, based on the comparison between radar and rain gauge data from a definite time domain [18]. Several studies can be found describing the results of the second category of methods, with approaches suitable for different time scales: seasonal (e.g., [19][20][21]), monthly (e.g., [22][23][24]), or multi-daily (e.g., [25]).
To adapt the estimation to a single event, different methodologies have been proposed, that can be referred to as assimilation techniques of radar and gauge information to derive bias-corrected precipitation fields. Most of them stem from the pioneering work of Brandes [26] that involves the computation of a calibration matrix for each event, based on the smoothing of both the original radar field and of the subsequent correction factors. The key idea of this work has been seminal for a wide range of radar adjustment methods accounting for the spatial variability of the a coefficient of the Z-R relation in an operative context (see e.g., [27]). As a conceptual evolution of the idea proposed in [26] recent methods propose adjustments of the relation (1) with Z-R data collected from moving windows at a sub-daily scale (e.g., [28,29]). More recent studies also report on the combination of this kind of approach with appropriate interpolation techniques to significantly increase the estimation quality [30].
Some authors ( [31,32]) suggest another approach to increase the estimation quality, abandoning the classic radar-rain gauge comparison, and adopting different merging techniques, based on the use of geostatistical interpolators like kriging with external drift and co-kriging. More complex methodologies involve the use of advanced techniques such as statistical objective analysis [33] or Kalman filters (e.g., [34,35]). Most of these methods are very time consuming and often not suited to operational contexts [31].
Among the most recently proposed methods, the main advancement achieved in the operational efficiency relates to the use of non-linear calibration (i.e., methods which allow a and b to vary simultaneously). As demonstrated in various studies (e.g., [36]) those methods systematically outperform linear methods (where only the a coefficient is let to vary) due to the high variability of the b coefficient in time and space ( [9,37]).
While most of the studies cited use either a spatial or a temporal point of view, we address the problem from both viewpoints at the same time, through a non-linear calibration procedure amenable to explain the variability of both the a and b coefficients simultaneously in space and time. We propose a simple and operationally-ready methodology that produces an accurate radar-based estimation of rainfall intensity, that is adaptive in time and independent of local conditions. The quasi real-time calibration method proposed implies that Z-R relationships are consistent with the evolution of the event, while the use of a spatially adaptive approach avoids the use of interpolation techniques, making the technique amenable in large areas with complex orography. The aim is to make the methodology suitable for a systematic and unsupervised operational application.

The Adaptive in Time and Space Estimation Technique (ATS)
The proposed methodology aims to account for the spatio-temporal variability of the Z-R relation by means of an adaptive estimation of the two coefficients of Equation (1). The methodology, called Adaptive in Time and Space estimation technique (ATS) relies on the definition of calibration domains limited in time (Section 2.1) and space (Section 2.2).
In the following the rationale for the selection of the domains and the related calibration procedures in time and space are presented.

Definition of the Time Domain
To face the variability of the Z-R relation in time we propose to adopt a quasi real-time calibration window adapted from [29]: Coefficients a and b are estimated for each time step (t i ) considering the Z-R pairs belonging to a calibration window of duration d, [t i−d ; t i ], as shown in Figure 1. Such an approach allows one to follow the evolution of the event in quasi real-time, calibrating a different power-law relationship (1) at each time step. The inconvenience of this approach is that the high amount of noise can lead to highly variable estimates of the a and b coefficients due to the small size of the estimation domains.
To overcome this problem we propose a methodology involving a systematic "cleaning" of the data, using a reflectivity threshold variable in time to exclude from the estimation domain Z-R pairs with rainfall rate at or near zero. Working at a low rainfall rate scale, small differences between observed amounts and uncalibrated radar data can lead to spuriously large or small calibration factors [26]. Furthermore, the threshold is used to represents the amount of noise of the radar measurements due to instrumental and sampling uncertainly. The threshold is also used to discriminate between the presence and absence of rainfall at the considered time step: the rainfall rate is systematically set null for the locations characterized by under-threshold reflectivity values.
The threshold is calibrated for each t i by selecting the rain gauges that have registered no rainfall at the instant t i−1 . The threshold value is defined as the quantile (with cumulative probability q) of the empirical distribution of the reflectivity data associated to the selected rain gauges. The choice of the optimal q value is discussed in Section 2.3.
After some preliminary test, we found that for a 10 minute resolution, a time window of 60 minutes provides sufficient information for a robust Z-R calibration. We decided not to let the width of the time window vary, to reduce the degrees of freedom of the procedure.

Definition of the Spatial Domains and Parameters Estimation
To account for the spatial variability of reflectivity, the calibration domain is further confined, for each location p j , to the N nearest rain gauges, using only above-threshold reflectivity values. In other words, for each pixel of the gridded study area, a specific calibration domain of N station is defined, N being therefore an important pre-determinate parameter in the procedure. The selection of the optimal N is widely discussed in Section 2.3. If the valid number of Z-R pairs (i.e., the above-threshold pairs into the whole estimation domain) are lesser than N, we attempt an estimation with the available pairs. If the number of Z-R pairs is not sufficient for a robust estimation (i.e., if the system does not converge after 400 iterations), a backup regional relationship is used (see Section 3.2 and Appendix A). Considering that the threshold value and the spatial distribution of rainfall change continuously, the procedure requires to re-consider the number and position of valid stations for each location p j at each time step t i .
Minimizing the sum of the squared differences between the observed and the estimated rainfall for each local domain and in each time step we obtain a pair of optimal a and b coefficients. The estimated rainfallR is evaluated according to the following relationship, obtained from Equation (1), cleaning the recorded reflectivity with the threshold value Z th : whereâ andb are the estimated coefficients of the Z-R relationship, Z * (dBZ) and Z t h * (dBZ) the radar reflectivity and the reflectivity threshold respectively, both expressed in decibels (Z * = 10 · log 10 Z).
As the definition of the local domains refers to the nearest rain gauge, local domains referred to near locations partially overlap one another, granting a smooth variation of the coefficients and the continuity of the radar QPE. The estimation procedure is carried out in the Z-R plane by non-linear regression. The adopted optimization algorithm is a subspace trust region method and is based on the interior-reflective Newton method [38]. The optimized coefficients are estimated iteratively using as initial values the "static" coefficients calibrated at the regional scale (see Section 3.2). Calibration and verification procedures are carried out in cross-validation mode for each considered rain gauge, i.e., excluding one rain gauge at a time from the evaluation of the Z-R relationship and then comparing the estimated rainfall depth with the actual measurement at the excluded rain gauge.

Calibration Procedure
The methodology is characterized by two calibration parameters: N and q, indicating the number of rain gauges in the local domain and the probability used for the definition of the zero-rainfall threshold respectively. First, we select an adequate number of severe precipitation events for the calibration. The rationale for the selection of the events is described in Section 3.2. On this dataset, in order to define the best N and q parameter values we proceed by maximizing the estimation efficiency at both hourly and event scale, considering jointly the absolute estimation error ε abs , related to the hourly scale: and the BIAS, related to the event scale: where R i,j andR i,j are the rainfall observed and estimated using Equation (1) respectively at t i in p j ; n is the number of rain gauges and d the number of records per rain gauge.
The calibration of N and q is therefore carried out by exploring the parameter space 0 ≤ q ≤ 0.9 and 0 ≤ N ≤ 80 (q=0 corresponding to the application of the methodology without any threshold) and minimizing the performance index I 3 , defined, for each event (ev), as: where: and in which min(ε abs ) and min(BIAS) are the minima of ε abs and BIAS, evaluated in the q-N parameter space. Optimal q and N are then chosen for each event as the values corresponding to the I 3 index minima.

Case Study
The area of interest is located in North-Western Italy, almost entirely within the borders of the Piemonte region. The area is characterized by a complex orography as shown by Figure 2. The presence of orographic barriers importantly affect radar detection quality, due to the presence of echoes from the ground, occlusions of radar beam, orographic precipitation under beam, etc. The "Bric della Croce" radar is a dual-polarization C-band doppler with digital receiver located at 736 m above sea level on the hills near Torino (Geographical coordinates: 7.734E 45.035N). ARPA Piemonte (the Regional Agency for Environmental Protection) stores reflectivity maps on a Cartesian grid of 250 × 250 km, with a resolution of 500 m in space and 10 minutes in time. The raw radar product is de-cluttered with a fuzzy logic algorithm [39] and the horizontal corrected reflectivity is projected on a 2D map by considering, at each point, the radar beam with lowest elevation available. This information is obtained comparing the terrain elevation for each cell of the gridded domain (extracted from a 500 × 500 m Digital Terrain Model) with the radar beam elevation. From October to May a VPR (Vertical Profile of Reflectivity) correction is applied. The mean VPR is calculated by averaging the vertical profiles falling into a neighbourhood with 50 km radius from the radar site, during the previous 1-hour period.
ARPA Piemonte also manages the rain gauge network, which consists of 378 tipping bucket gauges distributed over the Piemonte region (25,000 km 2 ). Rain data have a 10 minutes resolution in time and a minimum threshold of rainfall detection of 0.2 mm/10 min. This entails that the minimum intensity that a rain gauge can perceive is 12 mm/h. Intensity measures lower than 12 mm/h, with temporal resolution of 1 min, are hence the result of post-analysis algorithms [40].
As the aim of this work is the definition of an adaptive calibration method particular attention should be paid to the quality of input data.
The magnitude of the rain gauge errors is highly dependent on the local rainfall intensity and the timescale. For example, at moderate rainfall intensities of 10 mm/h [41] report relative standard errors of 4.9% for the 5-min and 2.9% for 15-min timescales. The available 10-minutes time resolution therefore results in a good compromise between relatively low standard errors and a good robustness at the hourly scale.
Much has already been said on the quality of the radar data in the previous sections. Supplementary analysis of the radar product quality, referred to the study case, is performed in Section 3.2, after the Selection of the events.

Definition of the Set of Events and Estimation of the Static Coefficients
The criteria that lead to the definition of the set of rainfall events to be used in this analysis comply with the need of representing the most intense precipitation events recently occurred in Piemonte. Two different categories of phenomena can be defined: (i) Convective events, typically very localized and characterized by high intensities and short durations, and (ii) stratiform events, covering wider areas and characterized by lower intensities and longer durations. We select the most intense events occurred in the region from 2003 to 2008, with different criteria for convective and stratiform events, based on the duration and spatial distribution of rainfall fields. Convective events are identified as those registering the annual maximum value (of hourly duration) in, at least, 22 rain gauges. Stratiform events, instead, are selected when the extension of the area registering the highest annual areal precipitation in 24/48 hours equals, at least, 1/4 of the regional area.
The selected events (18) are listed in Table 1. In order to limit the amount of noise into the calibration procedure, we define the radar visibility map, to exclude the areas characterized by low quality reflectivity records, due to different source of errors (beam block, attenuation, etc.). To this aim we define the relation between the relative error at each gauged location and the beam height. A beam height of 4000 m a.s.l. is then used as threshold value to define the high visibility area, to be included in the calibration domain.
To further evaluate the quality of the used radar products we report in Table 2 the rate between invalid radar record at gauge location. Events 1 and 6 report anomalous high numbers of invalid radar records. We decided not to exclude the events from the case study to define the impact of these anomalies on the results.  (1) is then calibrated over the high-visibility area to define a regional static (in time and space) relationship. This relation can be used as a backup-equation when the dynamic approach cannot be used. The regional calibration procedure, reported in the Appendix A, leads to the definition of the relationship (A1). Moreover, during the regional calibration phase, we decided to exclude events 17 and 18 from the case study, due to the probable presence of snow that could lead to mis-calibrations over the no-snowing areas (see the Appendix A for details).
Estimates obtained for the selected events with the regional relationship (A1) systematically outperform those obtained with commonly used relations (e.g., [2]). However, limitations due to the use of a single relation on a vast and complex territory, leading, in some cases, to poor reconstruction of the precipitation volumes are yet to be overcome. Figure 3 shows an example of the presence of over/underestimation clusters as a function of the distance of the location from the radar. Figure 3. Comparison between observed precipitation R cum and rainfall estimated with the regional formulaR cum for the event occurred on 10/31-11/1/2003 at the event scale. (The grey scale refers to rain gauge distance from radar in km).

Calibration of the ATS Technique
As stated in Section 2.1, the width d of the calibration window is set to one hour, in order to limit the influence of the temporal variability in the adaptive search for optimal parameter values. The use of Z-R data recorded with a 10 minutes frequency grants a good estimation robustness (5 samples for each window, each characterized by a number of observations equal to the number of considered rain gauges).
The adaptive calibration is carried out by evaluating, for each event, the pattern of the I 3,ev index (Equation (5)) in the N-q parameter space. The minimization of this index allows to define the N-q pair that maximizes the estimation quality at both the hourly and the event scale.
The mean of the event patterns, each normalized to the maxima of I 3,ev for the event, is plotted in Figure 4. The minimum of the I 3 index in the N-q parameter space identifies the optimal pair of values for the whole set of events. To evaluate the presence of different behaviours among the different event typologies, the calibration procedure is carried out at first on the whole set (Figure 4a), then on the convective ones (Figure 4b) and finally on the stratiform ones (Figure 4c).
For convective events, due to the high spatial variability of rainfall fields, a rather high q threshold is required, in order to exclude "false positive" values. The number of rain gauges in the local domain (N) seems to play a less important role, if a minimum of 5 rain gauges (corresponding to 25 Z-R pairs) is ensured. Below this value, the quality of the estimations rapidly declines, as shown in the panel (b) of Figure 4.
For stratiform events the number of rain gauges to be included in the domain becomes crucial and should be appropriately set (usually between 10 and 30) in order to avoid including in the estimation procedure stations with low information content. In this case the choice of the quantile q is less influential (see panel (c)), due to the low variability of the reflectivity values.   We finally verify that the use of different N-q pairs for each category of events does not improve significantly the estimation quality. Therefore, in order to propose a methodology robust and easily applicable in real-time, we suggest to use as optimal values, both for convective and stratiform events, N = 20 and q = 85%, corresponding to the best global mean value of the I 3 index (Figure 4a).
These values are calibrated on the dataset under analysis but we are confident that they can be applied even in analogous spatial domains, whenever a local calibration can not be pursued. The limit of N = 20 in the spatial domain can grant a good robustness at the adopted time resolutions, and the relatively high quantile q for the definition of the threshold makes the methodology amenable to different kind of events. To support this statement, we show in Table 3 that the adopted quantile lets the threshold vary consistently even in our case study, without significant loss of performance with respect to the event-variable q values displayed by panels (b) and (c) of Figure 4.

Results and Discussion
The ATS procedure is applied to the 16 events listed in Table 1 using the q-N pair defined in the previous section. The performances of the procedure are compared with two static methodologies: the Joss-Waldvogel formula [2], routinely used for rainfall estimation over the study area, and the regional relation (Equation (A1)). We also compare our procedure with the adaptive methodology proposed by Brandes [26], that entails the estimation of a corrective factor at each rain gauge site, with a radar-gauge comparison carried out in real time. All factors are then interpolated on the whole radar field with the procedure described in [31]. Operatively, we adopt as first-attempt relation the regional one (Equation (A1)). To avoid the problems related to the low rainfall rate scale, we consider as valid only the pairs for which both recorded and estimated rainfall exceed 2.5 mm.
The quality indicators are defined at both the hourly and the event scale. Figure 5a,b show the coefficients of determination obtained with the four methodologies for all the events. Figure 5c shows the BIAS values.
The results confirm the validity of the proposed methodology that allows us to obtain lower volumetric errors for almost all of the events, with respect to both the static formulations. In addition, the correlation coefficients, not subjected to optimization during the calibration phase, show a general improvement in the estimation quality at both the hourly and the event scale.
As for the comparison with the adaptive methodology, the ATS technique shows remarkable improvements at the event scale, due to the contribution of the adaptive threshold to discriminate between presence and absence of rainfall, that leads to a more affordable estimation of the cumulative rainfall at the event scale. The improvement is greater for stratiform events, thanks to the spatial uniformity of these events, allowing an easier identification of the local spatial domains. For convective events the improvements are less marked, as the spatial variability of rainfall fields often leads the ATS technique to work as a regional estimator, similar to the one adopted in [26].
Events 1 and 6 are the only exceptions to the generally good performances. They are characterized by a generalized deterioration of the estimation quality, partly attributable to the poorer quality of the available data. The number of invalid radar record is indeed higher in events 1 and 6 than in all the other events, as reported in Table 2. This reduces the quality of results, as the efficiency of the proposed method, that involves a dynamic calibration, is quite sensitive to the quality of the input data.
To underline the validity of the proposed technique, in Figure 6 we compare the measured precipitation with the regional estimate (top graphs) and with the ATS estimate (bottom graphs) for 3 different events.
Graphs in column (a) and (b) show the comparison between the cumulative rainfall obtained with the regional relationship and the one obtained with the ATS technique, for the events 15 and 13, both considered at the event scale. The greater ability of the ATS technique to retrace the event is clearly apparent. The scatter plots in column (c) refer to the event 2, where data is aggregated at the hourly scale. In this case there is a deterioration of the determination coefficient, essentially due to the underestimation of a (limited) number of high values (that are underestimated also by the regional relationship).  Figure 5. Comparison between the coefficients of determination obtained with the Joss-Waldvogel formula, the regional formula, the methodology proposed in [26] and the ATS technique (a) at the event scale and (b) at the hourly scale. Correlation coefficient R 2 falling outside the range (0,1) have not been reported. (c) reports the comparison between the BIAS obtained with the four above-mentioned methodologies at the event scale.
The results demonstrate that the proposed technique is particularly suitable for stratiform events. The wide spatial scale allows for an easy identification of the local spatial domain. Figure 6b shows the clearly improved estimation quality obtainable with the ATS technique with respect to a unique Z-R relation for a large-scale event.
The increase in the estimation quality is less significant for convective events, where the localized nature of the high values of reflectivity signal contrasts with the uneven spatial distribution of the rain gauges, preventing the identification of a uniform and numerically robust spatial domain. In these cases, to reach a sufficient number of representative station our procedure would require to use distant Z-R pairs, that can be poorly representative of the event core. Moreover, due to the rapid evolution in time, an hourly calibration window can be still too wide, leading to the use of non-representative Z-R pairs. Given that the methodology has different performances for different spatial scales of the events, its validity for convective rainstorms needs to be further assessed. A more dense rain gauge network and the availability of radar-rainfall data with higher temporal resolution may facilitate the use of the adaptive approach, increasing the estimation quality, also for convective events. Figure 6. Comparison between observed precipitation and rainfall estimated with the regional formula (top graphs) and with the ATS technique (bottom graphs). The comparison is made at the event scale for event 15 (column a) and 13 (column b). Comparison in (column c) refers to event 3 at the hourly scale. (The grey scale refers to the distance from radar, in km).

Conclusions
A new, operationally oriented, procedure for the estimation of areal rainfall from weather radar is proposed. The methodology, called Adaptive in Time and Space (ATS), makes use of confined spatial and temporal domains for a quasi real-time calibration of the relation between radar reflectivity and rainfall rate at a local scale. By doing so we take into account the spatial and temporal variability of the Z-R relation, making the technique suitable for systematic operational use, regardless of local conditions, characteristics of the radar, sampling methodologies and spatio-temporal distribution of the events under analysis. The application of the ATS technique to a set of 16 severe rainfall events in the North-Western region of Italy provides a general improvement of rainfall estimation quality. A significant improvement in obtained particularly in the reconstruction of rainfall patterns of stratiform events from weather radar.
An accurate pre-processing of both radar and rainfall measurements is required in order to maximize the quality of the estimation. To make the methodology applicable in real-time, this operation has to be carried out in continuous time. This would require the implementation of automated procedures for evaluating the quality of input data and of output estimates, for immediate detection of anomalous values.
Further refinements to the proposed methodology are possible, such as varying the parameters of the calibration windows in time (e.g., time window amplitude, numbers of rain gauge into the local domain, percentile to be used to define the threshold), according to the evolution of the event. In order to make the procedure more robust and accurate, it could be also helpful to use multiple regression techniques during the calibration phase, by considering other radar variables (e.g., specific differential phase, differential reflectivity, etc.).
Future developments will address increasing the number of the events and the spatial domain, both to address the model performance deterioration for convective events and to test the validity of the ATS technique in different climatic areas, as well as to checking how different radar sampling methodologies (i.e., different beam height, different signal polarization, etc.) would affect the quality of the estimations.

Conflicts of Interest
The authors declare no conflict of interest.
A. Determination of a regional Z-R static relationship In this appendix we briefly describe how we calibrate a static relation valid for the whole regional territory. To this aim, Equation (1) is reconsidered on the whole acceptance area in order to obtain regional estimates of the a and b coefficients. The calibration procedure is targeted at minimizing the absolute estimation error ε abs (3) in the parameter space (1< a <1000, 1< b <4). In order to reduce the processing time, preliminary "data binning" was applied; i.e., radar reflectivity data were grouped into classes, with minimum width equal to the resolution of the radar data (0.5 dBZ). The median value was chosen as representative of each class. Classes with less than 10 items (i.e., the tails of the distribution) were assembled, to increase the robustness of the estimators. To each reflectivity class Z, a value of R equal to the average of the corresponding measured rainfall was associated. Once identified the pair of parameters (a,b) characterized by the minimum standard deviation min(ε abs ), the "sub-minimum area", i.e., the area where the absolute error is at most twice of min(ε abs ), is defined for each event ( Figure A1a). The a-b calibration space is then limited to the "sub-minimum area". On this domain we evaluate the BIAS for each a-b pair and for each event, according to Equation (4) ( Figure A1b).
In order to define the optimal parameter pair, the absolute error ε abs (3) and the BIAS (4) are combined for each event to define the index I 3,ev (5). This index takes into account the quality of the estimation both at the hourly and event-scale.
We obtain the optimal pairs of a and b values, reported in Figure A2, by minimizing I 3,ev for 8 convective and 10 stratiform events (see Table 1). Events 17 and 18, occurred in December 2003 and December 2008, are therefore excluded from the subsequent processing since they appear isolated from the others, probably due to the presence of snow. No distinction between convective and stratiform emerges from Figure A2.
The global optimal values of a and b are then obtained by applying the described procedure to the Z-R pairs of the whole set of events (i.e., merging convective and stratiform events).
The regional static relation reads: where the (a,b) pair is represented by the big black dot in Figure A2.  Figure A2. Optimal (a-b) values of the Z-R relation for each single event (grey symbols, according to event type), for the different event categories (black rhombus and cross) and for the whole set of events (black dot).