Impact of radar data assimilation for the simulation of a heavy rainfall case in central Italy using WRF – 3 DVAR

The aim of this study is to investigate the role of the assimilation of Doppler weather radar (DWR) data in a mesoscale model for the forecast of a heavy rainfall event that occurred in Italy in the urban area of Rome from 19 to 22 May 2008. For this purpose, radar reflectivity and radial velocity acquired from Monte Midia Doppler radar are assimilated into the Weather Research Forecasting (WRF) model, version 3.4.1. The general goal is to improve the quantitative precipitation forecasts (QPF): with this aim, several experiments are performed using the three-dimensional variational (3DVAR) technique. Moreover, sensitivity tests to outer loops are performed to include non-linearity in the observation operators. In order to identify the best initial conditions (ICs), statistical indicators such as forecast accuracy, frequency bias, false alarm rate and equitable threat score for the accumulated precipitation are used. The results show that the assimilation of DWR data has a large impact on both the position of convective cells and on the rainfall forecast of the analyzed event. A positive impact is also found if they are ingested together with conventional observations. Sensitivity to the use of two or three outer loops is also found if DWR data are assimilated together with conventional data.


Introduction
The initial conditions (ICs) are a key term for a successful forecast performed using a high-resolution numerical weather prediction (NWP).Generally, localized mesoscale features are not well represented in the ICs.The assimilation of local observations into the ICs may produce a forecast improvement.In the last 80 years, different methods of data assimilation have been investigated: the successive corrections method (SCM), the optimal interpolation (OI), the variational methods 3DVAR and 4DVAR, and the Kalman filter (KF) are some notable examples.During the last decade, high-resolution mesoscale models initialized using variational data assimilation techniques (3DVAR/4DVAR) are being increasingly applied to study meteorological phenomena (Kalnay, 2003).One of the reasons for the variational analysis becoming more and more popular is the ability to directly incorporate some non-conventional observations such as satellite radiance, radar reflectivity and radial velocity into numerical models (Kalnay, 2003;Barker et al., 2004), through the use of a proper operator.
Doppler weather radar (DWR) data may improve weather analyses and forecasts because of their high temporal and spatial resolution.In the last decade, high-resolution data together with a sophisticated technique of data assimilation have been chosen to improve the predictability of both convective cells and mesoscale convective systems.Furthermore, the assimilation of radar reflectivity and radial velocity has shown the potential for very short-range numerical weather prediction of rapidly developing convective systems.It is well known that reflectivity is related to the number of falling drops per unit volume, and it depends on the number and size of hydrometeors, whereas the vertical component of radial velocity contains information on vertical atmospheric motions.Both are important for the triggering and development of convection.However, radial velocity's contribution to its vertical component is variable, and it depends on the elevation of the radar antenna and variations in the atmospheric refractive index.The latter might produce variations in the radar ray paths with respect to those expected under standard atmosphere conditions.In this work, a standard atmosphere is assumed, and this means that radar ray paths propagate in a straight line (Bech et al., 2003).
The first radar data assimilation system for the storm scale was developed based on the four-dimensional variational data assimilation (4DVAR) technique and a boundary layer fluid dynamics model for the retrieval of the threedimensional wind and temperature (Sun et al., 1991).This system, known as VDRAS (the Variational Doppler Radar Analysis System), was later expanded to include microphysical retrieval, as well as short-term forecasts initialized by these retrieved fields (Sun andCrook, 1997, 1998).Another variational-based radar data assimilation system was developed by Gao et al. (2004) using a three-dimensional variational data assimilation (3DVAR) technique in the framework of the ARPS (Advanced Research and Prediction System, Xue et al., 2003) model.A so-called 3.5-dimensional variational radar data assimilation based on the navy's COAMPS (Coupled Ocean/Atmosphere Mesoscale Prediction System) was developed and verified through a number of studies (Zhao et al., 2006;Xu et al., 2010).These variational systems showed great potential in the use of radar observations for initializing high-resolution numerical models through several case studies and real-time demonstrations (Sun et al., 2010;Xue et al., 2010).
These results motivated the development of a radar data assimilation scheme in the WRFVAR variational data assimilation system of the ARW-WRF (Advanced Research Weather Research and Forecasting) community model.
The operators for radial velocity (Xiao et al., 2005) and reflectivity (Xiao et al., 2007;Xiao and Sun, 2007) data were added to the three-dimensional variational data assimilation system (Barker et al., 2003(Barker et al., , 2004;;Skamarock et al., 2008) developed at the National Center for Atmospheric Research (NCAR) laboratories for both the fifthgeneration Penn State/NCAR Mesoscale Model (MM5) and the Weather Research Forecasting (WRF) model.WRFVAR includes both 3DVAR and 4DVAR components.The radar DA scheme was first developed for WRF-3DVAR (Xiao et al., 2005(Xiao et al., , 2007) ) and recently expanded to 4DVAR (Wang et al., 2013;Sun and Wang, 2012).Xiao et al. (2005) assimilated radial velocities from a single Doppler radar into MM5 using the 3DVAR scheme for a heavy rainfall event.The vertical velocity increments were included via Richardson's balance equation, and an observation operator for radial velocity was developed.The results suggested that the scheme for the radial velocity assimilation is stable and robust in a cycling mode using high-frequency radar data.Moreover, continuous assimilation with 3 h update cycles was important for improving heavy rainfall.
A radar reflectivity data assimilation scheme was also developed within the MM5-3DVAR system, as described by Xiao et al. (2007).They showed that the intensity and track of Typhoon Rusa (2002) were better predicted through the combined assimilation of both radar radial velocity and reflectivity.Xiao and Sun (2007) used the WRF model and its 3DVAR data assimilation module to ingest data from a network of Doppler radars in order to study a squall-line convective system.They found that the assimilation of both radial velocity and reflectivity improved the quantitative precipitation forecast (QPF) skills with respect to using only one of the two radar variables.Further improvements to QPF were obtained by the assimilation of more than one radar site using the cycling procedure.
A few studies on the impact of radar data assimilation were also performed for the European area using different forecast models: AROME (Application of Research to Operations at MEsoscale) in France and HIRLAM (HIgh Resolution Limited Area Model) in Sweden.Montmerle and Faccani (2009) applied the 3DVAR assimilation technique for assimilating radial velocities from Doppler radars of the French ARAMIS (Application Radar la Météorologie In-fraSynoptique) network.They found a positive impact of radar wind data on the analyses and on precipitation forecasts.Lindskog et al. (2004) also assimilated radial wind velocity.They used the velocity-azimuth display (VAD) technique that provides vertical profiles of horizontal winds from the Doppler radar raw radial wind and the Doppler radar radial wind superobservations (SOs), which is based on spatial averaging.They found improvements in the wind and temperature forecasts of the low and middle troposphere using either VAD or SOs data.
The aim of this paper is to assess the impact of the assimilation of radar radial velocity and reflectivity data on the precipitation forecast, by using 3DVAR and WRF numerical models.For this purpose, a heavy rainfall case that occurred in central Italy is analyzed: the Aniene event (Rome, 19-22 May 2008).To explore the impact of radar data assimilation, several numerical experiments using different ICs are performed using the WRF model.A comparison between experiments with and without radar data is performed; moreover, to investigate the role of non-linearity in the observation operators, the number of outer loops is increased.Therefore, the novelty of this study resides in the assimilation of radar data in a complex orography area such as the Mediterranean one, where the forecast is still a challenging task.This paper is organized as follows.In Sect.2, the case study and the radar data used for the assimilation are described.A brief explanation of the WRF-3DVAR system and the radar data operator is presented in Sect.3. The 3DVAR experiments and the corresponding forecast results are discussed in Sects.4 and 5, respectively.The impact of outer loops is analyzed in Sect.6, whereas summary and conclusions are given in Sect.7.

A heavy rainfall case: the meteorological situation
From 19 to 22 May 2008, a heavy rainfall event occurred in the urban area of Rome.In the first hours of 19 May, a cyclonic circulation on the southern Mediterranean Sea (Fig. 1b) associated with an intrusion of cold air from Scandinavia (Fig. 1a) caused instability on the Italian peninsula.On 20 May, a deep low-pressure system developed in the Genoa Gulf (Fig. 1d), causing severe weather conditions.Southwesterly flow started to blow over the Tyrrhenian Sea, advecting warm and humid air toward the area of Rome (Fig. 1c), further destabilizing the atmosphere.
To understand the meteorological evolution of the event better, the surface rainfall intensity (SRI) product (in mm h −1 ) estimated by the Italian national radar network (Vulpiani et al., 2008) is shown in Fig. 2. In the early afternoon of 20 May, the precipitation is detected both in the urban area of Rome and over the Tyrrhenian Sea (RM, Fig. 2a).Later, local showers are observed in the southeastern part of RM and over the Tyrrhenian Sea (Fig. 2b).Moderate precipitation lasted until the evening around Rome (Fig. 2c), whereas at 22:00 UTC, the rainfall moved eastward, reaching L'Aquila (AQ, Fig. 2d).

Brief description of WRF-3DVAR
The WRF-3DVAR system is based on the MM5-3DVAR system (Barker et al., 2004); further developments and progress can be found in Skamarock et al. (2008).The WRF-3DVAR system is based on the multivariate incremental formulation (Courtier et al., 1994), where the preconditioned control variables are stream function ψ, velocity potential χ ,  unbalanced pressure p u and total water mixing ratio q t .The aim of the three-dimensional variational approach is to produce the best compromise between an a priori estimation of the analysis field and observations through the iterative solution that minimizes a cost function J .Most leading assimilation schemes do not perform the minimization process in the model space, but they use a transformed or control variable space that is the space allowed for the corrections to the background.This new space is chosen to have a special and desirable property when the background field is represented in this space: its errors are uncorrelated, and variances are of unit size (the problem is said to be preconditioned).
The cost function for 3DVAR is where x b is the generic variable of an a priori state (first guess), y o is the observation, and H is the operator that converts the model state to the observation space.This cost function J measures the distance of a field x from the observations y o and from the background x b : these distances are scaled through the matrices R and B, the observation and the background error covariance matrices, respectively.A correct evaluation of the error covariance matrices, both B and R, is crucial to a good-quality final analysis.The observation covariance error matrix R is usually diagonal, because the correlations between different instruments are assumed equal to zero.The background error correlation matrix B determines the spreading of information among observation locations, providing balanced information in data-void regions.That is why it plays an important role mainly in datasparse areas.A straightforward estimation of B is not possible, as the correlations between the variables have never been observed.Therefore, it has to be estimated using a statistical method, such as the National Meteorological Center (NMC) method (Parrish and Derber, 1992) or the ensemble one (Fisher et al., 1999).The first method is commonly used for the evaluation of B in the WRF-3DVAR system.It averages differences, valid at the same time t 0 , between two forecasts, one of them starting several hours later than the other (e.g., month-long series of 24 h minus 12 h forecasts valid at the same time).The only requirement is that the time t 0 is large enough to avoid any problem related to the model spin-up.A detailed description of the 3DVAR system can be found in Barker et al. (2003Barker et al. ( , 2004)).

Radar data assimilation methodology
Doppler radar data contain important high-resolution meteorological features.Radial velocity produces information on the atmospheric motions that is important for the onset of convection; moreover, it is well known that radar reflectivity is a measurement of the amount of precipitating hydrometeors (rain, snow, etc.).However, due to the complexity of coding, such variables are not included in most of the data assimilation schemes.
On the contrary, the WRF-3DVAR operator for the assimilation of Doppler radar data accounts for both reflectivity and the vertical velocity component of radial velocity.Vertical increments are estimated, including a balance equation based on Richardson (1922).This is the so-called linearized Richardson equation, which combines a continuity equation, an adiabatic thermodynamic equation, and a hydrostatic relation: it becomes important when radar data are included in the analysis.
Moreover, the total water q t is used as the moisture control variable instead of the pseudo relative humidity for the assimilation of radar reflectivity (Xiao et al., 2007).The "pseudo" relative humidity is the water vapor mixing ratio divided by its saturated value in the background state.To assimilate radar radial velocity and reflectivity associated with warm rain, vertical velocity and the microphysical parameters of cloud water and rainwater must be added as control variables.It is well known that radar observations can be affected by several sources of errors, mainly due to ground clutter, attenuation and radio interferences.Particularly, weather radar operating in complex orography may be affected by a significant beam blockage that can strongly degrade the monitoring capabilities and accordingly the rainfall estimation at the ground.
For this reason, a preliminary procedure to correct acquired radar data is applied before the assimilation proce-  filter.Partial beam blockage is corrected by adopting a compensating technique (Fulton et al., 1998), while attenuation is mitigated by means of rain path-integrated attenuation (PIA) techniques (Picciotti et al., 2006).Once corrected, both reflectivity and velocity data are ready for ingestion into 3DVAR.Typical errors of 5.0 m s −1 for radial velocity and 5.0 dBZ for reflectivity are assumed.The assimilation window has been set to ±5 min at the analysis time.Finally, a total number of 344 160 radar observations were ingested into the model.

Radar observation operators: radial velocity and reflectivity
To assimilate radar observations, the following observation terms are added to the existing cost function: where J old is used to represent the existing cost function before radar data assimilation is developed.The variables V r stand for the radial velocity and Z for the reflectivity factor.The superscript "o" indicates the observations.The symbols R −1 v and R −1 z stand for observation error covariance matrices for radial velocity and reflectivity, respectively (Sun and Wang, 2013).Note that the summation over the observation time levels k is not needed in the case of 3DVAR.As explained previously, the observation operator H in the cost function (1) links the model variables in a model coordinate to the observation variables in an observation space.For the radar radial velocity, this linkage is formulated with the three-dimensional wind field (u, v, and w), the hydrometeor fall speed or terminal velocity v t , and the distance r between the location of a data point and the radar antenna: where (x, y, z) represents the location of the observation point and (x i , y i , z i ) represents the location of the radar station.
Following the algorithm of Sun and Crook (1998), the terminal velocity is given by where "a" is a correction factor defined as follows: where p is the base-state pressure and p 0 is the pressure at the ground.
The formulation of the reflectivity operator is not as straightforward, because it depends on the assumption of drop size distribution in a microphysical parameterization scheme and the classification of hydrometeors.To assimilate radar reflectivity directly, the total water mixing ratio q t was chosen as a control variable, and a warm rain process was introduced (Dudhia, 1989) into the WRF-3DVAR system.This allowed the increments of moist variables linked to the hydrometeors to be produced, such as the water vapor mixing ratio q v , the cloud water mixing ratio q c , and the rainwater mixing ratio q r .Once the 3DVAR system produces the q c and q r increments, the setup of the observation operator for the assimilation of reflectivity is obtained.
Following Sun and Crook (1997), Xiao et al. (2007) used the following relation for WRF-3DVAR: where Z is the radar reflectivity expressed in dBZ, ρ is the air density (kg m −3 ) and q r is the rain water mixing ratio (g kg −1 ).This relation was derived analytically by assuming a Marshall-Palmer distribution of raindrop size with intercept N 0 = 8×10 6 mm −4 , and no contribution of ice phases to the reflectivity is accounted for.Equation ( 6) was used in cost function (2) to assimilate reflectivity in the WRF-3VAR developed by Xiao et al. (2007).

Model setup
WRF model version 3.4.1 (Wang et al., 2011) is used for this study.Two domains (D0i), running independently, are used (Fig. 4).The outer domain (D01) has a resolution of 12 km with a number of horizontal grid points of 263 × 185.It is initialized using ECMWF analysis, and the boundary conditions are upgraded every 6 h.The inner domain (D02) has a horizontal grid spacing of 3 km (445 × 449), and it is initialized using WRF-ARW output at 12 km.37 eta vertical levels are used for all domains.The model top is at 50 hPa.The new Thompson et al. scheme (Thompson et al., 2008) for microphysics, the Kain-Fritsch scheme (Kain, 2004) for cumulus convection, for D01 only, the Mellor, Yamada, Nakanishi and Niino level 2.5 (Nakanishi and Niino, 2006) for planetary boundary layers, the rapid radiative transfer model (Mlawer et al., 1997) and the Goddard (Chou and Suarez, 1994) schemes for longwave and shortwave radiation, respectively, are used for all the experiments.The Noah land surface model (LSM; the successor to the OSU (Oregon State University) LSM described by Chen and Dudhia, 2001) for land surfaces is used.
Besides the first guess x b and the observations y o , as explained in Sect.3.1, another important input for WRF-3DVAR is the background error covariance matrix B. This is computed using the NMC method, and a domain-specific B has been generated (CV5 option) using the empirical orthogonal function (EOF) to represent the vertical covariance.
To estimate the NMC-based error statistics, two forecasts have been performed every day for a period of 1 week (Sugimoto et al., 2009) starting on 15 May 2008: a 24 h forecast (starting at 00:00 UTC) and a 12 h forecast (starting at 12:00 UTC) valid at the same time.The differences between the two forecasts at t +24 and t +12 are used to calculate the domain-averaged error statistics.

Experimental design
Several experiments are performed with the aim of improving ICs.The GTS (Global Telecommunication System) conventional observations -SYNOP (surface synoptic observations) and TEMP (upper-level temperature, humidity and winds) -are used with and without non-conventional radar data.The assimilation window for conventional data is set to ±1 h, whereas the one for radar data is set to ±5 min at the analysis time.To verify the impact of radar data, a total of four WRF experiments are performed for the highresolution domain (3 km): in the control experiment (CTL), no data assimilation is performed; in the CON experiment, the assimilation of SYNOP and TEMP observations is performed at 00:00 UTC, 20 May 2008; in the RAD experiment, the assimilation of Monte Midia radar data (both reflectivity and radial velocity) is performed; in the ALL experiment, both conventional data (SYNOP and TEMP) and radar data are ingested.A specific background error covariance matrix calculated for the highest-resolution domain is used for CON, RAD and ALL.For these above simulations, one outer loop is used as the default.Two additional experiments (ALL_2OL and ALL_3OL) are carried out similarly to ALL, but using two and three outer loops, respectively, during the assimilation procedure.The multiple outer loop strategy (Rizvi et al., 2008) allows for the inclusion of nonlinearity in the observation operators and for the assessment of the influence of observations entering for each cycle.The non-linear problem is solved iteratively as a sequence of linear problems by running more than one analysis outer loop, so the assimilation system is able to ingest more observations.

Initial conditions
To evaluate the impact of 3DVAR on the initial conditions, a comparison between the model and the observed vertical sounding at Pratica di Mare (PDM) at 00:00 UTC on 20 May 2008 is performed.This comparison is carried out using a few meteorological indexes, such as convective  available potential energy (CAPE), useful for assessing the instability of the atmosphere and therefore the possibility of triggering convection, a lifted index (LI), which measures the severity of the thunderstorm, and the K index, which gives an indication of the thunderstorm occurrence.The CAPE at Pratica di Mare suggests a weakly unstable environment and that severe thunderstorms are unlikely to occur (small LI).No significant differences are found between the experiments, except for the LI listed in Table 2.The CAPE is significantly underestimated for all simulations, whereas the LI is overestimated.The WRF K index is in good agreement with the observed one.
The observed and WRF radio soundings at Pratica Di Mare at 00:00 UTC, 20 May 2008 (Figs. 5 and 6, respectively) clearly show no large differences between the WRF simulations (CTL, CON, RAD and ALL) below 500hPa.On the other hand, the assimilation of conventional data slightly improves the vertical profile of the dew point temperature.Both RAD and ALL produce a slightly dryer PBL than the observed one.The sharp variation in the dew point content as well as the small temperature inversion between 600 and 300 hPa are completely missed by all the WRF simulations.The strong wind shear between 500 and 400 hPa and the one at 300 hPa are both underestimated by WRF.Finally, above 300 hPa, a still drier atmosphere than the one observed is produced by WRF.

CTL
No data assimilation, using WRF-ARW at 12 km to initialize the highest-resolution domain at 3 km.CON SYNOP and TEMP data assimilation at 00:00 UTC on 20 May 2008 at 3 km resolution, using a previous WRF forecast lasting 24 h as FG.RAD Radar data assimilation at 00:00 UTC on 20 May 2008 at 3 km resolution, using a previous WRF forecast lasting 24 h as FG.ALL Conventional and radar data assimilation at 00:00 UTC on 20 May 2008 at 3 km resolution, using a previous WRF forecast lasting 24 h as FG.

Impact on radar reflectivity and rainfall
The impact of Monte Midia radar data assimilation on the forecast is evaluated by comparing the observed reflectivity and the estimated rain rate with those produced by the WRF experiments.

Impact on the precipitation forecast
The VMI (vertical maximum intensity) reflectivity of the Monte Midia radar at 14:00 UTC on 20 May clearly shows the precipitation in central Italy (Fig. 7a) and on the seaside near the Lazio coast.Precipitation can also be inferred from the cloudiness derived from the Geostationary MSG (Meteosat Second Generation) satellite images (Fig. 7b, c), where a "banana-shaped" cloud structure is also clearly shown over the central Tyrrhenian Sea.Aligned cells are also detected inland near Rieti (Fig. 7a, red circle).The Monte Midia radar highlights this last feature: a reflectivity of approximately 50-55 dBZ at 14:00 UTC associated with a few small convective cells oriented north-south between Rome (RM) and Rieti (RI) (solid red oval in Fig. 7a) is clearly detected by the radar.Another area of moderate precipitation (about 40 dBZ on the border between Lazio and Abruzzo) was also detected (red arrow in Fig. 7a).
The CTL simulation clearly shows (Fig. 8a) a large area of reflectivity over the sea that is overestimated both in amount and the affected area.Moreover, an attempt to reproduce the small convective cells between Rome and Rieti is also found, but these small cells are wrongly oriented and the reflectivity is underestimated.The assimilation of conventional data (Fig. 8b) removes part of the north-south oriented area of reflectivity north of Viterbo (VT), making it closer to the observed one.The CON simulation produces an enlargement of the area covered by the cells near Rieti, but they are not yet correctly aligned.The RAD simulation (Fig. 8c) partly reduces the areal distribution of the reflectivity over the seaside, showing an attempt to separate the two elongated areas of reflectivity as they can be deduced by the satellite image (Fig. 7b, c).Moreover, RAD both reduces the convection in the Viterbo area and correctly generates aligned cells north of Rome near Rieti (Fig. 8c, red oval), as the Monte Midia radar detected (Fig. 7a).The assimilation of radar and conventional data (ALL, Fig. 8d) shows a distribution of the reflectivity over the Tyrrhenian Sea even closer to the one deduced by the satellite images than by RAD, but an over-reduction of the structure and reflectivity of the cells north of Rieti is also found.The previous analysis suggests that the assimilation of radar reflectivity and radial velocity largely impacts the forecast even if only one radar is assimilated.Moreover, the assimilation of either the radar data only or the radar and convectional data together reduces the differences between the forecast and the observation.
Similar to the reflectivity, a comparison between the observed (Fig. 9) and the WRF 12 h accumulated rainfall (Fig. 10) is now performed.
The estimated surface rainfall intensity (mm h −1 ) by the national radar network in central Italy on 20 May 2008 (Fig. 2) clearly shows precipitation over the sea near the coast of Lazio.The "banana-shaped" structure detected by the satellite is clearly reproduced by the rainfall structure at the end of the day (Fig. 2d).The control simulation (CTL, Fig. 10a) correctly reproduces the structure of the precipitation over the sea, but clearly shows an overestimation of the precipitation both over the sea and inland on the western side of the Apennines ridge.On the other hand, the precipitation is correctly restrained on the western side of the Apennines.The assimilation of conventional data (CON, Fig. 10b) clearly corrects the precipitation forecast by producing a reduction in both the areal distribution and the amount of the precipitation, both on the sea and inland.Moreover, an increase in the precipitation in the urban area of Rome is found.The RAD simulation (Fig. 10c) correctly produces a reduction in the accumulated precipitation over the sea and in the urban area of Rome; on the other hand, no large differences are found for the precipitation inland between CTL, CON and RAD.This last one is the sole simulation able to produce a red spot in the Rieti area (Fig. 10c) like the observed one (Fig. 9) clearly related to the aligned cells discussed previously.Similarly to CTL and CON, ALL (Fig. 10d) does not produce the red spot of precipitation in the Rieti area, but the precipitation between RM and FR is correctly reproduced, as for CTL and RAD.
In summary, the precipitation forecast also clearly shows an impact of the radar data assimilation by producing a precipitation pattern closer to the observed one.

Statistical evaluation
In order to compare the experiments carried out for the Aniene event objectively, four statistical indicators are used (Wilks, 1995): forecast accuracy (ACC), false alarm rate (FAR), frequency bias (FBIAS) and equitable threat score (ETS).
The ACC index shows the accuracy of the forecast; a perfect forecast has ACC = 1.FAR estimates the forecast frequency failures; a perfect forecast has FAR = 0. FBIAS gives information on the correctness of a precipitation forecast: values greater than 1 indicate an overestimation in the number of forecast events; the perfect value is FBIAS = 1.ETS represents the fraction of the events reproduced correctly, taking into account random hit chance; the best value is 1.The ETS index may have values lower than or near zero if the forecasts are wrong, so values of 0.5/0.6 are fairly good.
Figure 11 shows the results for the previous indexes for the 12 h accumulated rainfall as a function of different thresholds.The WRF simulations show no large differences between them, but ACC clearly shows improvements for CON (Fig. 11, green line) with respect to other experiments for thresholds between 5 and 30 mm/12 h; accordingly, the FAR index for CON shows values lower than the others for the same thresholds.Moreover, CON produces the best values for ETS at all thresholds, even though it rapidly degrades for the higher ones.Vice versa, FBIAS gives the worst values for CON between 5 and 20 mm/12 h, whereas it improves for higher thresholds between 30 and 40 mm/12 h.Moreover, the values of FBIAS are close to the best value of 1 for thresholds lower than 30 mm/12 h for all the simulations, whereas this index has values greater than 1 for all the experiments.
These results do not support completely the previous finding that the RAD and ALL experiments produce a better forecast; this disagreement is probably produced by the rain gauge location.In fact, the highest density of surface stations is nearby Rome (Fig. 3), where the CON simulation produces rainfall in good agreement with the observation.

Impact of outer loops on precipitation
Based on the previous results, a set of sensitivity tests for WRF-3DVAR is performed for the ALL experiment only.These new experiments are carried out using two and three outer loops during the assimilation procedure (Table 3).All the experiments performed using 3DVAR (CON, RAD and ALL) discussed previously are carried out using one outer loop.The multiple outer loop strategy (Rizvi et al., 2008) allows for the inclusion of non-linearity in the observation operators and for the assessment of the influence of observations entering for each cycle.The non-linear problem is solved iteratively as a sequence of linear problems by running more than one analysis outer loop, so the assimilation system is able to ingest more observations.
As already pointed out in Sect.5.2.1, ALL (Fig. 8d) produces a reduction in the differences between the forecast and the observed reflectivity.Using two outer loops during the assimilation procedure, a strong areal attenuation of the dBZ over Viterbo is obtained (Fig. 12a), but at the same time, the aligned cells near Rieti are not reproduced exactly.Adding one more loop (ALL_3OL, Fig. 12b), the reflectivity near Viterbo is increased again, whereas the cells north of Rieti are reproduced better: they are similar to what ALL produces and to what the Monte Midia radar detected.Finally, the "banana-shaped" structure over the central Tyrrhenian Sea is well defined.These results suggest that the increase in the number of outer loops may positively impact the forecast for heavy rainfall.
To evaluate the impact of outer loops better, ACC, FAR, FBIAS and ETS are used again.From the ACC index graph, no large differences are highlighted between the three experiments.On the other hand, FBIAS shows small differences for thresholds between 10 and 30 mm/12 h for ALL_2OL and ALL_3OL (Fig. 13, upper right, green and violet curves, respectively): both have values close to 1.The ETS index produces better scores for ALL_3OL (violet curve, Fig. 13, bottom left) for thresholds lower than 30 mm/12 h.Finally, FAR gives a quite good response for the two experiments if using more than one outer loop (green and violet curves, Fig. 13, bottom right) for very low thresholds between 0 and 5 mm/12 h.
The previous results confirm that the increase in the number of outer loops may positively impact the forecast for heavy rainfall.

Summary and conclusions
This study investigates the impact of radar data assimilation on the WRF forecast of a heavy rainfall event that occurred in the urban area of Rome.In order to improve the model's initial conditions, high-resolution information is used: radar radial velocity and reflectivity and conventional observations are assimilated using the WRF-3DVAR scheme.Several experiments, including sensitivity tests to multiple outer loops, are performed to assess the impact of radar data assimilation.The comparison between the observed and model reflectivity allows the assessment of a better performance of both the simulation-assimilating radar data alone (RAD) and together with conventional observations (ALL).If using more than one outer loop, for the ALL experiment only, a positive impact is obtained for two and three outer loops (ALL_2OL and ALL_3OL, respectively).Statistical indicators are used to compare the model simulations objectively.The statistical results do not support the previous finding, mostly because the surface stations are concentrated near the urban area of Rome.Obviously, the assimilation of these stations (CON) drives the results toward them.Statistical indicators are also used to investigate the outer loop impact: in this case, statistics confirm the model results.This is because the multiple outer loop technique allows the assimilation of a larger number of observations progressively into WRF-3DVAR.
In conclusion, the assimilation of radar radial velocity and reflectivity data shows a positive impact on precipitation forecasts, also if they are ingested together with conventional observations and if the outer loop strategy is used.
Finally, to evaluate the previous results properly, it has also to be considered that the WRF-3DVAR radar data assimilation scheme has a few limitations: 1. the lack of the ice phase in the 3DVAR microphysics; 2. the estimation of the B matrix for high-resolution data assimilation can be improved by using the ensemble approach or by the combined use of an ensemble Kalman filter.
Moreover, the lack of radar coverage on the Tyrrhenian coastline from where the system approaches it is also a limiting factor in the assimilation process.Finally, the technique of multiple outer loops will be investigated further by tuning length-scale and observation error parameters, and thinning of radar data has to be undertaken either to reduce the observation-error spatial correlation or the computational cost of the assimilation (Montmerle and Faccani, 2009).Moreover, following Wang et al. (2013), the new observation operator that, instead of assimilating radar reflectivity directly, assimilates retrieved rainwater and water vapor derived from radar reflectivity, could be applied.Wang et al. (2013) showed that one of the problems of assimilating the reflectivity is in the use of a linearized Z − q r (reflectivity-rainwater) equation as the observation operator and of the warm-rain partition scheme.On the other hand, the assimilation of the estimated water vapor is expected to provide a favorable environment that supports convection.
In the future, additional case studies will be carried out to quantify the percentage of success of WRF-3DVAR in assimilating DWR data to improve precipitation forecasts, as well as using data from several operative radars located in central Italy, also including dual-polarization systems.

Figure 1 .
Figure 1.European Center for Medium-Range Weather Forecasts (ECMWF) analyses at 06:00 UTC on 19 May 2008: (a) temperature and geopotential height, at 500 hPa, and (b) mean sea level pressure and surface wind in m s −1 ; ECMWF analyses at 06:00 UTC on 20 May 2008: (c) temperature and geopotential height, at 500 hPa, and (d) mean sea level pressure and surface wind in m s −1 .

Figure 2 .
Figure 2.Estimated surface rainfall intensity (mm h) by the national radar network in central Italy on 20 May 2008: (a) SRI at 14:00 UTC; (b) SRI at 16:00 UTC; (c) SRI at 20:00 UTC; (d) SRI at 22:00 UTC.Only the highlighted areas were covered by the national radar network in 2008 (courtesy of the Italian Civil Protection Department).

Figure 3 .
Figure 3.The panel shows the location of the Monte Midia radar between the Abruzzo (to the right) and Lazio (to the left) regions, its coverage (gray circle), and the rain gauge network (red crosses) used for calibration of the estimated rain.

Figure 4 .Figure 5 .
Figure 4. Model configuration for WRF simulations using two domains: (a) overlapping of 12 and 3 km; (b) the high-resolution D02 ( x = 3 km) includes the Monte Midia radar (red dot in the figure).

Figure 7 .
Figure 7. (a) VMI (dBZ) from the Monte Midia radar at 14:00 UTC, 20 May 2008 (courtesy of Centro Funzionale of the Abruzzo region).The solid red circle and the red arrow indicate the precipitation patterns selected for the analysis.(b) Geostationary MSG map taken at 14:00 UTC at the visible channel.(c) Geostationary MSG map taken at 14:00 UTC at the infrared channel.The "banana-shaped" cloud structure is clearly evident over the central Tyrrhenian Sea.

Figure 11 .
Figure 11.ACC (top left), FBIAS (top right), ETS (bottom left) and FAR (bottom right) as a function of threshold.The red line shows CTL, the green line CON, the violet line RAD, and the cyan line ALL.

Figure 13 .
Figure 13.Outer loop sensitivity: ACC, FBIAS, ETS and FAR for ALL, ALL_2OL and ALL_3OL.The red line refers to one loop, the green line to two loops, and the violet line to three loops.

Table 1 .
List of experiments and related description.

Table 2 .
Values of the CAPE, LI, and K indexes for PDM sounding and those produced by the model at the same time (00:00 UTC, 20 May 2008) and location (41 • 65 N, 12 • 43 E).

Table 3 .
Description of the experiments performed using 3DVAR with multiple outer loops.Expt DescriptionALL_2OL As expt ALL, setting 2 outer loops during the assimilation procedure.ALL_3OL As expt ALL, setting 3 outer loops during the assimilation procedure.