Comparison of Different Quantitative Precipitation Estimation Methods Based on a Severe Rainfall Event in Tuscany, Italy, November 2023

1. Introduction

Quantitative precipitation estimation (QPE) is a critical element in meteorology and hydrology, making accurate QPE methods essential for numerous applications, including flood forecasting. QPE methods typically employ data from two types of instrument: rain gauges and weather radars. Given the specific advantages and limitations of each of them, the information about rainfall provided by these instruments are somehow complementary. Therefore, their combination represents a solid ground for achieving accurate precipitation estimation.

Rain gauges measure precipitation directly, providing highly accurate data on the amount of rainfall at a specific locations. Given their precision, these instruments are used for calibrating and validating other precipitation estimation methods. However, rain gauge measurements may be affected by errors, both systematic and random, depending on the atmospheric conditions and the type of rain gauge [1,2]. Taking into account the point-wise nature of rain gauges, the rainfall measurements need to be spatially interpolated in order to generate a 2D rainfall map; for this purpose, deterministic (e.g., inverse distance weighting [3]) or geostatistical (e.g., kriging [4]) techniques are typically used (see [5,6] for a comparison). The nature of these instruments represents also their main limitation: in fact, a dense gauges network is required to capture the spatial variability of the rainfall patterns, especially in case of severe small-scale phenomena.

Weather radars, on the contrary, can monitor precipitation over large areas with high temporal and spatial resolution. Radars emit microwave pulses and measure the power backscattered from ensemble of hydrometeors, providing a practically continuous monitoring of the spatial distribution of rainfall an of its intensity [7]. Despite their broad coverage and excellent resolution, both in time and space, radar-based QPE is subject to several sources of error, including signal attenuation due to rainfall itself, ground clutter and beam blockage [8]. Moreover, radars do not provide direct rainfall measures, thus a conversion is required to obtain the QPE. This conversion is typically made resorting to the so-called Z–R relationship [9]:

$$Z=a{R}^b,$$

(1)

which relates the reflectivity factor Z (expressed in ${\mathrm{mm}}^6·$${\mathrm{m}}^{-3}$), a parameter derived from the backscattered power measured by the radar, and the rainfall rate R (expressed in mm·${\mathrm{h}}^{-1}$). Unfortunately, the choice of the coefficients a and b is critical. In fact, both Z and R depend on the drop size distribution (DSD), which is influenced by the type of event, the geographical area and the climatic conditions [8,10,11]. As a result, a and b are highly variable [12], and using the Z-R relation with a priori-fixed coefficients can lead to significant rainfall estimation errors [13].

Due to the complementary features and issues of rain gauges and weather radars, integrating their data can be usefully exploited to mitigate their limitations and to enhance the QPE accuracy. For this reason, merging rain gauge and radar data has been investigated for several years and various methods, with different approaches, have been developed (see for instance the review in [14]). One criterion with which the different approaches can be categorized is based on how the two data types are employed during the fusion procedure [14]. Methods belonging to the geostatistical interpolation category use the rainfall point measurements provided by rain gauges as primary quantitative information source, and the radar rainfall field as auxiliary spatial information to improve the accuracy of the estimation. These methods derive mostly from the kriging technique and the most employed is the kriging with external drift (KED), in which the external auxiliary variable (i.e, the radar-based QPE) is incorporated in the interpolation process of the primary variable (i.e., rainfall data provided by rain gauges) [4]. This category includes also the conditional merging (CM), in which the kriging-based QPE is corrected by adding an error obtained from radar rainfall observations [15]. Different merging approaches consist in adjusting the radar-based QPE by exploiting the rain gauge measurements. In this kind of methods – typically called radar bias adjustment methods because they address the issue of systematic bias in radar data – a multiplicative or additive correction factor obtained by comparing radar estimates with rain gauge measures is applied to the radar rainfall field. The simplest method in this category is the mean field bias (MFB) adjustment, in which a spatially constant multiplicative correction factor is used [9]. Other adjustment methods, such as that proposed by Brandes (Brandes spatial adjustment, BSA) [16], consider a correction factor that varies in space, so as to take into account of nonuniform bias across the radar coverage area. Finally, merging methods employing a space-time varying Z–R relationship have been proposed, such as that presented in [17] that is based on the dynamic adaptation of the coefficients a and b. This method, named space-time adaptive coefficient conversion (STACC), has proven its effectiveness for QPE especially when rainfall events with localized intense peaks affect areas characterized by a sparse rain gauge network [18].

In this study, we present an evaluation of different QPE methods based on data collected during a severe rainfall event in Tuscany, Italy. Methods relying solely on rain gauge data, solely on radar data, and on their integration, are considered. The paper is structured as follows: in Section 2, we first describe the event examined in this study and the employed radar and rain gauge datasets; then, we present the different QPE methods considered, as well as the procedure adopted to evaluate and compare them. The experimental results are shown and discussed in Section 3. Finally, in Section 4, we draw some concluding remarks.

2. Materials and Methods

2.1. Meteorological Event Description

The information that follows, concerning the meteorological phenomenon considered in this study, comes from the meteorological report provided by the LaMMA (Environmental Modeling and Monitoring Laboratory for Sustainable Development) Consortium [19], to which the reader is referred for further details.

In the afternoon of November 2, 2023, in Tuscany, the movement of a relatively stationary instability line (Figure 1) ahead of a cold front facilitated the development of several thunderstorms. These storms were intensified by an energetically favorable atmosphere (Figure 2) and highly humid low-level conditions, leading to an intense localized precipitation. Specifically, starting at 15:30 UTC, the convergence of strong winds in the southern parts of the region with the gradual influx of winds affecting the northern parts, combined with the advance of cooler air aloft and favorable conditions in the lower levels, led to the formation of a thunderstorm line. This system remained nearly stationary until 20:30 UTC, extending from the provinces of Livorno and Pisa to the provinces of Pistoia, Prato, and Florence. The location, intensity, and stationary nature of this thunderstorm system exhibited nonlinear characteristics and presented significant forecasting challenges. As a result of the heavy precipitation, some rivers and creeks flooded, causing extensive damages to the environment and affecting local communities.

2.2. Datasets

2.2.1. Radar Dataset

The weather radar data used in this study were collected by the Italian National weather radar network. Managed by the Italian National Civil Protection Department (DPCN), the network consists of 26 radar systems (22 C-band radars and 4 X-band radars, most of them dual-polarized) distributed to cover the entire national territory. Figure 3 shows the radar systems of the network that cover Tuscany. Specifically, these systems are C-band dual-polarized radars located in Monte Crocione (Tuscany, lat: ${43.96}^{\circ }$, lon: ${10.61}^{\circ }$, height: 1017 m a.s.l.), Monte Settepani (Liguria, lat: ${44.25}^{\circ }$, lon: ${8.20}^{\circ }$, height: 1390 m a.s.l.), Gattatico (Emilia Romagna, lat: ${44.79}^{\circ }$, lon: ${10.50}^{\circ }$, height: 34 m a.s.l.), San Pietro Capofiume (Emilia Romagna, lat: ${44.66}^{\circ }$, lon: ${11.62}^{\circ }$, height: 10 m a.s.l.), and Monte Serano (Umbria, lat: ${42.86}^{\circ }$, lon: ${12.79}^{\circ }$, height: 1428 m a.s.l.).

The radars of the network perform synchronous scans every 5 minutes. The raw data from all systems are processed by the DPCN (procedures are described in [20]), which afterward provides national-scale composite products. The radar reflectivity composite products include Constant Altitude Plan Position Indicator (CAPPI) – namely, the reflectivity at a specific altitude (from 1000 m up to 7000 m) obtained by interpolating data from scans made at multiple elevation angles – and Vertical Maximum Intensity (VMI), that is the maximum reflectivity value detected vertically above each ground pixel.

When facing the problem of merging radar and rain gauge data for QPE, the height of radar observations must be carefully considered, as it has an impact on estimation accuracy. In fact, several factors are involved in the choice of the most appropriate height for the reflectivity product, including local orography and topography, which may cause beam blockage and ground clutter, as well as vertical variations in precipitation resulting from evaporation, wind shear, and changes in hydrometeor type [8]. Considering the specific situation of Tuscany, we utilized the CAPPI radar reflectivity at an altitude of 2000 m to balance the various issues and reach a good compromise. Figure 4 shows the chosen product, which is available with a time step of 10 minutes on November 2, 2023. The 2D maps of reflectivity Z have a spatial resolution of 1×1 ${\mathrm{km}}^2$ (when mapped onto the x-y Cartesian coordinate system).

2.2.2. Rain Gauge Dataset

In this study, we used rainfall data from the rain gauge network of Tuscany, which is illustrated in Figure 5. The network is composed of 425 tipping-bucket rain gauges located across the entire 22,987 ${\mathrm{km}}^2$ area of Tuscany, which means approximately one gauge every 54 ${\mathrm{km}}^2$. Near real-time measurements of accumulated rainfall at 15-minute intervals, recorded by all available stations, are made accessible online through a dedicated platform [21].

By means of the rain gauge data, we defined the “study area” through the following procedure. First, we selected the time interval from 15:30 UTC to 20:30 UTC (i.e., when the considered rainfall event was most intense). Next, we calculated the convex hull encompassing all rain gauges that recorded an average rainfall rate of at least 1 mm·${\mathrm{h}}^{-1}$ during that interval. The study area was then defined as the region enclosed by the convex hull also shown in Figure 5.

2.3. QPE Methods

2.3.1. Radar-Based Method

Estimating rainfall by using solely weather radar data is typically done through the empirical relationship expressed in (1). In this work, we computed the rainfall field (in mm) as follows: first, the available 2D Z maps were averaged in space over a 3×3 window; then, they were converted into R maps with a fixed relationship; finally, the R maps were accumulated in time.

The relation employed for the conversion is $Z=300R^{1.4}$ ($R=0.017Z^{0.714}$), which is utilized in the QPE algorithms of the Multi-Radar Multi-Sensor (MRMS) system for convective precipitation [22]. This relation was preferred to the classic Marshall-Palmer relation ($Z=200R^{1.6}$) [23], which is more appropriate for stratiform phenomena. This QPE method is referred to as ZR in the following.

2.3.2. Rain Gauge-Based Method

Creating a 2D rainfall map from rain gauge data involves using spatial interpolation techniques to estimate rainfall values at locations which are not directly covered by the rain gauge network. In this regard, we employed the ordinary kriging (OK) technique, which aims at producing an estimate at a given location by considering the values of nearby observations and their spatial distribution. The estimate is a weighted average of the known values, with weights derived from the spatial correlation pattern to minimize the variance of the estimation error (see [4] for a detailed description of OK).

A key component in kriging is the semivariogram, which quantifies the spatial autocorrelation of the data. Thus, it is necessary to calculate the empirical semivariogram from the observed rainfall values and then fit a theoretical model (e.g., spherical, exponential, Gaussian). In this study, we assumed a spherical model and used the open-source library GSTools [24] for semivariogram computation and kriging implementation.

2.3.3. Kriging With External Drift (KED)

KED is a fusion method for estimating rainfall that enhances traditional kriging by integrating external information (drift) provided by radar into the spatial interpolation process of rain gauge data (see [4] for an exhaustive explanation about how KED works). As to the implementation of KED, in this study we used the GSTools library setting a spherical semivariogram model. As external drift variable, we considered the radar-based QPE obtained by means of the procedure in Section 2.3.1.

2.3.4. Conditional Merging (CM)

Kriging with radar-based error correction, also known as conditional merging (CM), is a technique that refines the kriging-based rainfall field by exploiting the superior spatial information provided by radars [15]. In CM, the rainfall field obtained by using kriging interpolation of rain gauge measures is corrected by adding an error field derived from radar data. The error field is calculated as follow: first, the radar-based rainfall field is estimated; next, the radar rainfall estimates at each rain gauge location are interpolated through kriging; finally, the resulting interpolated rainfall field is subtracted from the original radar-based rainfall field to derive the error field. In this work, we utilized OK (Section 2.3.2) and the rainfall map derived from ZR (Section 2.3.1) as the kriging technique and the radar-based rainfall field, respectively.

2.3.5. Mean Field Bias Adjustment (MFB)

MFB is a radar bias adjustment method that consists in applying a space-constant multiplicative adjustment factor to the radar-based QPE. The adjustment factor (AF) is derived from the rainfall estimated by the radar at rain gauge locations and the rainfall measured by the gauges [9], and is given by:

$$\mathrm{AF}={\displaystyle \frac{{\sum }_{i=1}^N{G}_i}{{\sum }_{i=1}^N{R}_i}},$$

(2)

where N is the number of available rain gauges, $G_i$ is the rainfall measured by ith rain gauge and $R_i$ is the rainfall estimated by radar at the ith rain gauge location.

For the implementation of MFB, we included in (2) only those rain gauges for which both the measured rainfall and the corresponding radar estimate exceeded 1 mm. The resulting adjustment factor was then applied to the rainfall map estimated through ZR (Section 2.3.1).

2.3.6. Brandes Spatial Adjustment (BSA)

BSA is an alternative radar bias adjustment method that, rather than applying a constant adjustment factor across space, uses a spatially varying correction. In the BSA method, a local adjustment factor is first computed at each rain gauge location by comparing the rainfall estimated by the radar with the rainfall measured by the gauge. These local adjustment factors are then spatially interpolated throughout the area of interest [16]. Specifically, to calculate the adjustment factor, we followed the procedure detailed in [25], excluding rain gauges for which either the measured rainfall or the corresponding radar estimate was less than 1 mm.

2.3.7. Space-Time Conversion Coefficients (STACC)

STACC is a rain gauge and radar data fusion method designed to employ a Z–R relationship that varies both in time and space. This approach dynamically adapts the coefficients in (1) according to the spatial and temporal evolution of the observed phenomenon, ensuring a more accurate QPE compared to using a constant a priori-selected relation. In the STACC method, the relationship used to convert the space-time average of Z into the time average of R is obtained through a two-step procedure consisting in local calibration and spatial interpolation. During the local calibration step, the optimal relationship coefficients are determined at each rain gauge location by applying a linear regression technique to the averaged values of rain gauge and radar data. In the successive step, these optimal coefficients are spatially interpolated across the study area. A more in-depth description of the STACC method can be found in [17].

As concerns the implementation of the STACC, we utilized a weighted least square regression methods that accounts for the uncertainties of the two variables [26]. In particular, we set the temporal variability (i.e., the temporal standard deviation) of rain gauge and radar data as the uncertainties. The local coefficients obtained were assumed to be acceptable and employed in the subsequent spatial interpolation step if their values fell within an established range (typically found in the literature [27], $a\in [30,1000]$ and $b\in [0.8,2]$). During the spatial interpolation step, the inverse distance weighting technique [28] was used.

2.4. Performance Analysis

Verifying the accuracy of QPE methods is challenging due to the unknown nature of the true rainfall field. Typically, this challenge is addressed by exploiting the high accuracy of rain gauge measurements, comparing the estimated rainfall values to those observed. Due to unavailability of an independent gauge network for the verification (i.e., a network different from that employed for the estimation), a leave-one-out cross-validation (LOOCV) is carried out to analyse the performance of the various QPE methods. The LOOCV procedure consists of the following steps:

• Remove one rain gauge from the dataset.
• Use a QPE method to estimate the rainfall at the location of the removed rain gauge.
• Compare the estimated rainfall y to the actual value x measured by the gauge.
• Repeat the above steps for all the available rain gauges.

To ensure a comprehensive analysis, this procedure is repeated throughout a designated observation period. Let T be the total observation time, during which the analysis is performed; assume that T is divided into several time slots of duration $T_s$, where $T_s$ is a submultiple of T, yielding a total of $N_t=T/T_s$ time slots. For each time slot, the previously described LOOCV is applied to the $N_g$ available rain gauges. Therefore, for each ith rain gauge, $i=1,2,\dots ,N_g$, and for each jth time slot, $j=1,2,\dots ,N_t$, we denote with $x_{i,j}$ the actual rainfall amount measured by the ith rain gauge during the jth time slot and with $y_{i,j}$ its estimate obtained with a QPE method by excluding the ith rain gauge from the dataset.

Finally, the accuracy of each QPE method is assessed by calculating the mean absolute error (MAE), the root mean square error (RMSE), and the coefficient of determination ($r^2$), defined as:

$$\mathrm{MAE}={\displaystyle \frac{1}{N_g·N_t}}\sum _{i=1}^{N_g}\sum _{j=1}^{N_t}\left|y_{i,j}-x_{i,j}\right|,$$

(3)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N_g·N_t}}\sum _{i=1}^{N_g}\sum _{j=1}^{N_t}{(y_{i,j}-x_{i,j})}^2},$$

(4)

$$r^2=1-{\displaystyle \frac{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_t}{(x_{i,j}-y_{i,j})}^2}{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_t}{(x_{i,j}-\overline{x})}^2}},$$

(5)

where $\overline{x}$ in (5) is the arithmetic mean of $x_{i,j}$ ($i=1,2,\dots ,N_g$, $j=1,2,\dots ,N_t$).

Another important aspect we would like to emphasize in this study is evaluating the different QPE methods in critical scenarios, such as when the rain gauge network fails to detect areas of very intense rainfall [18]. This critical situation may arise when a highly intense and localized precipitation event occurs, and the gauge network is not adequately dense to capture it effectively. Certainly, in such cases, meteorological radars are crucial for identifying high-risk areas. However, from a quantitative perspective, the lack of rain gauge data significantly impacts on the accuracy of rainfall estimates. In this study, the robustness of the various QPE methods in critical scenarios is assessed by applying the LOOCV approach to selected rain gauges in specific areas where precipitation is intense.

3. Experimental Results and Discussion

3.1. Leave-One-Out Cross-Validation (LOOCV)

The LOOCV procedure detailed in Section 2.4 was applied to data collected between 15:30 UTC and 20:30 UTC, spanning a total duration of $T=5$ h. This period was chosen as it corresponds to the highest intensity interval of the investigated event. Rain gauges placed on the vertices of the convex hull that defines the study area (see Section 2.2.2) were never removed in the LOOCV procedure. In the case of kriging-based methods (i.e., OK, KED, and CM), data from only the 50 nearest gauges were used to estimate rainfall at the removed rain gauge location to reduce computational cost. In the following, we differentiate the QPE methods according to the these categories: methods that rely on a single data source, either radar and rain gauges (i.e., ZR and OK); geostatistical fusion methods, that use rain gauge data as primary source of information and are based on geostatistical interpolation (i.e., KED and CM); non-geostatistical fusion methods, that use radar data as a primary source of information, by operating either a radar bias adjustment (MFB and BSA) or an adaptive use of the Z-R relationship (STACC). For the sake of simplicity, in some cases the last two categories are jointly referred to as fusion-based methods.

The results obtained for $T_s=30$ min are presented in Table 1. Among the methods that rely on a single data source, OK outperforms ZR according to all performance metrics. The good performance of OK can be attributed to the high density of the rain gauges network (see Section 2.2.2). This fact positively impacts also the results of the geostatistical fusion methods, with KED achieving the best overall performance among all considered methods. As to the remainder fusion-based methods (i.e., non-geostatistical), MFB demonstrates the worst performance (providing just slightly better metrics than ZR), whereas STACC produces the best results.

Since we are mostly interested to areas with significant levels of rainfall, the LOOCV procedure was repeated choosing exclusively those rain gauges that recorded an average rainfall rate greater than or equal to a specified threshold $\gamma$ over the interval T. Four different thresholds were selected: ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$, and ${\gamma }_4$, set to 5, 10, 15, and 20 mm·${\mathrm{h}}^{-1}$, respectively. Table 2 shows the values of the three performance metrics obtained for each method and threshold, while Figure 6 shows, for each method, the trends of these metrics vs. $\gamma$. From these results, we observe that ZR remains the least effective among all QPE methods, especially at the highest values of $\gamma$, for which the $r^2$ value drops to zero. Regarding the fusion-based methods, MFB again provides the poorest performance (with $r^2$ approaching zero), while KED continues to demonstrate the best results. As the threshold level increases, the improvement of the fusion-based methods over single-source ones becomes more pronounced, especially in terms of RMSE and $r^2$.

3.2. Analysis with Missing Data: Case 1

Figure 7 shows the radar reflectivity map, space-time averaged from 15:30 UTC to 16:30 UTC. The rain gauges labeled with numbers represent the ones that measured the highest rainfall levels during this time interval, with rain gauge #525 recording the greatest amount (56.80 mm); it is evident that these rain gauges are also located, as expected, in high radar reflectivity areas. The maps illustrated in Figure 8 display the rainfall accumulated over the same time interval, as estimated through the various QPE methods. In the maps obtained by MFB, BSA, and STACC the rain gauges that were actually used, i.e. those considered acceptable (see Section 2.3.5 and Section 2.3.7), are indicated by blue crosses.

To simulate and analyse a critical scenario (Section 2.4), we removed rain gauge #525 and estimated the rainfall y at its location by using the different QPE methods. We then compared these estimates to the actual measured value and computed the estimation error e to evaluate the accuracy of each method. The estimation errors are reported in Table 3. All the QPE methods underestimate the actual amount of rainfall, with the single-source estimation methods (OK and ZR) yielding, as expected, the largest errors (65.39% and 67.78%, respectively). The poor performance of these methods is inherited by methods that use geostatistical interpolation (KED and CM): in this case, the additional spatial information coming from radar is insufficient to compensate for the inaccuracy coming from rain gauges.

The ZR underestimation affects also the radar bias adjustment methods, i.e. MFB and BSA. For instance, in the BSA case, a large adjustment factor would be needed in the position of raingauge #525 to compensate for the radar underestimate. Instead, one gets the interpolation of smaller correction factors computed at the surrounding rain gauges. Similar considerations can be made for MFB.

The smallest error (27.73%) is achieved by the STACC method. Probably, the reason lies in the fusion approach of the STACC. When rain gauge #525 is removed, the Z–R relation coefficients at its location are determined by spatial interpolation of the coefficients calibrated at the surrounding rain gauge positions. Since the calibration is performed using linear regression, the resulting relationship, although based on smaller values of Z and R, can still be effectively applied to higher values of these quantities.

The simulation was repeated by removing also the other labeled rain gauges (i.e., #422, #829, #482, #1873) to further amplify the impact of missing data in the analysed high rainfall area. As shown in Table 4, removing these additional rain gauges negatively impacts all methods except for ZR, which depends only on radar data, and STACC. As expected, OK is significantly impacted, resulting in an error of 90.00 %, which consequently affects KED and CM results. MFB and BSA provide a slightly worse estimates, since, due to the aforementioned reasons, the adjustment factor is further underestimated. On the contrary, despite the removal of additional rain gauges, the STACC estimate error slightly decreases, confirming the low sensitivity of this fusion approach to the absence of rain gauges in areas featuring high levels of precipitation.

For the sake of completeness, the accumulated rainfall maps obtained by removing the labeled rain gauges are shown in Figure 9 (the map obtained by means of the ZR method is not included because it does not differ from the one shown in Figure 8). By comparing Figure 8 with Figure 9, it can be observed how the removal of the rain gauges affects the rainfall estimation results, both in terms of quantitative precision and spatial pattern.

3.3. Analysis with Missing Data: Case 2

Figure 10 shows the space-time averaged radar reflectivity from 16:30 UTC to 17:30 UTC. The rain gauges labeled with numbers represent those that measured the highest rainfall levels during this time interval, with rain gauges #3980 and #1835 – which are very close to each other – recording the greatest amount (both 44.00 mm). These rain gauges are situated in an area with very high reflectivity, the highest on the map (about 50 dBZ). The maps illustrated in Figure 11 display the rainfall accumulated over the time interval, as estimated by the various QPE methods. As in Figure 8 and Figure 9, the blue crosses in the maps related to MFB, BSA, and STACC indicate the acceptable rain gauges that were actually used.

In a similar way as what presented in Section 3.2, we simulated a critical scenario by removing rain gauges #3980 and #1835. The performance of each QPE method was then assessed by calculating the rainfall amount y at their locations and the error e between these estimates and the actual values. Table 5 presents the estimation errors obtained from these simulations. In this case, OK again tends to underestimate the actual level of rainfall for both rain gauges, but, unlike the previous simulation, the ZR method leads to an overestimate due to the very high radar reflectivity associated with the rain gauge locations and the chosen Z–R relationship; this influences on the results relative to KED, CM, MFB, and BSA. In KED and CM, the additional spatial information provided by the radar helps to correct the underestimation from the rain gauges, leading to very good results for KED, and results for CM similar to ZR. Concerning MFB and STACC, the overestimation of the ZR method is further amplified by the application of the adjustment factor, resulting in high errors (more than 100% for MFB). As in the previous simulation, the STACC method yields the best results, with errors of 5.52% and 1.91% for rain gauges #1835 and #3980, respectively, demonstrating again to be the most robust among all the analysed methods.

We also repeated the simulation by removing rain gauge #1830 as well, thereby increasing the effect of missing measurements. The results in Table 6 show that all the methods, except for ZR, lead to large errors, intensifying either the underestimation (for OK) or the overestimation (for all the other methods). In particular, KED seems to be the most affected method by the removal of rain gauge #1830, showing the largest increase in terms of relative error (e.g., from 11.55% to 97.80% for rain gauge #1835). Despite a loss of performance, the STACC method still provides good results and confirms to be the most robust method in this simulation.

As we did for the previous simulation, we derived the accumulated rainfall maps after removing the specified rain gauges. The results are shown in Figure 12.

4. Conclusions

In this paper, a comprehensive comparison of different QPE methods based on rain gauge and weather radar data has been presented. To this aim, data collected during a severe rainfall event occurred in Tuscany, Italy were employed.

We performed a LOOCV to assess the accuracy and the reliability of each considered method. The results have shown that relying exclusively on radar data leads to quantitatively inaccurate predictions, as expected. Additionally, using the MFB method – which integrates rain gauge data by simply applying a constant adjustment factor to the radar-based rainfall field – does not significantly enhance the estimation accuracy. The results have also demonstrated that a sufficiently dense rain gauge network, such as the one in Tuscany, can produce good results even when it is utilized alone, without the radar data, as is done in the OK method. In general, however, fusion methods, like KED, CM, BSA, and STACC, allows better results than OK (with reference to the used performance metrics MAE, RMSE, and $r^2$) to be achieved, especially when considering rain gauges that recorded a certain level of rainfall. Among these methods, KED, which belongs to the geostatistical interpolation fusion method category and thus benefits from the dense rain gauge network, has proven to be the most effective QPE method in this context, providing the best performance in terms of all the metrics.

A second type of analysis was also carried out, trying to highlight the robustness of the QPE methods in critical situations, in which high rainfall areas might occur in regions low-populated in terms of rain gauges. This analysis was performed by removing some rain gauges from the data available to the QPE algorithms and evaluating the estimation. This analysis disclosed, as expected, the main drawback of methods relying solely on rain gauge data and pointed out the essential role of spatial information provided by radar. Among the fusion methods, STACC demonstrated the highest robustness under these challenging conditions, yielding the smallest estimation errors.

In general, the study has highlighted the advantages of merging rain gauge data with weather radar observations to improve QPE. In fact, through the integration of these two data sources, we have achieved more accurate rainfall estimates compared to using either data source independently. The findings have offered further evidence of the validity of fusion techniques for both operational meteorology and hydrological applications, such as flood forecasting, where precise precipitation measurements are essential.