Investigating 3D and 4D Variational Rapid-Update-Cycling Assimilation of Weather Radar Reflectivity for a Heavy Rain Event in Central Italy

The precipitation forecast over the Mediterranean basin is still a challenge because of the complex orographic region which amplifies the need for local observation to correctly initialize the forecast. In this context the data assimilation techniques 10 play a key role in improving the initial conditions and consequently the timing and position of precipitation forecast. For the first time, the ability of a cycling 4D-Var to reproduce a heavy rain event in central Italy, as well as to provide a comparison with the largely used cycling 3D-Var, is evaluated in this study. The radar reflectivity measured by the Italian ground radar network is assimilated in the WRF model to simulate an event occurred on May 3, 2018 in central Italy. In order to evaluate the impact of data assimilation, several simulations are objectively compared by means of a Fraction Skill Score (FSS), which 15 is calculated for several threshold values, and a Receiver Operating Characteristic (ROC) curve. The results suggest that both assimilation methods in cycling mode improve the 1, 3 and 6-hourly quantitative precipitation estimation. More specifically, the cycling 4D-Var with a warm start initialization shows the highest FSS values in the first hours of simulation both with light and heavy precipitation. Finally, the ROC curve confirms the benefit of 4D-Var: the area under the curve is 0.91 compared to the 0.88 of control experiment without data assimilation. 20


Introduction
Nowadays, the high-resolution rainfall forecast from Numerical Weather Prediction (NWP) models is essential for several applications. It is used by civil protection agencies to contrast the hydrological risks and safeguard people during severe weather event, by disaster management agencies to prepare emergency interventions as well as by public event managers, private enterprises, and common people to plan their daily activities. Recently, the development of more accurate 25 parametrization of physical processes allowed a significant progress in NWP at high resolution but the prediction of exact location, timing and intensity of a convective event is still a challenge (Stensrud et al., 2009;Yano et al., 2018;Mass et al., 2002;Torcasio et al., 2021).
The Mediterranean basin is prone to flash flood and heavy precipitation events. One of the most relevant peculiarities of this area is the presence of mountain ranges in proximity of coastal areas that lift the airflow, favouring the condensation and the 30 triggering of convection. The Italian territory, which is characterized by a complex orography with two relevant mountain chains (Alpes and Apennines) and steep, urbanized, small catchments, is one of the most exposed Mediterranean area to the hydrogeological risks. About 90% of Italian municipalities are susceptible to floods which have caused 466 deaths between 1990 and 2006 alone and over 19 billion of Euro of damages (Llasat et al., 2010).
The Mediterranean region has also been identified as hotspot for climate change because of its high sensitivity to Global 35 Greenhouse Gas (GHG) concentrations (Giorgi, 2006). The heat stress growing, expected in future years, will contribute to a temperature rise, as well as to increase the water vapour amount in the atmosphere. These factors will lead to an intensification of intense and extreme precipitation events (Willet et al., 2008;Giorgi, 2011). In this context, an accurate prediction and observation of rainfall will be crucial to prevent economic and social damages. The precipitation forecasts from the Numerical Weather Prediction (NWP) models is strongly dependent on the initial state, that dominates the evolution of the prognostic 40 variables and consequently the development of precipitation. In NWP model convective rain events are usually associated with low forecast probabilities due to the high spatial variability of precipitation and uncertainties in convective initiation, hence, small errors in the initial conditions produce a significant shift in the position and intensity of convective events.
Nowadays, the large availability of high frequency (both in space and time) meteorological data, remote sensing observations and in situ measurements, has encouraged many operational centres to use data assimilation techniques for improving the 45 accuracy of initial state. More specifically, the assimilation of ground radar reflectivity and radial velocity with threedimensional variational (3D-Var) method proves good results in terms of quantitative precipitation forecast (QPF) for several case studiesy in United States and Korea (Xiao and Sun, 2007;Lee et al., 2010;Ha et al., 2011). The assimilation of radar data with 3D-Var confirms positive results also in Europe by using Advanced Regional Prediction System (ARPS) and Application of Research to Operations at Mesoscale (AROME) models for two heavy rainfall cases in Croatia and France (Stanesic and 50 Brewster, 2016;Caumont et al. 2009) respectively, as well as by using WRF model (Schwitalla et al., 2014), for the simulation of a convective event during the Convective and Orographically induced Precipitation Study (COPS).
The benefits of assimilation of radar reflectivity are also confirmed for two flash flood events in central (Maiello et al., 2014; and northern Italy (Lagasio et al., 2019) with WRF model, and in combination with lightning (Federico et., 2019;Torcasio et al., 2021). Gastaldo et al. (2018), instead, assimilated reflectivity volumes using a local ensemble transform Kalman 55 filter (LETKF) implemented in the convection-permitting model of the COnsortium for Small-scale MOdelling (COSMO) operating at Hydro-Meteo-Climate Service of the Emilia-Romagna Region (Arpae-SIMC) pointing out the positive impact of radar assimilation in quantitative precipitation forecast (QPF) accuracy both when a latent heat nudging (LHN) is applied or not. Finally, a new work made by Gastaldo et al. (2021) confirms the positive impact, up to 7 hours, of radar assimilation with LETKF in COSMO model, especially in convective cases, replacing the LHN scheme. Recently, a few works showed a positive 60 impact in the prediction of intense precipitations using four-dimensional variational (4D-Var) with radar data and conventional observations if compared to 3D-Var;, they concern the simulation of a cyclonic event in Antarctic region (Chu et al., 2013) and a squall line over US Great Plains (Sun and Wang, 2013), respectively. In Europe, the 4D-Var method proved good performance, improving the QPF and reducing the precipitation spinup time (Mazzarella et al., 2017;, but due to its high computational cost, it is scarcely applied except in big operational centres. 65 The low predictability and the high spatial and temporal variability of convective precipitation pattern requires a rapid upd ate of initial state (analysis) to reduce the errors and to ensure a well-balanced and physically consistent initial conditions. Newly, significant efforts have been made by scientific community to improve the temporal and spatial resolution of remote sensing data by weather radar or satellite-borne sensors, and this has enabled the first experiments with an update frequency equal or less than 3 hours. In this context, several weather centres have adopted a cycling assimilation with a 3-hourly update frequency 70 by using 3D-Var with promising results (Ballard et al., 2012;Gao et al., 2004;Stephan et al., 2008;Wang et al. 2013a;Caumont et al. 2009). Only the High-Resolution Rapid Refresh (HRRR) system, developed by NCAR, assimilates radar data with subhourly frequency over USA, but it does not use a variational method (Smith et al., 2008).
The cycling assimilation with 4D-Var is still a challenge because of the high demand of computational resources. A first attempt to apply the cycling 4D-Var in operational mode was made during the London Olympic Game in 2012 using an NWP-75 based nowcast system (Ballard et al., 2016). This work shows the advantages of 4D-Var that ingests more observations than 3D-Var, estimating with a greater accuracy the initial state of the atmosphere. In addition, the weather radar reflectivity over the whole Italian territory, previously used in a LETKF assimilation scheme (Gastaldo et al., 2018) with promising results, are now assimilated for the two variational methods.
This study aims at: i) assessing the performance of 1-hour cycling assimilation with 3D-Var and 4D-Var methods by using 80 WRF model; ii) evaluating the impact of radar reflectivity mosaic, acquired by the Italian radar network, in cycling assimilation with variational methods; iii) quantifying the improvements of assimilation techniques in terms of QPF for a complex orography region in the Mediterranean basin. In this regard, a heavy rain case occurred in Central Italy on May 3 2018 is used and several experiments are carried out to quantify the impact of the two assimilation methods in cycling mode. WithTo the aim of identifying the best configuration in terms of QPF, two different statistical methods are applied: a filtering 85 (neighbourhood) technique and a Receiver Operating Characteristic (ROC) statistical indicator. The keys novelty of this paper resides in: i) the application of hourly cycling 4D-Var assimilation with WRF model; ii) the comparison between the two variational assimilation methods 3D/4D-Var in cycling mode; iii) the assessment of the precipitation forecast skill of cycling 34D/43D-Var assimilation methods in an orographically complex region.
The paper is organized as follows. The case study is described in Section 2, while the assimilated dataset is presented in Section 90 3. Section 4 contains the methodological approach: more specifically a brief overview of the assimilation and verification methods, while the model's setup and the design of experiments are described in Section 5. The results are summarized in Section 6, lastly, the conclusions are discussed in Section 7.
On 2 May 2018, the synoptic scenario was characterized by an upper level through at tear off stage extended from Central 95 Europe to northern Africa whose evolution is slowed by the presence of two anticyclonic circulations over Atlantic Ocean and Russia respectively (Fig. 1). The upper level through, with a northwest-southeast oriented axis, brought cool and dry air from the Artic region towards the Mediterranean basin, and advected a moist and warm south-easterly flow over the Adriatic region.
In the following 24 hours the upper level through evolved in a cut off low, producing a surfacedeep low-pressure system (992 hPa) on the western side of Sicily (Fig. 2a). The surface depression slowly moved north-eastward, dissipating its energy over 100 the next 12 hours.

105
The mesoscale conditions over Italian peninsula showed a strong and moist south-easterly flow at upper and middle atmospheric levels, whereas a convergency line between the north-eastern (Bora) and south-eastern (Sirocco) winds developed at low levels ( Fig. 2b). This produced heavy precipitations on May 3 over the central Italy Adriatic region which was enhanced by the complex orography of this area. Indeed, the highest peaks of the Apennines mountains (Gran Sasso and Majella, 2912 m and 2793 m MSL, respectively) near the coast further increased the atmospheric instability, promoting the triggering of 110 convective cells.
the rain gauge station located in Fano a Corno, which measured 152 mm. The precipitation also affected the coastal sector with amounts of about 73 mm at Pescara and Pineto and 42 mm at Vasto (Fig. 3).

Datasets Description
To the aim of assessing the impact of cycling 3D-Var and 4D-Var, the composite radar reflectivity provided by the Italian radar network (Vulpiani et al., 2008) are assimilated in WRF for the numerical simulations. 125 The Italian ground radar network includes 20 C-band and 3 X-band radar, managed by 11 regional administrations. Among these, 7 C-band and the 3 X-band systems (all with dual-polarization capability) are managed directly by the Italian Civil Protection Department (DPC), which is also the developer and distributor of the national precipitation product. The processing architecture is basically composed of two main steps: firstly, the radar measurements are locally processed by a unique softw are system, secondly all the products are centralized to generate the national level products. 130 Different sources of error can affect the radar measurements (Collier, 1996): non-weather returns (clutter), partial beam blocking, beam broadening at increasing distances, vertical variability of precipitation, Radio Local Area Network (RLAN) interferences and rain path attenuation. Due to the morphology of the Italian territory, the uncertainty can be mainly associated to the orography-related effects, especially in southern Italy where the radar coverage as well as the radar overlapping is poor. The DPC processing system aims at identifying most of the uncertainty sources in order to compensate them, whenever it is 135 possible, before estimating precipitation. A point-by-point description of the operational radar processing chain can be found in (Petracca et al., 2018). After the processing, some composite products are generated in real time by DPC and disseminated at national and regional level. Among these, the reflectivity Constant Altitude Plan Position Indicator (CAPPI) at 2000, 3000 and 5000 m MSL, which cover the whole Italian territory, are assimilated into WRF model.
It is worth mentioning that the CAPPI gives a horizontal cross-section of the data at constant altitude, it is a two-dimensional 140 areal representation extracted from three-dimensional radar volume scan data.
Moreover, a thinning is applied to the CAPPI reflectivity data to ensure uncorrelated observation errors (Chang et al., 2014;Liu and Rabier, 2003) and to reduce the computational complexity, especially for 4D-Var method. CAPPI data are thinned over the 3km domain grid (described in section 5.2) by using an 'ad hoc' procedure. The rainfall data to assess the impact of 3D/4D-Var methods in cycling mode are provided by the rain-gauge network of the Italian Civil Protection Department (DPC) 145 composed by roughly 3000 tipping bucket gauges with the minimum detectable rain of 0.2 mm. A careful quality check is applied to the data before using them in the statistical analysis (Hanachi et al., 2014). More specifically, the following actions are carried out: • control of rain gauges with the same name but different coordinates; • removal of data associated with rain gauges without valid coordinates; 150 • removal of duplicate data; • identification of anomalous data (for example very different values respect to the surrounding rain gauges).

Variational Data Assimilation
The Weather Research and Forecasting (WRF) model (Skamarock et al., 2019) is used for the numerical simulations, while the CAPPI radar data are assimilated by using the WRF data assimilation system (WRFDA, Barker et al., 2012). This system 155 contains the algorithms required by 3D/4D-Var variational assimilation methods which are described in the following subsections. In this study the 4D-Var is employed in cycling mode to reduce the computational cost and compared with the widely used 3D-Var again with a cycling update.

3D-Var and 4D-Var methods 160
Nowadays, the WRF 3D-Var (Barker et al., 2012) variational assimilation method is widely applied in meteorological and oceanographic modelling to assimilate a large variety of observations and to generate reliable initial conditions. The 3D -Var is a mathematical technique that combines observations and a short-range forecast (first guess) through the minimization of a cost function. The goal of this method is to reduce the misfit between the observation and the background fields to obtain the best estimate of the true state of the atmosphere. In general, the cost function with respect an atmospheric state variable x is 165 defined as follow: where B and R are the background and observation error matrices, respectively; y0 represents the observations, x0 the background field or first guess and finally H is the forward observation operator that converts data from model space to observation space. 170 The 3D-Var method has the advantage of being computationally cheap even if it misses the time dependency, hence all the observations that are acquired during the assimilation window are considered at its central time.
The WRF 4D-Var (Huang et al., 2009) take the time variable into account using a numerical weather forecast as dynamical constraint. More specifically, the method computes the model trajectory that reduce the misfit with the observations distribu ted in the assimilation window. The initial atmospheric state is determined by minimizing the following cost function (2): 175 The assimilation window is divided in k discrete sub windows, where x0 is the analysis state vector at time k0; Hk and Mk are the nonlinear forward and observation models, respectively. Finally, B and R are, as already mentioned, the background and observation error matrices.
The 4D-Var has the advantage of assimilating the observation at its exact time compared to 3D-Var, ensuring a more consistent 180 physics and dynamical balance of the initial conditions. Given the high computational complexity, the incremental approach proposed by Courtier et al. (1994) has been adopted. The tangent linear and adjoint model with a simplified physics are used in the inner loop to increase the computational efficiency. Despite that the application of 4D-Var in operational mode is still limited to the major weather centres only.
The B matrix is a key component in the assimilation processes because it weights the errors in the background field adjusting 185 the spatial spreading of observational information. In this context, the National Meteorological Center (NMC) method (Parrish and Derber, 1992), widely use in the data assimilation community, has been adopted for this work. The B matrix is estimated evaluating the difference between forecasts valid at the same time, but one of them is initialized 12 hours after the other. More in detail, the new method (Wang et al., 2013b) considers u, v, temperature, pseudo relative humidity, and surface pressure as control variables, in contrast with the old one which utilizes the stream function, the unbalanced velocity potential and the 190 unbalanced temperature. Recent works using the new method, shows a slightly benefit in the precipitation forecast skill as well as an improving of performance when radar data are assimilated (Wang et al., 2013a(Wang et al., , 2013cSun et al., 2016). Thus, for this work the B matrix is computed over a period of 2 weeks from 1 May to 15 May 2018 by using this new method. According to the values provided in Mazzarella et al. (2020), the var_scaling, and len_scaling parameters which adjust the influence of BE matrix over the background field. Len_scaling controls the spatial decorrelation for the following five variables: unbalanced 195 velocity potential, unbalanced temperature, pseudo-relative humidity, unbalanced surface pressure and stream function. The use of a len_scaling factor of 0.5 reduces the variable perturbation length scale by 50%, ensuring that the water vapour increments are comparable with the weather radar range; therefore, this value has been adopted for the simulations.
The no-echo option, developed by Min and Kim (2016), allows the assimilation of null-echo within the radar observation range. This information reduces the excessive humidity and the contents of the following hydrometeors: snow, graupel, and 200 ha formattato: Tipo di carattere: Grassetto ha formattato: Tipo di carattere: Grassetto ha formattato: Tipo di carattere: Grassetto rainwater based on radar reflectivity, improving the convective precipitation predictability. In addition, a recent study (Lee et al., 2020) confirms the benefit of this option for the simulation of three summer convective events over the Seoul metropolit an area. For this reason, the no-echo option is used for this study.

Radar observation operator
The assimilation of CAPPI radar reflectivity for both 3D/4D-Var assimilation methods has been performed through the indirect 205 method (Wang et al. 2013a;Gao and Stensrud 2012). The new approach converts the observed reflectivity in the three hydrometeor mixing ratios in contrast to the direct method (Xiao et al. 2007) which only uses the rainwater mixing ratio.
The forward reflectivity operator considers the contribution of snow, rain and hail components and it is represented as: where Ze is the equivalent reflectivity factor, varies linearly between 0 at Tb = −5 °C and 1 at Tb = 5 °C, and Tb is the 210 background temperature from an NWP.
The rain contribution to the equivalent reflectivity (Smith et al. 1975) is calculated as follow: For the hail component the formulation based on (Lin et al. 1983;Gilmore et al. 2004) is adopted: 220 ( ℎ ) = 4.33 × 10 10 ( ℎ ) 1.75 Some microphysics schemes do not produce hail component, but WRFDA code recognises them and uses the qh variable as a graupel species qg (Lagasio et al., 2019).
Finally, the equivalent reflectivity (Ze), given by the sum of the three components Z(qs), Z(qr) and Z(qh), is converted in Z (dBZ) using: 225 = 10 log 10 The option, allowing the assimilation of in-cloud humidity from radar reflectivity, has been considered for this study (Wang et al., 2013a). The estimate of the saturated water vapor observation is produced considering the assumption that in cloud humidity is saturated. The in-cloud relative humidity assumed to be 100% where radar reflectivity is higher than a thresh old above cloud base, so that the estimated water vapor equals to the saturation water vapor that is calculated in Eq. (9) based on 230 the pressure and temperature of the background. In this paper the threshold is set to 25 dBZ (WRFDAdefault value). A full 10 description of the indirect assimilation method is provided in Wang et al. (2013a). In this case theThe forward observed operator H is defined by: where qv represents the specific humidity, rh the relative humidity, and qsh the saturated specific humidity of water vapor. It is 235 worth noticing that the assimilation of the in-cloud humidity is used in combination with the indirect assimilation. Thus, the numerical experiments also include the assimilation of the in-cloud humidity in addition to the hydrometeor species retrieved with the indirect method alone.Hence, this option allows the assimilation of vapor from in cloud region in addition to the three hydrometeor species calculated with the indirect method.

Model setup 240
The numerical simulations are carried out with the ARW-WRFV4.0 model (Skamarock et al., 2019). WRF is a mesoscale model, supported by National Center for Atmospheric Research (NCAR) and largely used by the atmospheric modelling community. The main features consist of a fully compressible Euler nonhydrostatic equations with hydrostatic option, staggered Arakawa C grid, and a highly portable and flexible code that optimize the use of computational resources. A two-way nesting configuration with two domains is used for this study: the mother domain with 379x431 grid points, covers 250 the Italian peninsula with a horizontal resolution of 3 km, while the inner domain (340x319 grid points) is centred over the Abruzzo region (central Italy) with a grid spacing of 1 km (Fig. 4a). To avoid compatibility issues with WRFDA, the vertical terrain following coordinates are used instead of terrain-following hybrid. Both domains adopt 40 vertical levels from the ground up 100 hPa. Because of the high spatial resolution, the convection is explicitly resolved. The same physical parameterizations used by the Center of Excellence in Telesensing of Environment and Model Prediction of Severe Events 255 (CETEMPS) meteorological-hydrological forecast system, are set for this work (Ferretti et al., 2014). More specifically, the microphysical processes are parameterized by using the WSM6 scheme (Hong and Lim, 2006), while the MYJ scheme (Janjic et al., 1994) is applied for the PBL. Shortwave and long wave radiation are considered through the Rapid Radiative Transfer Model (RRTM) scheme (Mlawer et al., 1997). Finally, the Noah land surface model is chosen to parametrize the land surface processes (Chen and Dudhia, 2001). The operational analysis/forecast cycle from the integrated forecast system (IFS) global 260 model of the European Centre for Medium-Range Weather Forecasts (ECMWF) with a spatial resolution of 0.1° x 0.1° are used for this work and the boundary conditions are update every 3 hours.

Design of experiments
A total of five experiments are carried out to evaluate the impact of 3D/4D-Var in cycling mode. All simulations started at 0000 UTC on 3 May 2018 and last for 24 hours. Both 3D-Var and 4D-Var are applied every hour in cycling mode from 0000 265 UTC to 0300 UTC assimilating the CAPPI reflectivity data at 2000, 3000 and 5000 m MSL. The same number of observations has been assimilated for both 3D-Var and 4D-Var simulations, considering a 10-minute assimilation window. More specifically the CAPPI are assimilated at 0000 UTC, 0010 UTC, 0100 UTC, 0110 UTC, 0200 UTC, 0210 UTC, 0300 UTC and at 0010 UTC, 0110 UTC, 0210 UTC and 0310 UTC (Fig. 4b, c, d). Moreover, two additional simulations are performed with the aim of evaluating the performance of cycling assimilation methods in warm start mode. For this purpose, a previously 270 numerical forecast, initialized six hours before, is used as background field. The same aforementioned CAPPI are assimilated in these experiments.
In table 1 a summary of all the experiments performed.

Verification methodologies
To the aim of evaluating the impact of both 3D-Var and 4D-Var in cycling mode, an objective comparison between the observed and forecast precipitation is performed. To this purpose, the rainfall data collected by DPC rain gauge network are (Jones,1999) is used to remap the rain data to the 3km domain grid. The statistical analysis is only performed in a restricte d area of the mother domain (Fig. 4a, red rectangle) that includes Lazio and Abruzzo regions (hereafter LA, Fig. 3), because these are the regions where relevant accumulated precipitations was recorded.
Finally, this work represents a preliminary study to investigate the reliability of the cycling assimilation before implementing it in an operational meteorological-hydrological chain like that operative at CETEMPS focused on the meteorological-285 hydrological forecast in the Abruzzo region.
To avoid the spatial limitations of using a point-by-point approach in the evaluation of quantitative precipitation forecast (Roberts, 2003) a filtering method is used. Based on the positive response in the recent literature Wang et al., 2013c;Romine et al., 2013;Gustafsson et al., 2018), the Fraction Skill Scores (FSS; Roberts and Lean, 2008), is computed for the precipitation assessment. A few accumulated periods are considered to the aim of investigating the forecast ability for 290 different kind of precipitation. For the 3-hourly accumulated precipitation a neighbourhood size of 3x3 grid points is used for this purpose.
In addition, the Receiver Operating Characteristic (ROC) statistical indicator is computed to strengthen the statistical anal ysis, evaluating how skilful are the cycling assimilation methods in terms of QPF. The ROC curve synthetizes the information obtained with different thresholds in one diagram, just comparing the hit rate against false alarm rate (Swets , 1973). Finally, 295 the area under the ROC curve (AUC) which is largely used to quantify the skill of a forecast system (Mason, 1982;Buizza and Palmer, 1998;Storer et al., 2019;Ferretti et al., 2020), is calculated for each simulation.
Before presenting the results, to the aim of clarifying the impact of reflectivity CAPPI at 2, 3 and 5 km in cycling 3D/4D -Var methods, the analysis increments (analysis minus first guess) are calculated at the end of the last assimilation cycle ( 031200 UTC). The assimilation of radar reflectivity mainly impacts on water vapour and cloud hydrometeors variables other than on 300 temperature and wind components. Therefore, the analysis increments for qr (rainwater mixing ratio, Fig. 5) and qv (water vapor mixing ratio, not shown), which best represent the added value of the radar reflectivity assimilation, are considered. For this purpose, the vertical level 15 (about 2000 m above ground) is selected because the influence of radar data is more relevant.
A qualitative comparison between the two assimilation methods points out that the 4D-Var is more impactful in terms of qr both with warm (Fig. 5a) and cold start (Fig. 5c) compared to 3D-Var (Figs. 5b and 5d). Furthermore, the larger analysis 305 increments of qr are along the Adriatic Sea near Abruzzo coastline in agreement with the assimilated CAPPI maps, which show high reflectivity values in this area. The distribution of the analysis increments for the qvvapour distribution(not shown here) is similar to the previous analysis increments (qr in Fig. 5), but the magnitude becomes larger than qr. The 4D-Var simulations produce greater qv increments compared to 3D-Var, in line with the next result.

Results 315
A preliminary qualitative comparison is performed with the aim of explaining the impact of cycling assimilation in short -term rainfall prediction. To this purpose, the differences between forecast and observed precipitation fields have been calculated for each simulation considering the rainfall amounts from 0900 UTC to 1200 UTC, when the FSS shows the highest values and the gap between the different simulations is more meaningful.
The results (Fig. 6) confirm the positive impact of cycling assimilation: both methods reduce the underestimation of the rainfall 320 (blue area) over the mountain area at the border between Lazio and Abruzzo regions. In this context, the 4D-Var and 3D-Var experiments with a warm start initialization show the best performances in this area (Figs. 6c and 6d), improving the precipitation forecast accuracy also compared to the CTL (Fig. 6e). Conversely, the two simulations in cold mode overestimate the rainfall along the coastal area of the Abruzzo region, even though they partially mitigate the error in the internal areas ( Figs. 6a and 6b). 325

Statistical Evaluation
To the aim of evaluating the performance of 3D-Var and 4D-Var assimilation methods in cycling mode, the FSS and ROC, 330 previously descripted in Section 5.2, have been calculated. The FSS is calculated in LA region, where significant precipitation occurred. To this purpose, several threshold values are used to evaluate the impact of cycling assimilation with light, moderate and heavy precipitation. Finally, the ROC and AUC are calculated to reinforce the statistical analysis and summarize the information obtained with the FSS.
To compare the 34D/43D-Var experiments in warm/cold start and their ability to reproduce the precipitation pattern, the 335 statistical index has been calculated considering the precipitation accumulated over three specified time periods: 1, 3 and 6 hours, respectively. The time series of FSS is presented in the sections below.

Hourly precipitation
The FSS for the 1 mm h -1 accumulated threshold (Fig. 7a) is calculated for all the experiments starting from 0600 UTC to 1500 UTC in LA region. The results are quite similar for all experiments in the first hour of simulation due to the small amount of 340 accumulated precipitation, even if CTL is starting from higher values than the other experiments. Later, all experiments display higher FSS values than CTL experiment (red dashed line), pointing out the positive feedback of cycling assimilation in the interval from 0700 UTC to 1200 UTC. Moreover, the simulations with warm initialization clearly show (Fig. 7a, blue and green lines) higher values demonstrating a larger impact on the precipitation forecast than those with cold start. On the other hand, CTL performs better in the last three hours of simulation, suggesting that the positive impact of radar reflectivity 345 assimilation is decreasing in time. The FSS computed for the 3 mm h -1 threshold (moderate precipitation) is shown in Fig. 7b.
The simulations with data assimilation confirm the improvements in term of QPF from 0700 UTC to 1200 UTC. In addition, the two experiments with 4D-Var (Fig. 7b, yellow and blue lines) performs better than 3D-Var, except in the first hours of simulation. This behaviour proves the positive impact of 4D-Var with moderate precipitation. However, similarly to the results for 1 mm h -1 threshold, the impact of radar reflectivity assimilation decreases toward the end of the simulation. Finally, the 350 FSS is calculated for 8 mm h -1 threshold (Fig. 7c) in conformity with World Meteorological Organisation (WMO-No. 544, 2003) that classifies the precipitation with a rain rate greater than 7.6 mm h-1 as heavy rain. According to the results for the 3 mm h -1 threshold, the FSS confirms the good performance of 4D-Var in cycling mode. More specifically the CYC4DAVAR_warm (Fig. 7c, blue line) shows the highest FSS values in the whole period with heavy precipitation, except for a very short period where the 4D-Var cold start is reaching higher values than the warm start. Conversely, the benefit of 355 3D-Var simulations is limited to the first hours of simulation.

360
The statistical analysis with hourly precipitation clearly shows the positive impact of assimilation in cycling mode with bot h methods in the initial hours of simulation. Moreover, the simulations in warm mode show the best performance compared to those in cold start, especially the CYC4DVAR_warm (Fig. 7c, blue line) when heavy precipitation occurred.

Three hourly precipitation
The statistical analysis is also carried out for the 3-hourly precipitation accumulation interval to deeper investigate the response 365 of cycling data assimilation for precipitation characterized by different mechanisms as for example convection or orographic uplifting. The FSS for 1 mm 3h -1 threshold (Fig. 8a) highlights the benefit of using 3D and 4D-Var with cold start for a short initial period (0600-0730 UTC), suggesting that the background field is probably more accurate than the previous forecast.
But the rate of increasing for both the warm start is larger than CTL and cold starts. Indeed, the FSS values for the warm start are higher than the cold starts up to 1200 UTC (Fig. 8a, blue and green lines). Later the improvement reduces, and CTL 370 performs better (Fig. 8a, red dashed line). The score is also computed for 3 mm 3h -1 threshold (Fig. 8b). Similarly to the results for 1 mm threshold the cycling 4D-Var in cold mode (CYC4DVAR_cold) improves the precipitation forecast at the initial time (0600 UTC to 0800 UTC), while the CYC4DVAR_warm mode confirms the highestr values from 0900 to 1200 UTC as well as the poor performance in the first hour. Moreover, the reduction in FSS values for both 4D-Var and 3D-Var compared to CTL in the final hours of simulation period proves that the influence of cycling is restricted to a short time ran ge. 375 Finally, the FSS is also calculated for threshold value of 10 mm 3h -1 in LA region (Fig. 8c). The highest values are associated with CYC4DVAR_warm (blue line), although all experiments show an improvement compared to CTL in the first hours of simulation period. Later, all the experiments converge to the CTL but the values for FSS decrease due to the small accumulated rainfall.
The statistical analysis performed using the three threshold values (1, 3 and 10 mm), proves the positive impact of cycling 380 assimilation with radar reflectivity in the interval 0600-1500 UTC and confirms the benefit of 4D-Var compared to 3D-Var.
Moreover, the two simulations with a warm start initialization show a low impact at 0600 UTC. The reason for this is probably the accuracy of the initial conditions which is lower than the ones used for the cold start simulation, despite the ECMWF operational analysis are used to simulate the first 6 hours for the warm start simulations. The reason for this is probably the poor dynamic balance of the initial fields at the beginning of simulations in warm mode. 385

Six hourly precipitation
The FSS is also calculated using the 6-hourly precipitation in LA region using three threshold values: 10 mm 6h -1 , 15 mm6h -1 and 25 mm 6h -1 . The FSS calculated for the 10 mm 6h -1 threshold (Fig. 9a) proves the positive impact of CYC4DVAR_warm (blue line) compared to the 3D-Var and CTL simulations for the whole simulation interval. The results also suggest that the cycling assimilation with warm start performs better than experiments in cold start for both assimilation methods. 395 According to the results for the previous analysisthreshold, the CYC4DVAR_warm confirms the best performance in term of FSS for the 15 mm 6h -1 threshold (Fig. 9b). Both CYC4DVAR_cold and CYC3DVAR_warm display a positive impact in QPF during the whole simulation period. Conversely, the CYC3DVAR_cold shows a worsening at 1800 UTC when compared to CTL and the other experiments.
The FSS is also calculated for 25 mm 6h -1 (Fig. 9c) threshold to investigate the impact of cycling assimilation with larger 405 accumulated precipitation. The CYC4DVAR_warm shows the highest FSS values, pointing out the benefit of 4D-Var in warm start mode in the whole simulation interval compared to CTL. Nevertheless, also the other simulations with cycling assimilation show a positive impact in the whole interval, even though with a slight improvement and up to 1800 UTC. In conclusion, the time series of FSS points out the benefit of using a cycling assimilation for radar reflectivity. The warm start with 3D/4D-Var assimilation methods confirm the improvement in terms of QPF, for: 410 • the hourly accumulated precipitation, highlighting the good performance in the localization and timing of the onset of the precipitation and for very light precipitation; • the 3-hourly accumulated precipitation highlighting the improvements for convective cells and orographic precipitation; • the 6-hourly accumulated precipitation.
On the other hand, the FSS depends on threshold values, thus a further statistical indicator has been calculated to the aim o f 415 endorse the previous results. Hence, the ROC curve, which summarizes the result obtained with several thresholds in one plot allowing for an easy comparison, is built for this study.

Receiver operating characteristic (ROC)
The Receiver Operating Characteristic (ROC), which compare the probability of detection (POD) and false alarm rate (FAR), is calculated to evaluate how the simulations are skilful in precipitation forecast. The 12-hourly precipitation accumulated 420 from 0600 UTC to 1800 UTC on May 3, 2018 for the NWP experiments and the rain gauge network are used to build the curve ( Fig. 10) for the LA region. To investigate the ability of cycling assimilation in predicting rainfall with light, medium and heavy intensity the following precipitation threshold are chosen: 1 mm, 12 mm, 24 mm, 36 mm, 48 mm, and 54 mm. The curves show low POD values for 1mm threshold, and a worsening compared to the CTL. On the other hand, the benefit of cycling assimilation is clearly found with higher threshold values. In fact, the steepness of 3D/4D-Var curves is greater than 425 CTL, suggesting a good forecast skill with moderate and heavy precipitation, namely from 12 mm to 48 mm. In this regard, the CYC4DVAR_warm (green line) shows the best performance, while the CYC3DVAR_cold (red line) reduces its impact with high threshold values. In addition, the area under the curve (AUC) is also computed to objectively compare each curve.

Conclusions
In this paper the impact of cycling 3D-Var and 4D-Var variational assimilation methods on the forecast of a heavy precipitation event, occurred in an orographically complex region namely the Central Italy, is evaluated. The reflectivity CAPPI (Constant Altitude Plan Position Indicator) obtained by Italian radar network at 2, 3 and 5 km are assimilated every hour in cycling mode 440 into the WRF model. The comparison among the experiments is performed by using a filtering approach, the Fraction Skill Score (FSS), for the Lazio-Abruzzo regions, where relevant rainfall occurred on May 3, 2018. In this regard, the statistical analysis is performed considering 1, 3 and 6 hourly accumulated precipitation with three different threshold values in order to evaluate the benefit of cycling assimilation with light, moderate and heavy precipitation. Finally, the ROC curve is built, to further evaluate the reliability of cycling assimilation in precipitation forecast. 445 The FSS time series for the hourly precipitation highlights the positive impact of radar data for both 3D-Var and 4D-Var assimilation methods compared to CTL. The benefit of using a cycling assimilation is clearly shown in the results for both light and moderate precipitation. However, the impact reduces in the last hours of the simulation when all experiments converge to the CTL. Conversely, the poor amount of precipitation at the start time, reduce the impact of both assimilation assimilating reflectivity data. In this regard, the cycling 4D-Var has a greater impact than 3D-Var experiments and consequently higher FSS value. More specifically, the cold start initializations for cycling both 4D-Var and 3D-Var show an improvement in terms of QPF compared to CTL at beginning of analysis, while the experiments in warm start perform better after few hours. This behaviour is probably caused by a slightly unbalanced initial field for the warm start simulations.
For what concerns the 6-hourly precipitation, the FSS for 10 mm 6h -1 threshold confirms the improvement of warm start 455 simulations compared to CTL and cold initialization. The CYC4DVAR_warm clearly displays the greatest FSS values at 15 and 25 mm6h -1 thresholds, pointing out the positive impact of radar reflectivity. Also, the 4D-Var in cold start and the 3D-Var with a warm initialization produce an improvement in QPF although it is smaller than CYC4DVAR_warm. On the other hand, the CYC3DVAR_cold shows a worsening in FSS.
Finally, the ability of cycling assimilation to reproduce the 12 h precipitation field is evaluated by using the ROC and the area 460 under the curve (AUC). The curves are calculated considering the period from 0600 UTC to 1800 UTC on May 3, 2018.
because of the significant rainfall. The comparison between the simulations confirms that cycling 4D -Var in both warm and cold mode is the best technique, indeed the highest value of AUC=0.91 is obtained. The CTL shows lower steepness than the cycling 4D-Var and an AUC of 0.88. Finally, the AUCs for the two simulations with 3D-Var, respectively 0.87 for warm initialization and 0.89 for the cold, are lower than 4D-Var and comparable with CTL simulation. Therefore, the impact of 3D-465 Var over 12h accumulated precipitation is less clear.
In conclusion, the cycling assimilation with 3D-VAR and 4D-Var methods for this heavy rain event, improves the reliability of the precipitation forecast, even if the positive impact reduces in time. Therefore, to further investigate the impact of cycling assimilation with 3D/4D methods and to generalize the achieved results, a larger number of events should be considered.
Moreover, the two simulations with a warm start initialization produce good results in terms of FSS, but the differences are 470 small compared to the cold simulations that perform better at the initial time. This behaviour suggests that the precipitation spinup time decreases in cycling assimilation with cold start, while for warm initialization this is not true. In addition, the cycling 4D-Var with warm start (CYC4DVAR_warm) shows better performance than 3D-Var over all precipitation accumulation intervals considered for this study. Finally, the ROC cures and the AUC values also confirms the benefit of 4D-Var in warm start. 475 The huge computational cost of 4D-Var was already highlighted in Mazzarella et al., (2020), in fact a simulation with 1h assimilation window needed more than 6 hours. As a result of this, we have developed the idea to apply the 4D-Var in cycling mode with an assimilation window of 10 minutes, the results of which are discussed above. For what concerns the computation time, we calculated the time needed to perform the three cycles of assimilation for both assimilation metho ds. Specifically, the 3D-Var takes approximately 15-20 minutes whereas the 4D-Var, more computationally expensive, required ~2h. On the other 480 hand, the use of 4D-Var with an assimilation window of 3 hours, takes over 12 hours. Thus, the cycling approach significatively reduces the computation time and allows for using 4D-Var in small weather centres too. All numerical experiments are performed on the ECMWF's Cray HPC using 1080 computational cores. The next step of this work will be to assimilate the radial velocity to improve the accuracy of wind field and vertical velocity, thus the positioning of convective cells. This opportunity allows us to complete the assessment of weather radar assimilation 485 in a 4D-Var cycling data assimilation. In addition, the impact of a wider data assimilation windows in cycling 4D -Var could be tested, in combination with a strategy with more outer loops. These solutions allow the assimilating of more data and to take into account the non-linear effects, producing significant increments in the analysis field. For what concerns the computational time, the three assimilation cycles with 4D-Var requires about 3 hours using the CRAY XC40 cluster at ECMWF. On the other hand, the cycling 3D-Var is significantly faster and takes less than 1 hour. However, despite the 4D-490 Var requires greater computational resources, it may be used for operational purpose when using the assimilation window adopted for this work. Lastly, the results of this study are helpful to decide which cycling assimilation methods will be implemented in the operational CETEMPS meteorological-hydrogeological chain and if a nowcasting algorithm based on cycling WRF 4D-Var may be applied..

Data availability 495
Rain gauges data are provided by DEWETRA data portal (http://www.mydewetra.org/), the platform is accessible upon request to the Italian Civil Protection Department (DPC). Radar data are obtained from the Italian Civil Protection Department (DPC), access to these data are restricted and readers must request it through contacting the lead author.