Assimilating X- and S-Band Radar Data for a Heavy Precipitation Event in Italy

Abstract

During the night between 9 and 10 September 2017, multiple flash floods associated with a heavy-precipitation event affected the town of Livorno, located in Tuscany, Italy. Accumulated precipitation exceeding 200 mm in two hours was recorded. This rainfall intensity is associated with a return period of higher than 200 years. As a consequence, all the largest streams of the Livorno municipality flooded several areas of the town. We used the limited-area weather research and forecasting (WRF) model, in a convection-permitting setup, to reconstruct the extreme event leading to the flash floods. We evaluated possible forecasting improvements emerging from the assimilation of local ground stations and X- and S-band radar data into the WRF, using the configuration operational at the meteorological center of Tuscany region (LaMMA) at the time of the event. Simulations were verified against weather station observations, through an innovative method aimed at disentangling the positioning and intensity errors of precipitation forecasts. A more accurate description of the low-level flows and a better assessment of the atmospheric water vapor field showed how the assimilation of radar data can improve quantitative precipitation forecasts.

1. Introduction

Almost every year, during the fall, heavy precipitation events (HPEs) bring destruction and cause fatalities somewhere in the Western Mediterranean (WM) region [1,2,3,4]. HPEs can either persist for several days in large areas, resulting in extensive flooding (e.g., 1966 Arno, Italy [5,6] and 1994 Piedmont, Italy [7,8]), or manifest at the sub-daily scale and produce flash floods (e.g., 1999 Aude, France [9,10] and 2011 Liguria, Italy [11,12]).

Short HPEs often assume the form of V-shaped, quasi-stationary convective systems [13,14]. These storms are characterized by a back-building dynamics, where convective cells developing upstream of the affected area continuously replace dissipating cells, resulting in pulsating heavy rain [15,16]. This continuous cell replacement is frequently initiated offshore [17]. This structure becomes quasi-stationary due to the perduration of a favorable synoptic and mesoscale environment [18]. Important factors for the persistence and evolution of these precipitation phenomena can be found at the surface (e.g., moist and conditionally unstable air associated with warm sea surface temperatures [19,20]), within the boundary layer (e.g., convergence associated with a frontal system [21], induced by orography [11,13,20] or land-sea differences [19,22]), and throughout the troposphere (e.g., coupling between low and upper-level jet streaks [23]). When these conditions persist for several hours, the system may assume the characteristics of a mesoscale convective system, which is able to partially drive the mesoscale circulation [24,25,26].

The frequency of these HPEs over the WM region has already increased [27,28], although a significant reduction of annual rainfall has been observed in relatively recent years [29,30]. Climate projections support this trend, by predicting a further increase in the frequency of HPEs, and a decrease in average precipitation in the WM region [31]. In this framework, considering that in the period 1975–2001 the average mortality per flash flood was 5.6% of the affected population [32], the most threatening consequence of the increasing HPE frequency is the expected increase of human losses associated with the flash floods they produce. Monitoring, mitigating, predicting, and communicating the potential of such dramatic events are necessary for reducing mortality [33].

Several studies investigated improvements in the characterization of severe precipitation events arising from the assimilation of radar data. Some of these studies [34,35,36] compared the performances of ensemble-based or hybrid assimilation techniques [37,38] and concluded that they provide better analyses and forecasts than traditional three-dimensional variational (3D-Var) data assimilation. Other studies [39,40,41,42,43,44] demonstrated the effectiveness of the 3D-Var method in assimilating radar data while still having a relatively low computational cost. The vast majority of these studies used weather analyses as initial and boundary conditions and assimilated data which were measured during the events. This methodology, despite being optimal for accurate reconstruction of HPEs, is not suitable for evaluating operational improvements. For bridging the gap between the improvements described in previous works and operational improvements, a reasonable latency period consisting of observation collection, data assimilation, WRF run-time, and forecast communication for early warnings should be taken into account. Considering this latency by assimilating observations before—and not during—the storm allows for a quantification of the impact data assimilation has on the predictability of HPEs in a pseudo-operational framework.

This work focuses on a numerical reconstruction of the quasi-stationary convective system that, during the night between 9 and 10 September 2017, caused several flash floods in the town of Livorno, Italy, resulting in nine fatalities. The Livorno case was recently studied by [45] using the non-hydrostatic model RAMS@ISAC [46]. The authors evaluated the impact of assimilating lightning and radar reflectivity data in the short term (i.e., forecast length shorter than 3 h) and at the convection-permitting scale (i.e., grid spacing up to about 1 km). The authors underlined the paramount importance of the reflectivity data and found that the assimilated runs outperformed the control one (i.e., no data assimilated), by reducing the number of missed precipitation events at the expense of a higher number of false alarms. The authors acknowledged that not updating the initial and boundary conditions in the assimilation step represents a limit of their work, and they claimed that exploring this issue deserves further investigations. In a recent paper, Lagasio et al. [47] tested the effect of assimilating a broad range of remote sensed data into the WRF model. The authors found that information about the wind field, and to a minor extent, atmospheric water vapor content, is crucial for better localization of the rainfall peaks.

We present a simulation of the Livorno case using the WRF model in an operational-like configuration. We investigated the impact of assimilating radar and automatic weather station data for a refinement of the operational setup, which should lead to immediate benefits for any early warning system that relies on data produced by numerical weather models. In particular, the main goal of the work was to evaluate to what extent the information carried by a relatively small radar system, such as the X-band radar located at the Livorno harbor [48], can improve the predictions of the WRF model in the short term. In fact, in recent years, some European projects, such as those funded in the framework of the European Cross-Border Cooperation Programme Italy-France “Maritime,” dealt with the deployment of a network of X-band radars in the Tuscany region and surrounding areas. The explicit goal of such projects is to monitor potential severe weather occurring offshore, using radar instruments, an essential tool for nowcasting applications. The cross-border sharing of such relevant meteorological observations and the integration with existing tools and methodologies are intended to improve the forecasting and alerting capabilities of operational regional weather agencies. Furthermore, the uniqueness of our work leverages the operational forecasts issued twice a day at the weather service of the Tuscany regional (LaMMA). In fact, with such data, we computed a custom and long-term (about ten months) background error covariance matrix. The quality of the background error covariance matrix is crucial in a 3D-Var assimilation step [40], because its statistics determine how the observations spread in the model space and how each model grid point contains dynamically balanced increments. We claim that this represents a consistent improvement with respect to previous studies, which computed the background error covariance matrix for a shorter time period (less than three months [39,40,45,49,50]) or used the default matrix provided by the software developers [42]. Finally, to evaluate the results obtained when assimilating observations from weather stations and radars and compare such results with those obtained without assimilation, we developed an object-oriented verification method, aimed at evaluating both the intensity and position errors of rainfall predictions.

The paper is organized as follows: in Section 2 we describe the meteorological context of the Livorno case and the modeling setup implemented; in Section 3 we give the details about the method used to evaluate model outputs; in Section 4 we present the main results; in Section 5 we discuss the implications of the results and the limitations of this study; and we describe our vision of the path forward in Section 6.

2. Materials and Methods

2.1. Synoptic Conditions

During the first hours of the 9 September, a large trough deepening over the Eastern Atlantic Ocean approached the Mediterranean Sea. At 0000 UTC of 10 September 2017, the axis of the trough was oriented from the Scandinavian Peninsula to the Western Mediterranean Sea (see Figure 1a). The trough slowly moved eastward, causing the deepening of a low-pressure area (see Figure 1b) over the Ligurian sea and lee of the Italian Alps [51,52]. As a result, the pre-existing warm and humid air masses over the Tyrrhenian Sea interacted with the dry and relatively cold flow from France (see Figure 1b and Figure 7b in [45]), making the environment conducive for persistent precipitation systems. Furthermore, upper level divergence, sustained convective available potential energy values in the range 500–1000 J· kg${}^{-1}$, and strong wind shear (plots not shown) favored convective motions. A well-defined line of convergence of low-level winds over the Tyrrhenian Sea acted as a trigger to overcome the convective inhibition energy (see blue wind vectors in Figure 1b).

Following the classification of [54], the Livorno case belongs to the HPEs of short duration (less than 12 h) with spatial extents smaller than $50\times 50$ km${}^2$.

2.2. Observed Precipitation Measurements

The Livorno case can be described by using the data obtained from the weather stations managed by the Hydrological Service of Tuscany (www.sir.toscana.it, (accessed on 29 April 2021)). The network is composed by approximately 400 rain gauges, located over an area of about 23,000 km${}^2$, reporting conditions every 15 min [55]. In Figure 2a we report the precipitation map of the accumulated precipitation in the 6 h period ending on 0300 UTC of 10 September 2017. One precipitation maxima, located in the coastal area, and exceeding a cumulative rainfall amount of about 240 mm can be noticed, to the south of the Livorno township (indicated with the red triangle in Figure 2a). The extreme is located in a widespread area of total precipitation exceeding 100 mm, which covers the entire central coast of Tuscany.

The variability of both the duration and intensity of precipitation across the territory played a major role in the downstream hydrological responses. The sector south of the Livorno township, which was responsible for the flash floods in the urban areas, produced precipitation that was limited in time: 210.2 mm was registered in the 2 h period between 0045 UTC and 0245 UTC, 10 September (see Figure 2b reporting the data collected at the station of Valle Benedetta). This corresponds to 86% of the total event precipitation registered by this station. Another weather station, Quercianella, still located in the southern sector, measured similar values: 188.6 mm in 2 h, corresponding to 89% of the total precipitation during the event. Analogies in the maximum hourly rain rate are even stronger, since the first station measured 120.8 mm in 1 h, and the second 1 mm more. For both rain gauges, [57] estimated a return period of$>200$ years for the 1 and 3 h duration rainfall events. The sector north of the town of Livorno, represented in Figure 2b by the stations of Bocca d’Arno and Metato, was characterized by a peak rain rate between 50 and 70 mm in 1 h and a total event duration of 3 to 4 h. This precipitation was mostly concentrated between 2000/2100 UTC of 9 September and 0000 UTC of 10 September. In the rest of the region (plots not shown), precipitation was either intense and brief (further south: Castellina Marittima rain gauge, return period $>200$ years for 1 and 3 h duration rainfall events [57]), or long lasting and less intense (north of Livorno: Coltano and Stagno rain gauges, return period $\simeq 90$ years for 12 h duration rainfall events [57]).

2.3. Radar Data

Several weather radars were operational during the rainfall event, belonging to different networks: the Italian national [58,59], the French national [60], and the Tuscany regional [48,61] ones. Among these available radars, two were selected for this work: the Aléria and the Livorno radars. The former is a Doppler S-band (the signal frequency is 2802 MHz) system located in Aléria in the central Eastern Corsica and was detecting the severe weather system from a distance of approximately 70–80 km. The latter is a non-Doppler, single polarization X-band radar (its frequency is 9410 MHz) located in the Livorno harbor, and it monitored the dynamics of the meteorological event moving towards the radar site. These systems have been designed and produced by different manufacturers and have different technical characteristics, particularly in regard to the operational frequencies, system dimensions, and architectures [48]. The Aléria radar (S-band) is characterized by low attenuation impacts on radio signals due to the atmosphere and rainfall, and consequently can reliably monitor rainfall fields up to a distance of approximately 100–150 km. The Livorno radar system (X-band) is more sensitive to smaller rainfall drops and can provide a finer resolution, but with a limited range due to the strong effects of atmosphere and rainfall attenuation on radio signals. Moreover, the effects of path and wet radome attenuation due to the rain falling directly over the radar system were increased in this case by the strong precipitation occurring at the Livorno site. The calibration of the two radars is performed during maintenance procedures, using reference radiofrequency sources directly sent to the receiver, so that the receiver dynamics can be monitored and correctly taken into account in the radar constant. In this paper, reflectivity data were used for the assimilation into the WRF model (the Aléria radar provided also radial velocity data).

These two radar systems detected the rainfall event, since its origin was over the sea between North Corsica and the Ligurian Gulf, while approaching the Tuscany coasts. The Livorno radar localized the weather event and followed its dynamics with high spatial detail, as shown in Figure 3a. Then the longitudinal shape of the strong precipitation front, well detected when approaching the Tuscany coasts (Figure 3c), underestimated the reflectivity intensity of the storm, at least during its beginning and early development on the sea (Figure 3b,d). Later, when the storm reached Livorno, the Aléria radar detected the strong rainfall and its persistence over the area (Figure 3f), whereas the Livorno radar system underestimated the intensity of precipitation, presumably due to a significant path and radome attenuation (Figure 3e). These two radars, in association with the Italian weather radar network [58,59], were used by [61] to characterize the Livorno event and were found to provide valuable information for the quantitative precipitation estimation.

2.4. Satellite and Lighting Data

The deep convection of the Livorno event was detected by the Rapid Scan High Rate (5-min) mode of the Spinning Enhanced Visible and Infrared Imager on board the Meteosat Second Generation satellite [62]. Satellite images (maps not shown) indicate a cloud top brightness temperature below −65 ${}^{\circ }$C, a temperature that can be found at approximately 12.5 km height, according to the atmospheric sounding measured at 0000 UTC on 10 September at Pratica di Mare (data not shown), a sounding station located approximately 250 km south of the area of interest. Further and in-depth analyses of the Livorno case using satellite-based observations are available in [63].

The deep convection of the Livorno event produced (see Figure 4a) an average of four lightning strikes per second over northern Tuscany in the six hours between 2150 UTC of 9 September and 0350 UTC of 10 September; this number is approximately 10% of the global flash rate [64].

For those rain gauges that recorded the greatest amounts of precipitation (Valle Benedetta and Quercianella), corresponding lightning time series graphs are shown in Figure 4b,c as histograms of the number of strokes per hour within ±0.1 degree of the two weather stations centers, from 2200 UTC of 9 September to 0400 UTC of 10 September 2017. Totals of 9015 and 7177 strikes were detected during the 6 h period for Quercianella and Valle Benedetta stations, respectively. The highest number of flashes was recorded at 0100 UTC on 10 September for both Quercianella and Valle Benedetta (3053 and 2528, respectively).

2.5. Modeling Setup

The numerical model used in this work to simulate the 9–10 September Livorno case was the WRF model [65]. The Mesoscale and Microscale Meteorology (MMM) Laboratory at the National Center for Atmospheric Research (NCAR) has led the development of the WRF model since its inception in the late 1990s. It is a fully compressible, Eulerian, non-hydrostatic mesoscale model, designed to provide accurate numerical weather forecasts both for research activities and for practical operations. In this work, we implemented the Advanced Research WRF (ARW) version of the model updated to version 4.0 (June 2018). The model’s dynamics, equations, and numerical schemes implemented in the WRF-ARW core are fully described in [65,66]; the model’s physics, including the different options available, is described in [67]. In

Table 1. Basic settings of the WRF model simulations.

Parameter	Value
Rows × Columns	440 × 400
Grid spacing	3 km
Vertical levels	50
Time step	12 s
Cumulus convection	explicit (no parameterisation)
Micro-physics option	Thompson scheme [68]
Boundary-layer option	Yonsei University [69]
Land-surface option	Unified Noah model [70]
Radiation option	Rapid radiative transfer model [71]
Turbulence option	Yonsei University+2D Smagorisnski [69]

we give a short summary of the basic model settings used in this study: the domain of integration is indicated in Figure 1 with the red rectangle; the grid spacing is 3 km; a stretched vertical grid with 50 vertical levels up to a height of 100 hPa was applied, with 12 levels in the lowest 1500 m. The time step was 12 s, and we stress the fact that the model was set with the explicit treatment of convective processes. In Table 1, we also give references to the micro-physics, boundary-layer, land-surface, radiation, and turbulence schemes adopted. All the WRF simulations (with and without assimilation) performed in this study share the same basic set of physical parameterizations listed in Table 1. This configuration mirrors the one that was operational at the weather service of the Tuscany regional (see www.lamma.toscana.it, (accessed on 29 April 2021)) at the time of the Livorno event, although in this study we used a more updated WRF model version (4.0 instead of 3.9).

To start the WRF simulations, initial and boundary conditions were derived from the ECMWF-IFS global model data operational at the time of the event (model cycle 43r3). The spectral resolution was TCO1279, (where CO means that a cubic-octahedral grid was used), which corresponds to a grid spacing of approximately 9 km; the number of vertical levels was 137 and the top of the atmosphere was set to 0.01 hPa.

The control forecast was initialized at 1200 UTC on 9 September 2017, and the forecast length was set to 15 h (ending time was 0300 UTC, 10 September). To test the impact of assimilating conventional and radar data, we utilized the WRF data assimilation (WRFDA) system [72]. The WRFDA software was used to perform 3 h cycling 3D-Var data assimilation using the rapid update cycle approach, similarly to what is found in [73]. We implemented an assimilation step every 3 h starting from 1200 UTC on 9 September up to 1800 UTC of the same day with an assimilation window of ±30 min; i.e., all the observations between $t-30$ minutes and $t+30$ minutes were assumed to be valid at the analysis time t. The 3D-Var implementation of the WRFDA package relies on previous developments designed for the fifth-generation Penn State/NCAR Mesoscale Model [74]. Although detailed descriptions of the 3D-Var method can be found in [41,75,76], the technique can be summarized as follows: the basic idea is to estimate the optimal state of the atmosphere (in the model space) by using the observations available, a short-term forecast (often referred as first guess or background), and information about error statistics on both the observations and the model. Let $t\in \mathbb{R}$ be the time of the analysis and let $\mathbf{x}=\mathbf{x}\left(t\right)\in {\mathbb{R}}^n$ be the model analysis at time t. Using an iterative process, the 3D-Var method looks for the minimum value of the cost function $J\left(\mathbf{x}\right)$, defined as:

$$J\left(\mathbf{x}\right)=\frac{1}{2}\left\{{(\mathbf{x}-{\mathbf{x}}_{\mathbf{b}})}^T{\mathbf{B}}^{-1}(\mathbf{x}-{\mathbf{x}}_{\mathbf{b}})+{({\mathbf{y}}_{\mathbf{o}}-H\left(\mathbf{x}\right))}^T{\mathbf{R}}^{-1}({\mathbf{y}}_{\mathbf{o}}-H\left(\mathbf{x}\right))\right\},$$

(1)

where, following the notations of [77], ${\mathbf{x}}_{\mathbf{b}}\in {\mathbb{R}}^n$ is the background model state at time t, $\mathbf{B}\in {\mathbb{R}}^{n\times n}$ is the covariance matrix of the background errors, ${\mathbf{y}}_{\mathbf{o}}\in {\mathbb{R}}^p$ is the vector of observations ($p<n$) at time t, $H:{\mathbb{R}}^n\mapsto {\mathbb{R}}^p$ is the observation operator that transforms the variables from the model state to the observation space, and $\mathbf{R}\in {\mathbb{R}}^{p\times p}$ is the covariance matrix of observation errors.

Accurate estimates of the $\mathbf{B}$ and $\mathbf{R}$ matrices determine the quality of the analysis. As regards the $\mathbf{R}$ matrix, default values provided by the WRFDA software (see User’s Guide [65]) were used for diagonal elements, whereas off-diagonal’s elements were set to zeros. In fact, correlations between different instruments are usually assumed as null. Although some studies demonstrated that including inter-correlations in the $\mathbf{R}$ matrix may provide better analyses [78], this mainly holds when satellite radiance data are considered [79,80]; thus, we claim that the assumption of the $\mathbf{R}$’s structure was valid in our study. We computed the background error correlation matrix $\mathbf{B}$ by means of the National Meteorological Center (NMC) method [81], which estimates the value of the elements of the matrix statistically, by averaging the differences between two short-term forecasts valid at the same time but initialized one shortly after the other (we used 12 h later). Although the NMC method is described as a very crude first step by its authors, it is still widely used nowadays. The $\mathbf{B}$ matrix we used is the result of the NMC method applied for ten months (from December 2018 to October 2019).

Following [82], the observation operator H for the radial velocity $V_r$ is defined by the relationship:

$$V_r=u\frac{x-x_i}{r_i}+v\frac{y-y_i}{r_i}+(w-v_T)\frac{z-z_i}{r_i},$$

(2)

where $(u,v,w)$ are the modeled wind components, $(x_i,y_i,z_i)$ are the coordinates of the radar antenna, $(x,y,z)$ are the coordinates of the radar observation, $r_i$ is the distance between $(x,y,z)$ and $(x_i,y_i,z_i)$, and $v_T$ is the mass-weighted terminal velocity of the precipitation. To calculate $v_T$, the default formula [83] was used—that is:

$$v_T=5.4·{(p_0/\overline{p})}^{0.4}·q_r^{0.125},$$

(3)

where $p_0$ is the surface pressure, $\overline{p}$ is the base-state pressure, and $q_r$ (unit g· kg${}^{-1}$) is the model predicted rainwater mixing ratio. Let Z be the reflectivity data expressed in dBZ; then the nonlinear $Z-q_r$ relationship is defined by [82]:

$$Z=c_1+c_2·{log}_{10}\left(\rho q_r\right),$$

(4)

where $c_1=43.1$ and $c_2=17.5$ are constants and $\rho$ (unit kg· m${}^{-3}$) is the air density. Following the performances achieved by [43], radar reflectivity data were assimilated using the indirect technique proposed by [76], that is, by inverting the $Z-q_r$ relation given in Equation (4) and assimilating the observed rainwater mixing ratio estimated from reflectivity values.

Before assimilation, radar reflectivity data undergo to a thinning procedure. This consists of eliminating ground/sea clutter and georeferencing radar data onto the model grid; then for each model grid point the procedure selects the maximum reflectivity value (and the corresponding radial velocity) among radar observations belonging to that grid point. Finally, to reduce the spin-up in the first hours of integration, we implemented the digital diabatic filter initialization [84] at each starting time of the WRF simulationsw.

To assess the impact of assimilating different types of observations, three experimental forecasts were deployed, and their results are compared with those obtained with the control forecast (named with the code C), which did not assimilate any observational dataset. The first experimental forecast (code S) assimilated data from ground weather stations; the second (code $S+R$) assimilated radar data in addition to ground weather station data. Conventional data assimilated in both experiments were: pressure, temperature, and relative humidity at the 2 m above ground level; wind speed and direction at the 10 m above ground level. The total number of weather stations reporting valid data was 1036 for the assimilation step performed at 1500 UTC on 9 September and 1051 for the 1800 UTC assimilation. Remote sensed data assimilated were reflectivity data from the X- and S-band radars and radial velocity data from the Aléria S-band radar. The locations of weather stations and radars are indicated with the blue points and red triangles, respectively, in Figure 5.

To investigate the role of sea surface temperatures (SSTs) in the Livorno case, we performed a further numerical experiment (code $S+R+M$), which shared the same settings and design of the $S+R$ experiment, but included a simplified ocean model that modified its status (SSTs and ocean mixed layer depth) thanks to the interface heat, momentum, and mass fluxes. In fact, the interaction between the atmosphere and the ocean can play a key role in the formation and intensification of extreme atmospheric events [85] through the energy fluxes, which can impact the generation of HPEs by modifying the structure of the boundary layer, the distribution of wind fields, and therefore the position of the convergence line [86,87]. Historically, numerical weather models implemented the SSTs as a static boundary, taken from the global models or from low-resolution satellite datasets. However, when following this approach, it is difficult to reproduce complex energy feedback between air and sea [88], particularly in coastal areas. In order to limit this impact, without slowing down the calculation time with the coupling of an ocean model, we proceeded by using a “slab ocean,” also known as a simple ocean mixed layer [89], which updates the SSTs every hour, according to the energy fluxes at the air–sea interface. This approach is based on a simplified model that aims at implementing the SSTs and ocean mixed layer depth into the WRF model. This one-dimensional approach was initialized with SST data, a mixed layer depth value of 40 m, and a lapse rate 0.14 K · m${}^{-1}$. These data were provided by the Copernicus Marine Services (CMEMS) numerical model. WRF modifies the SST values according to the energy and momentum fluxes at the interface, not taking into account the ocean dynamics, but making each grid point evolve as a function of the energy budget at the surface.

Table 2. Codes assigned to name the forecasts and descriptions of the data used in each assimilated forecast.

Forecast Code	Data Assimilated
C	none (control run)
S	Convectional data from weather stations:
	pressure, 2-m temperature
	2-m relative humidity
	10-m wind speed and direction
$S+R$	as in S plus reflectivity data from X- and S-band radars,
	and radial velocity data from the S-band radar
$S+R+M$	as in $S+R$
	(it differs from the $S+R$ experiment because
	it implements a simplified marine model)

shows a summary of the various forecasts and the codes adopted to name them. Figure 6 shows a scheme of the modeling setup.

3. Quantitative Precipitation Forecast Verification

When evaluating the performances of different forecasts, a quantitative evaluation of the spatial agreement between predicted and observed values is crucial. A direct numerical comparison can be misleading, especially for a variable whose spatial distribution is highly concentrated in a small area, as happens during HPEs. In fact, any forecast that correctly predicts the occurrence of highly localized heavy rain, may incur the so called “double penalty” error [90] if it places the event in a nearby area, producing, for example, a root mean squared error ($RMSE$, see [91]) higher than another forecast which completely misses the prediction. To overcome such a limitation, object-based verification methods have been developed by the scientific community [92], and currently, software packages for practical applications exist [93]. These methods exhibit some drawbacks, such as the smoothing and filtering the observations undergo, and the large number of parameters whose settings are somewhat arbitrary.

Considering the features of the phenomenon under examination, we have devised a simple and robust ad hoc verification method, which does not alter in any way the recorded rainfall values. No arbitrary or subjective parameters are required, except for the ones needed to define the area where the event occurred, namely, the center and radius of the circle by which the area of interest is inscribed. For the flash flood under examination, we chose a circle with a radius of 0.4 degrees in the longitude–latitude forecast domain and centered at the point $P_C$ with coordinates $(x_{C,lon},x_{C,lat})=(10.{64}^{\circ }\mathrm{E},43.{75}^{\circ }\mathrm{N})$, about 30 km northeast of Livorno. This area includes 91 rain gauges located at points $P_i$ ($i=1,\dots ,91$), which measured the highest rainfall rate in the 6 h period ending on 0300 UTC, 10 September (see Figure 2). By applying affine transformations (more precisely, rigid roto-translations preserving the Euclidean distance; these transformations are usually denoted as the $SE\left(2\right)$ group) to the 91 position vectors ${\mathbf{x}}_i=(P_i-P_C)\in {\mathbb{R}}^2$, we searched for the transformation while minimizing the $RMSE$ between the cumulative rain $r_i$ measured at station $P_i$ and the $f\left({\mathbf{y}}_i\right)$ predicted by the forecast model f at the transformed points ${\mathbf{y}}_i\in {\mathbb{R}}^2$ (see Figure 7). Since we applied a roto-translation, any point ${\mathbf{y}}_i$ is determined by:

$${\mathbf{y}}_i=R_{\varphi }\left({\mathbf{x}}_i\right)+{\mathbf{x}}_0,$$

(5)

where $\varphi$ is the rotation angle, $R_{\varphi }\in {\mathbb{R}}^{2\times 2}$ is the rotation operator, and has the form:

$$R_{\varphi }=\left[\begin{array}{cc}cos\varphi & -sin\varphi \\ sin\varphi & cos\varphi \end{array}\right],$$

(6)

and ${\mathbf{x}}_0=(x_{0,lon},x_{0,lat})\in {\mathbb{R}}^2$ is the translation vector. All the parameters of the roto-translation are defined in the longitude–latitude space. We note that the distance, in km, between $P_C$ and the center of any circle translated by a vector ${\mathbf{x}}_0$ is given by:

$$d_0\simeq R\sqrt{{cos}^2\left(x_{C,lat}\right)·x_{0,lon}^2+x_{0,lat}^2},$$

(7)

where $x_{0,lon}$, $x_{0,lat}$, and $x_{C,lat}$ are expressed in radians, and R is the Earth’s radius (approximately 6370 km). The minimization of the positive defined error function

$$RMSE=\sqrt{\langle {(f\left({\mathbf{y}}_i\right)-r_i)}^2\rangle },$$

(8)

where $\langle ·\rangle$ stands for the average operator over the index i, determines the particular set of parameters for which the agreement between predicted and observed cumulative rain is the highest possible, at least in terms of $RMSE$. In this way, the location and intensity errors are disentangled, the former of which is given by the affine transformation parameter’s amplitude ($\varphi$ and ${\mathbf{x}}_0$), and the latter by the minimized $RMSE$.

To contend with the validation method described above, standard verification scores were computed. Observed rainfall data registered at the 91 rain gauges were compared with corresponding modeled data, which were extracted at the 91 nearest grid points containing rain gauge locations. Standard verification statistics considered were: $RMSE$ (defined in Equation (8)), mean error ($ME$), and the multiplicative bias ($mbias$). Let $r_i$ and $f\left(P_i\right)$ be the observed and modeled rainfall values at location $P_i$, respectively; then $ME$ and $mbias$ can be defined as follows (see also [91]):

$$\begin{array}{cc}\hfill ME=& \langle f\left(P_i\right)-r_i\rangle ,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill mbias=& \frac{\langle f\left(P_i\right)\rangle }{\langle r_i\rangle }.\hfill \end{array}$$

(10)

Furthermore, to summarize multiple aspects of the model’s performance in a single diagram, we made use of the performance diagrams (see [94]). Such diagrams plot four skill measures of dichotomous (yes/no) forecasts, probability of detection ($POD$), success ratio ($SR$), frequency bias ($bias$), and critical success index ($CSI$), using the $2\times 2$ contingency table shown in

Table 3. The $2\times 2$ contingency table.

		Event Observed
		yes	no
Event Forecast	yes	A	B
	no	C	D

.

The four skill measures are defined as follows:

$$\begin{array}{cc}\hfill POD=& \frac{A}{A+C},\\ \hfill SR=& 1-\frac{B}{A+B},\\ \hfill bias=& \frac{A+B}{A+C},\\ \hfill CSI=& \frac{A}{A+B+C}.\end{array}$$

4. Results

In Figure 8, we show with the shaded colours the quantitative precipitation forecast (QPF) in the 6 h period ending at 0300 UTC, 20 September for the four experiments listed in Table 2.

In each panel, the gray circle indicates the area containing the 91 rain gauges (closed gray points) that registered the heaviest precipitation. The visual inspection of Figure 8 indicates that the C and S forecasts provided QPF values close to zero in the Livorno area. On the other hand, both numerical experiments predicted QPF peaks up to 100–120 mm in an area a few tens of kilometers (approximately 40–60 km) to the north of the area of interest, where rain gauges registered accumulated precipitation of up to 30–40 mm (see Figure 2a). Both the $S+R$ and $S+R+M$ forecasts predicted high QPF values inland of the Livorno area; however, QPF maxima (in the ranges 125–150 and 80–100 mm for the $S+R$ and $S+R+M$ forecasts, respectively) largely underestimated actual observed values, which reached approximately 200–230 mm (see Figure 2b).

In

Table 4. Verification scores, as defined in Equations (8)–(10), for the four experimental runs.

Forecast Code	$\mathit{RMSE}$	$\mathit{ME}$	$\mathit{mbias}$
C	65.8	−50.3	<${10}^{-3}$
S	65.9	−50.5	<${10}^{-3}$
$S+R$	46.3	−2.5	0.95
$S+R+M$	40.4	−2.9	0.94

, we show the verification skill scores defined in Equations (8)–(10). The scores were computed using observed rainfall registered at the 91 rain gauges shown in the gray circle of Figure 8 and the corresponding modeled precipitation extracted at the grid points closest to the rain gauge locations.

The $S+R$ and $S+R+M$ forecasts exhibited the more satisfactory scores; in fact, the $RMSE$s were approximately 40/45 mm, whereas C and S errors were higher and reached approximately 65 mm. The $ME$s and biases demonstrate the underestimation of all the four predictions; underestimations were remarkable for C and S ($ME$ ≃ −50 mm and $mbias<{10}^{-3}$). On the other hand, the $S+R$ and $S+R+M$ forecasts provided more skillful scores, with the $mbias$ being close to the perfect score of 1 and $ME$s being approximately −2.5/−3 mm.

In Figure 9, we show the performance diagram for selected rainfall thresholds corresponding to approximately the 20th, 40th, 60th, and 80th percentiles of the observational dataset. The plots show that, for precipitation thresholds equal to 17, 34, and 49 mm, the $S+R$ and $S+R+M$ forecasts behaved similarly and outperformed the C and S forecasts, which had no skill, since both matching points lie close to the bottom-right corner. As the precipitation threshold equaled 83 mm, which approximately corresponds to the 80th percentile of the observations, the $S+R$ forecast did not provide any valuable information, since the matching point lies close to the bottom-left corner. The $S+R+M$ forecast has some skill as regards the $bias$ (approximately 0.5) and the $SR$ score (approximately 0.3), whereas the performances of $CSI$ (≃10) and $POD$ (<0.2) were limited. Skill scores for precipitation thresholds greater than the 80th percentile of the observations were poor or null for all the forecasts (performance diagrams not shown).

In each panel of Figure 8, we also show the result of the roto-translation of the gray circle minimizing the $RMSE$ between the observed and predicted rainfall data (the transformed circle and rain gauge locations are indicated with the black color). The angle of rotation $\varphi$ was measured, by definition, counterclockwise from north. For each experimental run, the parameters of the $RMSE$ minimizing roto-translation are shown in

Table 5. Parameters of the roto-translation minimizing the $RMSE$ between observed and predicted values extracted at the 91 rain gauge locations located within the gray and black circles, respectively, shown in Figure 8. $RMS{E}_{min}$ is the error after the minimizing roto-translation; $x_{0,lon}$ and $x_{0,lat}$ are the coordinates of the translation vector ${\mathbf{x}}_0$ defined in Equation (5); $d_0$ indicates the distance, in km, between the centers of the black and gray circles (i.e, the dashed segment of Figure 8) and is given in Equation (7). The angle $\varphi$ indicates the rotation amplitude with respect to north, counterclockwise.

Forecast Code	${\mathit{RMSE}}_{\mathit{min}}$	${\mathit{x}}_{0,\mathit{lon}}$	${\mathit{x}}_{0,\mathit{lat}}$	${\mathit{d}}_0$	$\mathit{\varphi }$
C	29.7	−0.36	+0.87	100	${343}^{\circ }$
S	29.0	−0.30	+0.75	86	${333}^{\circ }$
$S+R$	29.9	−0.24	−0.18	27	${196}^{\circ }$
$S+R+M$	30.1	+0.48	+0.41	60	${351}^{\circ }$

.

After the roto-translations, all four experiments exhibited similar $RMS{E}_{min}$; experiment S provided the lowest value (29.0 mm) and experiment $S+R+M$ had the largest one (30.1 mm). However, the parameter $d_0$ indicates that the $S+R$ experiments displaced the transformed black circle by approximately 27 km, whereas the other experiments provided higher displacements, ranging from 60 to 100 km. We stress the fact that not only was the displacement lower when radar data were assimilated, but this was also the case when radar were oriented along the axis of the perturbation hitting the Livorno township (from southwest to northeast). The angle of rotation $\varphi$ reported in Table 5 spans from ${196}^{\circ }$ for $S+R$ to ${351}^{\circ }$ for $S+R+M$.

To further evaluate the impact of assimilating conventional and radar data, in Figure 10 we show the 10 m wind speed and direction at 0000 UTC on 10 September for the four numerical experiments. All four forecasts reconstructed a well defined convergence line over the Ligurian Sea between southerly and westerly winds, which is responsible for triggering convective rainfall. However, in C experiment, that convergence line was positioned approximately 70 km to the north of the Livorno area, causing the wrong displacement of precipitation. The impact of the conventional data (experiment S) was negligible and did not modify the position of the convergence line with respect to the control run C (see Figure 10a,c). On the other hand, assimilating radar data substantially modified the 10 m wind speed and direction forecast, producing a better localization of the convergence line, positioned close to the area that registered the heaviest precipitation (see panel (c) in Figure 10). The implementation of the simplified ocean model (see panel (d) in Figure 10) generated prefrontal winds (southern part of the convergence line) more intense by about 8 m/s, and a local pressure decrease of about 1.5 hPa (map not shown); as a consequence, the front-genesis evolved more quickly with greater intensity in the event area.

To understand how the precipitation field is modified in the experimental runs, in Figure 11, we show the specific humidity averaged over the entire column of atmosphere at 0000 UTC on 10 September. We stress the fact that these maps are the results of both the assimilation procedure and the model evolution in the first hours of integration. The $S+R$ run (panel (c) in Figure 11) clearly shows a tongue of high water vapor intensity extending from the northern tip of Corsica Island to the northern Tuscany coasts (and further inland), which was not simulated, nor was it in the control run C (panel (a) in Figure 11) or in the assimilated run S (panel (b) in Figure 11). The addition of the simplified marine model to the $S+R$ experiment (panel (d) in Figure 11) moved this line of high water vapor content northwards, close to the Ligurian coasts, slightly increasing its intensity, with values greater than 6.3 g· kg${}^{-1}$.

5. Discussion

As regards the infamous rainfall event that occurred in the Livorno area (Central Italy) during the night between 9 and 10 September 2017, we evaluated possible improvements arising from the assimilation of both conventional and radar data in the setup of the WRF model, which was operational at the meteorological agency of the Tuscany region (LaMMA) at the time of the event. In addition, we implemented a simplified ocean model, having a relatively low computational cost, to evaluate any further improvements due to a better description of energy exchanges between air and sea, which are known to play a crucial role in the formation and intensification of HPEs [85,86,87]. Such experimental forecasts were compared against a control run, which did not assimilate any data, nor implement the simplified ocean model.

Since no reliable gridded observed rainfall dataset was available for the Livorno case, as ground-truth data, we used the observations collected at the rain gauges managed by the Hydrological Service of Tuscany [55]. To assess the accuracy of QPFs, we introduced a novel method (see the sketch in Figure 7), aimed at a joint evaluation of both the intensity and the position errors of the forecasts, while preserving the recorded rainfall data (i.e., no interpolation or filtering of the observations was performed to comply with the model data). To further validate the predictions, we also computed standard verification scores and constructed performance diagrams for predefined precipitation thresholds.

Figure 8a demonstrates that the control simulation C was not able to correctly predict rainfall peaks in the Livorno area; this is not new information and confirms the findings of similar studies [43,45,63]. However, the results shown in Figure 8 and Table 5 suggest that the $S+R$ forecast is more valuable. In fact, looking at the position error, the $d_0$ parameter is the lowest, indicating the optimal positioning of the transformed black circle close to the original gray circle. We note that the parameter $d_0$ for the $S+R$ forecast is about one third the length of $d_0$ for the S forecast and about one quarter the length of $d_0$ for the C control run, indicating the relatively higher accuracy of $S+R$ with respect to C and S. For all the forecasts, with the exception of the $S+R$ case, the $RMSE$ minimization transformation moved the circle including the representative rain observation points inland, next to the area with the highest predicted rain and applied translation with slight rotation (the parameter $\varphi$ is close to ${360}^{\circ }$). In the $R+S$ case, on the contrary, the translation moved the circle towards the seaside with a rotation of $\varphi ={196}^{\circ }$ that fliped coastline and inland observation points, reversing north with south. Indeed, the $S+R$ forecast located the heaviest rain inland, instead of near the coast as actually happened, and further south. Consequently, the minimization process tends to rotate the circle of stations by about 180${}^{\circ }$, flipping coastline and inland rain gauges, and northern with southern ones. The translation towards the seaside indicates that the heaviest rain predicted inland was even heavier than the maximum values recorded on the coast. As regards the intensity error, the lowest $RMS{E}_{min}$ value was obtained for the S forecast; however, all the forecasts provided comparable intensity errors (at least for the selected rain gauges).

Similar conclusions can be drawn looking at the results shown in Table 4 and Figure 9. Both the $S+R$ and $S+R+M$ forecasts outperformed C and S, yielding better scores. However, in the absence of any assessment regarding the position errors, it is difficult to evaluate which forecast between $S+R$ and $S+R+M$ provided the more valuable predictions. We deduce that the position/intensity error information shown in Table 5 and Figure 8 confirms and complements the information obtained from standard verification skill scores shown.

Our conclusions agree with those of [45], who found that assimilating radar (and lighting) data provides better QPFs by changing the water vapor mixing ratio during the assimilation step. Incidentally, we note that the performance diagram shown in Figure 9b agrees well with Figure 16f in [45] and provides similar, or slightly better, results, supporting common conclusions. However, we acknowledge that a direct comparison should be treated with care, because the radar data and assimilation techniques were different, and the period over which rainfall data were accumulated was shorter (3 vs. 6 h).

Our findings further confirm [39,41,42,43,45] that 3D-Var radar data assimilation is an effective method for improving QPFs by correcting the initial conditions of limited-area weather numerical models. In fact, concerning the assimilation of radial velocity data (available at the Aléria S-band radar only), Figure 10 suggests that remote sensed data are beneficial for better positioning of the low-level wind convergence line. Moreover, the indirect assimilation of radar reflectivity (available at both the Aléria and Livorno radars) and in turn the modification of water vapor mixing ratio (see Figure 11) are supposed to provide an environment that is conducive for convection; convection carries large contents of moisture higher into the atmosphere (see panels (c) and (d) in Figure 11) and energetically favors cloud formation. Thus, water vapor’s spatial and temporal variations indirectly provide some indications on how the different model setups may impact the event prediction. Our findings agree with those of [41], who discussed how radial velocity assimilation has a large impact on wind velocity, whereas reflectivity data have a direct impact on hydrometeor analyses. We conclude that, since the origin of the Livorno event took place on the sea, the assimilation of conventional data, available only on land, had a low impact on the accuracy of QPF for the S experiment. On the other hand, assimilating radar data, which were available for the offshore stage, provided pertinent and crucial information regarding rainfall triggering and atmospheric moisture content. Furthermore, we claim that such remote sensed data modified the model dynamics in a way that persisted a few hours (up to 9 h) after the assimilation step (which ended at 1800 UTC on 9 September).

Given the outputs of the $S+R+M$ experiment (see panel (d) in Figure 8, Figure 10 and Figure 11), the forecasts were not dramatically affected by the use of the simplified ocean model. This is consistent with the conclusions drawn by [95], who found that the evolution of the SSTs during the model integration has a marginal impact on the short-range predictions (forecast length less than 18 h). We can deduce that, for the Livorno case, the exchange of both heat and water vapor between the air and sea did not play a crucial role, as also found by [47]. In fact, the authors argued that feeding SSTs estimates from Sentinel-3 satellite observations into the WRF model did not improve the Livorno QPF for a high precipitation rate. This happened possibly because at the kilometer scale, the SSTs determine the intensity of the warm low-level jet [88], which, in our case, was partially corrected by the ingestion of radial velocity radar data. In our $S+R+M$ forecast, it was observed that the heat southerly fluxes were more intense by about 10–15% if compared to $S+R$ data (map not shown). We found that, although the skill scores of $S+R$ and $S+R+M$ were similar (see Table 4 and Figure 9), the simplified ocean model was able to modify the precipitation field by reducing the QPF peaks close to the coast and thus producing a position error (see panels (c) and (d) in Figure 8 and Table 5).

We claim that assimilating radar data turned out to be effective because we used a covariance matrix of the background errors $\mathbf{B}$ that was the result of a long-term application of the NMC method (approximately ten months). In fact, to estimate $\mathbf{B}$, we used the operational forecasts issued twice a day at the regional meteorological service of Tuscany (LaMMA), which shared the same setup (in terms of resolution, number of vertical levels, physical parameterisations, etc.) as the runs presented herein. This is a remarkable improvement with respect to previous similar studies; to mention a few works, Refs. [39,43] employed a 1-week and a 1-month period, respectively, to compute the $\mathbf{B}$ matrix, whereas [45,50] applied the NMC method during the Hydrological cycle in the Mediterranean Experiment—First Special Observing Period (HyMeX-SOP1, which lasted for about two months in 2012, see [20]). It is also worth mentioning [42] who used the default matrix provided by the WRFDA system, which is produced with global data and its use for regional cases is sometimes discouraged [65]. Additionally, more recent works (see [96]) used a time period shorter than one month to compute the covariance matrices of the background errors. However, we acknowledge that further investigations are needed to assess the impact of a long-term application of the NMC method to build a high-level, quality $\mathbf{B}$ matrix and obtain better QPF data. We stress the fact that our study focused on a single event; to draw more general and statistically sound conclusions, the promising methodology and results have to be evaluated for different locations and similar events, and that is what we intend to do. However such experiments must rely on the availability of the same X- and S-band radar data used in this study; this limits the period to HPEs occurred since 2014 [48].

6. Conclusions

Although tremendous improvements were achieved in recent years, the operational forecasting of HPEs in the Western Mediterranean basin still remains a challenging task, particularly in the context of a warming climate. As concerns the Livorno case, the main findings of this work are:

• Assimilating reflectivity data from X- and S-band radars and radial velocity data from S-band radar significantly improves the descriptions of atmospheric humidity content and low-level winds, resulting in better QPFs.
• The application of a simplified ocean model, although modifying the low-level jet associated with the event, scarcely impacts the short-range forecast (length shorter than 12 h) of precipitation.
• The novel QPF verification method introduced in this paper, based on roto-translation RMSE-minimisation, confirmed and reinforced the results achieved with standard verification scores, thereby providing more information about the position error of the WRF simulations.

The conclusions we drew support the deployment of a dense network of relatively small radars in the coastal areas of the Western Mediterranean Sea. Such high-resolution (both spatial and temporal) data may have beneficial impacts on short-range numerical weather predictions. In fact, by extracting the maximum value from local observations, such as those collected by X-band radars which are not processed by international meteorological organizations, higher quality analyses can improve the description of finer scale features and thus provide better initial conditions to start limited-area weather models.