A meteo-hydrological modelling system for the reconstruction of river runoff : the case of the Ofanto river catchment

A meteo-hydrological modelling system has been designed for the reconstruction of long time series of rainfall and river runoff events. The modelling chain consists of the mesoscale meteorological model Weather Research and Forecasting WRF, the Land Surface Model NOAH-MP and the hydrology-hydraulics model WRF-Hydro. Two 3-months periods are reconstructed for winter 2011 and autumn 2013, containing heavy rainfall and river flooding events. Several sensitivity tests were performed along with an assessment of which tunable parameters, numerical choices and forcing data most impacted on the 5 modelling performance. The calibration of the experiments highlighted that the infiltration and aquifer coefficients should be considered as seasonally dependent. The WRF precipitation was validated by a comparison with rain gauges in the Ofanto basin. The WRF model was demonstrated to be sensitive to the initialization time and a spin-up of about 1.5 days was needed before the start of the major rainfall 10 events in order to improve the accuracy of the reconstruction. However, this was not sufficient and an optimal interpolation method was developed to correct the precipitation simulation. It is based on an Objective Analysis and a Least Square melding scheme, collectively named OA+LS. We demonstrated that the OA+LS method is a powerful tool to reduce the precipitation uncertainties and produce a lower error precipitation reconstruction that itself generates a better river discharge time series.The validation of the river streamflow showed promising statistical indices. 15 The final set-up of our meteo-hydrological modelling system was able to realistically reconstruct the local rainfall and the Ofanto hydrograph.


Introduction
The problem of reconstructing long time series of river runoff for climate impact studies and management purposes has been dealt with in different ways in the literature, from statistical and spatially lumped methods to the most recent deterministic and fully distributed reconstructions (Todini, 1995;Menabde and Sivapalan, 2001;Todini and Ciarapica, 2002;Rigon et al., 2006;Jones et al., 2008;Ruelland et al., 2008).
Here we show a deterministic modelling approach applied to a small river catchment in the Puglia region (Italy).The scientific community is making considerable efforts to increase the performance of high-resolution meteorological and hydrological models as well as to build up coupled modelling systems.However, modelling the spatial and temporal distribution of the water cycle, especially in case of extreme events, is both a research and operational challenge because the water cycle includes several processes which interact with each other and span a wide range of spatial and temporal scales.
G. Verri et al.: A meteo-hydrological modelling system for the reconstruction of river runoff Hydrology modelling has gained limited success in the past due to a lack of input data and to model shortcomings.The precipitation records are crucial for assimilation, initialisation and validation procedures, but they need to be collected with a fine spacing of stations.A rigorous and detailed representation of soil type and topography conditions is required as they locally modulate the synoptic-scale weather and affect the basin drainage.In addition, numerous and complex air-land and subsurface processes are involved in the local water cycle of river catchments, but they cannot be directly measured.
Relatively higher-resolution data on topography, land use, soil types and river routing have become available in the last decade.At the same time, the Quantitative Precipitation Forecast (QPF) has been significantly improved (Cuo et al., 2011).Ensemble approaches have been developed to embed the different sources of errors into a probabilistic forecasting (Buizza et al., 2008;Davolio et al., 2013).
Several open issues affect both the meteorological and hydrological modelling.QPF remains one of the most critical tasks for mesoscale meteorological models since the precipitation field is the end result of many multi-scales and interacting processes and is sensitive to topography, soil types and land use conditions (Krzysztofowicz, 1999;Hapuarachchi et al., 2011;Zappa et al., 2011).The grid spacing of mesoscale meteorological models does not allow us to fully resolve the scales of the single convective cells/systems (Moeng et al., 2007;Shin et al., 2013).Moreover, the quality of meteorological modelling is also critical for ensuring the quality of hydrological modelling as the uncertainties associated with the meteorological unknowns propagate into the hydrological models (Pappenberger et al., 2005;Zappa et al., 2010).
Finally, an additional source of uncertainty is due to the parameterisation of many of the physical processes involved in the water cycle, e.g.water infiltration through the soil column, groundwater drainage and the aquifer water storage (Krzysztofowicz, 2011;Gupta et al., 2005).These parameterisations imply many tunable coefficients which play an important role in modulating the soil water distribution.
Several hydrological models have been developed in the last few decades: from simple empirical models to shallow water (SW) systems with different levels of approximation, i.e. kinematic, diffusive or fully dynamics wave equations.
Two hydrological models that are currently extensively used for operational purposes are the Hydrologic Engineering Center's Hydrologic modelling System (HEC-HSM; US-ACE, 2015), developed by the US Army Corps of Engineers, and the TOPographic Kinematic APproximation and Integration (TOPKAPI; Todini and Ciarapica, 2002;Liu and Todini, 2002).The latter is currently employed for operational discharge forecasting by the Italian Civil Protection at ARPA-SIM (Agenzia Regionale Prevenzione e Ambiente -Servizio Idrometeorologico).HEC-HSM and TOPKAPI are both rainfall-runoff models and are based on the kinematic wave approximation of SW equations which solve the chan-nel streamflow.No additional representation of 2-D overland water flow is considered for neither the gravitational drainage nor the coupling with a land surface model.Nickovic et al. (2010) propose the HYdrology surface runoff PROgnostic Model (HYPROM), which to our knowledge is the only hydrological model that includes fully prognostic SW equations for representing both the 2-D overland water flow and the 1-D channel streamflow.In contrast, no representation of the subsurface physical processes is provided by HYPROM.
A more complex hydrology-hydraulics model is the fully distributed Weather Research and Forecasting Hydrological modelling (WRF-Hydro) system (Gochis et al., 2013), originally designed as the hydrological extension of the WRF atmospheric model (Skamarock et al., 2008).WRF-Hydro is based on the diffusive wave approximation for representing both the 2-D overland water flow and the 1-D river streamflow.In addition, the WRF-Hydro system fully solves the subsurface soil physics and is two-way coupled with the NOAH-MP land surface model (Niu et al., 2011).Due to all these valuable features we used this model in our study.
The aim of this study is to evaluate the ability of our integrated modelling system to simulate the water distribution in the catchment of the Ofanto river, which flows through Southern Italy.Our meteo-hydrological modelling system was set up for the first time in this basin.Two simulations were carried out over the winter 2011 and autumn 2013.Several rainfall events and dry periods characterise Southern Italy during the selected time ranges.The final goal of this work is a reliable reconstruction of the local rainfall and river runoff; thus the evaluation of the modelling performance was carried out, focusing on both precipitation and streamflow prediction.We describe both the precipitation reconstruction and the hydrograph results since mainly the former drive the quality of the hydrology of a river basin.Critical issues such as the precipitation modelling skill and the calibration of both NOAH-MP and WRF-Hydro models are discussed.A novel approach for the correction of the modelled precipitation is also presented.
The paper is organised as follows.Section 2 describes the study area.The modelling system and the experimental setup are presented in Sect.3. Section 4 discusses the modelling results.Conclusions are drawn in Sect. 5.

The study area
The basin of the Ofanto river, flowing through Southern Italy and ending in the Adriatic Sea, was chosen as case study in order to test the meteo-hydrological modelling chain we implemented over the central Mediterranean area (left picture in Fig. 1) with a focus on Southern Italy (right picture in Fig. 1).This is intended to be a relocatable case study as the final configuration of our meteo-hydrological modelling chain may be easily applied to investigate rainfall and runoff events in other study areas with similar physiographic characteristics.
The Ofanto river is a semi-perennial river, whose discharge is close to zero during the dry season but may significantly increase when heavy rain events occur and eventually cause the river to flood.The mean annual discharge at its outlet is around 15 m 3 s −1 ; minimum monthly climatology is 2.27 m 3 s −1 in August and reaches its monthly peak, 35 m 3 s −1 , in January.The local annual mean rainfall is about 720 mm with large precipitation gradients over small spatial and temporal scales; the annual mean temperature is around 14 • C (Romano et al., 2009).The watershed area (Fig. 2), covering Campania, Basilicata and Puglia in Southern Italy, is about 2790 km 2 , making it a medium-sized catchment (between 1000 and 10 000 km 2 ); the mean topography slope is 8 % and the total length is around 170 km, making it the second longest river in Southern Italy.The river source is located south of Torella dei Lombardi, a small village located at 715 m above the sea level.This is not the only source of the river; there are a few tributaries with a lower water volume which prevent the bed from becoming dry.The Ofanto basin consists of two distinct areas: the northeast and the southwest.The southwestern part, representing the upstream reach of the river, is mainly mountainous or hilly due to the Apennine range; the northeastern part, representing the downstream reach of the river, is a flat area which includes the flooding area of the river (see top panel of Fig. 2).The predominant soil type in the upstream subbasin is "loam" according to the United States Geological Survey (USGS) dataset, while a less pervious type, "clay loam", prevails in the downstream sub-basin (see Fig. 2c).The dot markers in the bottom panels in Fig. 2 indicate four monitoring points along the river network and by using an ArcGIS tool the Ofanto basin can be divided into four subbasins (Fig. 2b).In particular, the Calitri gauge separates the uppermost sub-basin, corresponding to a karst aquifer, from the three downstream sub-basins where a porous aquifer is located and favours the saltwater intrusion from the Adriatic Sea.
This case study is challenging in terms of both meteorological and hydrological modelling purposes.With regard to the hydrological modelling, this study deals with a mediumsized catchment with a "rain-runoff" response time which varies from several days up to a few hours, which makes the streamflow prediction highly demanding.Concerning the meteorological modelling, the case study is located in the Southern Italy, where several heavy rainfall and flash flood events have occurred in the last decades, triggered by lee cyclogenesis and convective instability (Federico et al., 2008(Federico et al., , 2009;;Moscatello et al., 2008;Miglietta et al., 2008;Mastrangelo et al., 2011).
In September 2000 a severe flood occurred in Soverato, a small town along the Ionian coast of Calabria (Bellecci et al., 2003).Another high-impact storm occurred over the Calabria region on 10-12 December 2003 (Federico et al., 2008).An intense weather storm hit Southern Italy and in particular the Puglia region on 12-14 November 2004 (Mastrangelo et al., 2011).A flash flood episode affected a small area in Puglia on 22 October 2005 (Miglietta et al., 2008).A tropical-like cyclone affected southeastern Italy on 26 September 2006 (Moscatello et al., 2008;Laviola et al., 2011).The small scales of motion involved, meso-β and meso-γ scales, make it difficult to numerically model these extreme events.
Two recent heavy rainfall events in 2011 and 2013 with consequent flooding of the Ofanto river are discussed in this paper.
3 The method 3.1 Experimental design of the meteo-hydrological modelling system A meteo-hydrological modelling system was implemented to reconstruct the Ofanto river runoff and is presented here for the first time.
The model chain consists of the meteorological mesoscale model WRF v3.6.1, the land surface model NOAH-MP and the hydrological model WRF-Hydro v2.NOAH-MP works as a sub-model of the WRF and WRF-Hydro system and is coupled in a two-way mode with both of them.The WRF and WRF-Hydro system are coupled one-way.It is worth highlighting that the one-way coupled meteorological and hydrological models have not been found to produce weaker performances than fully coupled systems.Senatore et al. (2015) recently proved that with the current state of development, the one-way or two-way couplings between WRF and WRF-Hydro show a comparable performance especially in terms of precipitation simulation.
The various modules of the modelling chain and how they face each other are shown in Fig. 3.
A detailed description of the equations and parameterisations which are relevant to the discussion of our results is provided in Appendix A.
The WRF model is a widely used mesoscale meteorological model which solves fully compressible, non-hydrostatic Euler equations.The model adopts a terrain-following hydrostatic pressure vertical coordinate.We chose 58 unevenly spaced levels and set the top of the model at 50 hPa.The two domains nested in two-way mode were considered as a coarse domain covering the central Mediterranean area with a horizontal resolution of 6 km and an inner domain over Southern Italy with a 2 km horizontal resolution (the domains are depicted in Fig. 1).
The two domain set-ups (Fig. 1) aim to capture the genesis and the development of the mesoscale cyclonic patterns responsible for the heavy rain events in the coarse domain; moreover, the finer grid mesh of the inner domain enables to reconstruct the local convection including the orography effects in the region of interest, i.e. the southeastern Italy.We tested different extensions and grid spacing of the coarse domain and we compared the two-domain approach with the one-domain-only set-up.We found that the two-way coupling mode improves the reconstruction of precipitation at local scales (not shown).
The analysis fields built by ECMWF-IFS (European Centre for Medium-Range Weather Forecasts -Integrated Forecasting System) with a 16 km horizontal resolution and 6 h frequency were adopted as WRF initial and boundary conditions.
A set of sensitivity tests (not shown) highlighted that the terrestrial datasets, i.e. the topography elevation and the land use categories, strongly affect the air-land fluxes and the near-surface atmospheric fields.Thus, the default USGS datasets with 800 m resolution for topography and land use have been upgraded with higher-resolution and more recent data: Corine 250 m land use categories and EU-DEM 30 m topography data both released by the European Environmental Agency (EEA); SRTM 90 m topography released by NOAA is adopted only over those regions (i.e.northern Africa) not covered by the EU-DEM dataset.
In addition, different numerical schemes for the parameterised atmospheric processes have been tested and compared by evaluating how they affect the simulation of the near-surface atmospheric fields.The final model configuration uses the RRTMG scheme (Iacono et al., 2008) for long-wave and short-wave radiation, the Monin-Obukhov scheme (Monin and Obukhov, 1954) for representing the surface sub-layer of the planetary boundary layer (PBL) and the Yonsei University scheme (YSU) is the non-local K-profile scheme (Hong and Lim, 2006) that represents the PBL mixed sub-layer.The microphysics was based on the Thompson double-moment six-class scheme (Thompson et al., 2008) for both domains.The cumulus-convection parameterization was based on the Kain-Fritsch scheme (Kain and Fritsch, 1993) in the coarse domain, while no convection scheme, meaning that the convection is assumed to have been solved explicitly, was found to perform better in the inner domain.Our sensitivity tests shows that the explicit convection works better than the convection parameterization in the inner domain, as its grid spacing is in the "convection-permitting" scale range (Prein et al., 2015).This is an expected result, largely documented by previous studies on severe convective weather forecasts: Done et al. (2004), Weisman et al. (2008), Kain et al. (2008) and Schwartz et al. (2009Schwartz et al. ( , 2010)), among the others.
Table 1 summarises the terrestrial datasets and parameterisation schemes that we adopted as a result of the sensitivity tests.
Overall our experimental design is based on the past studies of WRF for local rainfall events in the same region that stressed the two-way nesting: Miglietta et al. (2008), Moscatello et al. (2008), Federico et al. (2008), Laviola et al. (2011) and Mastrangelo et al. (2011), among the others.
The NOAH-MP land surface sub-model was used to solve the partition of the surface water into an infiltration rate and surface runoff together with the water content and temperature of four soil layers up to 1 m below the ground level.The surface runoff is represented as the excess of surface infiltration capacity, while subsurface runoff is the gravitational drainage at a 1 m soil depth, i.e. through the bottom of the solved soil column.The snow modelling is also active in Table 1.Terrestrial datasets and parameterization settings adopted over WRF coarse domain 1 (6 km grid spacing) and inner domain 2 (2 km grid spacing).NOAH-MP model: a multilayer snow pack, the snow albedo and the melting/refreezing capability are solved by NOAH-MP.Moreover, the evaporation component coming from the snow sublimation is added and the evaporation component coming from the canopy water is split into the rainfall and the snowfall terms.The ECMWF analyses used for computing the initial and boundary conditions provide also the accumulated snow depth at the ground level.For our case studies, the snowfall and the melting processes do not seem to play a crucial role.

WRF-ARW set-up
The added value of the whole WRF-Hydro system with respect to the NOAH-MP is the ability to laterally route both the surface and subsurface water flow, as well as representing their interaction.The surface runoff is routed by a 2-D SW system (Eqs.A6-A8) and is also refined on the catchment area: the channel network gains water inflow from both the surface runoff and the aquifer discharge, and the channel streamflow is solved by 1-D SW system (Eqs.A11-A12).
We set up the WRF-Hydro model over the WRF inner domain.A detailed catchment routing grid is computed using a GIS procedure starting with EU-DEM 30 m topography data.The catchment grid is reproduced with a high level of accuracy (Fig. 2d): the drainage directions are first drawn and the river network is then refined by identifying all the branches and the hierarchy of tributaries using Strahler's method (1952).We chose a grid spacing equal to 200 m, which is 10 times higher than the WRF and NOAH-MP spatial grid.The aquifer water storage is switched on and the aquifer grid is assumed to identically match the watershed grid; four sub-basins (Fig. 2b) are defined as the areas located upstream of the monitoring points set along the river network, which enabled us to customise the coefficients of the aquifer recharge/discharge over smaller areas.
Two simulations were performed over January-March 2011 (hereafter "Experiment 1") and November-December 2013 (hereafter "Experiment 2").The selected time windows included several rainfall events of different intensities.The strongest weather storms occurred on 1 March 2011 and 1 December 2013, hereafter referred to as "Event 1" and "Event 2", and were followed by the flooding of the Ofanto river.Various features of the experiments are summarised in Table 2.
Figure 4 shows the concatenation procedure we adopted for both experiments: a chain of 72 h long simulations was carried out and the reinitialisation option was chosen for WRF, while a restart option was adopted for WRF-Hydro.We chose a frequent reinitialisation strategy following previous studies carried out with regional climate models on seasonal and sub-seasonal scales (Qian et al., 2003;Koster et al., 2010;Lucas-Picher et al., 2013).These studies highlighted the benefits of working with a concatenation of short simulations rather than a standard continuous simulation.Above all, the reinitialisation mitigates the problems of systematic errors and numerical drift and thus improves the accuracy in reproducing the local scale precipitation.The hydraulics component of WRF-Hydro system is initialised with the NOAH MP overland and subsurface water flows that are dry at the initial time.Thus a spin-up period is required to laterally route the groundwater of the basin and to allow the river network to reach a steady state.Senatore et al. (2015) considered monthly spin-up for evaluating the WRF-Hydro results and we decided to follow the same strategy.

A two-step correction of the precipitation: objective analysis (OA) mapping and observation-model merging
The simulation of the localisation, amount and timing of precipitation is crucial for the reconstruction of a river runoff time series but uncertainties are large in mesoscale models, particularly due to unresolved meso-β and meso-γ scale processes.In our experiments the horizontal resolution of the WRF inner domain was 2 km.This is quite a high resolution for mesoscale modelling and considers the convection as explicitly resolved by the model.Sensitivity tests done at an early stage confirmed that the "explicit convection" worked better than any convection scheme.However, a resolution of a few kilometres is a "grey zone" for representing the con-vection, because at these scales the power spectrum of the turbulence reaches its peak (Moeng et al., 2007;Shin and Hong, 2013).This means that the WRF model does not fully reproduce the convective motions and consequently neither reproduces the rainfall events with very local scale features.
In order to increase the performance of the precipitation reconstruction by WRF, we developed a two-step correction algorithm.First we used an OA technique to address the statistical interpolation of the scattered precipitation observations.Mathematical details on the OA technique and the calibration of the OA parameters are provided in Appendix B. Secondly, a least square (LS) melding algorithm was developed in order to merge the OA optimal estimate with the modelled precipitation.The assimilated precipitation is thus given by the following correction formula: where P b is the modelled precipitation, P o is the OA optimal estimate and P a is the corrected precipitation.In addition, σ 2 o /σ 2 f is the normalised variance associated to the OA (formula B4 in Appendix B) with values in the [0,1] range, as shown by the bottom panel of Fig. B1, and σ 2 b /σ 2 f is the normalised variance associated with the modelled precipitation also with values in the [0,1] range.Inside the Ofanto basin, the model normalised variance was defined for each grid point as the standard deviation of the modelled precipitation divided by its maximum value, Hereafter, our two-step correction procedure based on OA mapping and LS formula will be referred to as the "OA+LS" method.
Overall, we found that the OA+LS method is robust for the correction of the precipitation field.Panels in Fig. 5 show the increased performance of the corrected precipitation with respect to the in situ rain-gauge data for 2 days inside Experiment 1 and Experiment 2. The left panels in Fig. 5 show the simulated precipitation field, while the right panels provide the corrected precipitation field.The daily observed precipitation is shown by the stations network.There are 27 raingauge stations managed by the Civil Protection of Puglia Region, which are regularly distributed over the whole catchment as shown in Fig. 2a.The records cover the whole simulation periods with a 30 min frequency.It should be noted that a quality control of the observations was performed and two stations (named Cerignola and Borgo Libertà) were removed from the validation of Experiment 2 as they differed too much from the surrounding stations.
The added value of the OA+LS method is further discussed and quantified in the next section.

Results and discussion
The evaluation of the performance of our modelling system focuses on two fields: the precipitation and the river streamflow.We will show the calibration of the tunable coefficients involved in the parameterisation schemes and discuss the validation of the modelling results.
Critical issues such as uncertainties in the precipitation and the required spin-up of the meteorological simulations are stressed.
This study highlights that an integrated modelling system, which includes both the surface and subsurface runoff along with the aquifer water storage, is required to obtain a reliable reconstruction of heavy rainfall and flooding in the Ofanto catchment.

Mesoscale hydro-meteorological features
Figures 6 and 7 provide the mesoscale maps of the two severe weather events occurred on 1 March 2011 (Event 1) and 1 December 2013 (Event 2).
The 500 hPa geopotential maps highlight how the upperlevel features affect the lower-level cyclogenesis.WRF maps for Event 1 show a strong trough of low pressure at 500 hPa centred over the western Mediterranean Sea (left panel in Fig. 6), which is due to a cold front (not shown) progressing eastward.At lower levels a strong synoptic wind, coming from the southeast and blowing over the warm Ionian Sea, reaches the Italian Peninsula (right panel in Fig. 6); a weak cyclonic pattern is centred over the Tyrrhenian Sea with an associated sea level pressure gradient of 10 hPa moving in the middle of the core (blue patch in Fig. 6a). Figure 7a shows the 500 hPa geopotential maps for Event 2: a weak trough covers the western Mediterranean Sea and the Atlas Mountains in the upper troposphere, with a small but deep core south of Sicily.This corresponds to a strong cyclonic circulation at a lower level (right panel Fig. 7) with a sea level pressure gradient reaching 16 hPa in the cyclone's eye.This cyclone is situated almost directly beneath the cutoff low in the 500 hPa height field and corresponds to a southerly winds carrying warm, moist air reaching Southern Italy and a colder wind developing downslope of the Balkans.
These cyclones triggered Event 1 and Event 2 over Southern Italy with heavy local rainfall and river flooding.As detailed in the Puglia Civil Protection reports, anomalous rainfall hit the Puglia region during both events.A precipitation peak of 186.9 mm day −1 was recorded on 1 March 2011 (Event 1) at the Quasano station in the middle of Puglia, exceeding its historical maximum value of 116 mm day −1

The sensitivity to the initialization time
The simulation of the precipitation field may suffer from a initialisation time that is particularly close or far from the occurrence of precipitation peak events.In the first case the model is unable to develop the mesoscale features required to trigger the local weather pattern, and in the second case the numerical drift may affect the simulation results (Fiori et al., 2014).We avoided the second case by using a concatenation procedure which consists of a chain of WRF 72 h long and reinitialised simulations.
In addition to Experiment 1 and Experiment 2 we performed extra WRF 72 h runs focusing on specific events to test the sensitivity of the simulated precipitation in relation to the initialization time: the panels of Fig. 8 highlight the differences between the 24 h cumulated precipitation on 18 February 2011 started 14 and 38 h before the rain peak of Event 1.The left panel shows the 24 h cumulate precipitation modelled by WRF with initial time set 14 h before the rain peak as it results from the chain of 72 h long simulations of Experiment 1, the right panel shows the 24 h WRF precipitation by choosing 38 h as the lead time.Overall, we found that the WRF ability to correctly reproduce rainfall events is lowered by a initialization time that is too close to the peak events.
We conclude that our WRF model would need to be reinitialised approximately 1.5 days earlier than the start of the heavy rain events to increase skill in the prediction of precipitation.For this reason as a future step we plan to develop a robust WRF ensemble, which consists of overlapping chains of 72 h simulations with a delayed start time.

Validation of the precipitation
For a comprehensive assessment of the QPF performance, we calculated the average bias, correlation coefficient (CORR) and coefficient of RMSE variation (CV(RMSE)) across all stations.it should be noted that the CV(RMSE) is computed as the root mean square difference between modelled and observed values, i.e.RMSE, divided by the model standard deviation.We believe that the CV(RMSE) indicator gives more rigorous information on the accuracy of the numerical results than the RMSE and the NRMSE (i.e. the RMSE divided by the mean observed value).This is because it weights the www.nat-hazards-earth-syst-sci.net/17/1741/2017/Nat.Hazards Earth Syst.Sci., 17, 1741-1761, 2017 model-observation scatter with respect to the variance of the model time series.This also makes the comparison between experiments performed over different time ranges with different extreme events more meaningful.Table 3 summarises the statistical indices: we consider the 24 h cumulated precipitation at the model grid points nearest to the basin gauge stations (top panel of Fig. 2).Experiment 2 shows a higher bias (computed as the modelled minus the observed value) than Experiment 1, as well as a better correlation and a lower or equal CV(RMSE).This is expected since Experiment 2 is characterised by an initial period (i.e.November 2013) of continuous rainfall and a second, almost dry, period, while a series of shorter rain events took place during Experiment 1, making the simulation of single events difficult.
Similarly to the studies by Yucel and Onen (2014) and Senatore et al. (2015) based on the WRF model, we found that our model set-up tends to overestimate the local rainfall (positive bias).The correction we propose based on OA+LS method strongly reduces this tendency and generates a weak underestimation.The correction procedure reduces CV(RMSE) by 84 % in Experiment 1 and by 64 % in Experiment 2 and increases CORR by 41 % in Experiment 1 and by 14 % in Experiment 2. Thus, our correction procedure proves to have a great impact on the precipitation modelling performance.The statistical indices we found are even better than those obtained by more sophisticated methods in similar studies with the same model.For example, the 3DVAR assimilation scheme by Yucel and Onen (2014) reduces the precipitation RMSE by 3.3 % with respect to the no correc- Overall, we conclude that in our experiments modelled and corrected WRF precipitation show a good agreement with the gauge stations in the Ofanto catchment in comparison with the results of similar studies (e.g.Fiori et al., 2014;Yucel and Onen, 2014;Senatore et al., 2015;Givati et al., 2016).The correlation of WRF precipitation reached by Yucel and Onen (2014) was 0.364 with a 3DVAR assimilation scheme.A relatively high correlation of WRF precipitation, i.e. 0.92, was found by Givati et al. (2016) by increasing the grid spacing up to 3 km.This result refers to only the WRF cumulated precipitation between +6 and +30 h with respect to the start time, and each simulation is 30 h long with the reinitialisation option.Our model set-up, in contrast, considers 72 h as the simulation range plus the reinitialisation option, and the validation is performed over the whole simulation range.

Calibration of NOAH-MP and WRF-Hydro tunable parameters
A calibration procedure of the tunable coefficients of both NOAH-MP and WRF-Hydro models was carried out in order to realistically reproduce the Ofanto hydrograph.As a first step we adopted an automated calibration procedure, based on the PEST software (Doherty, 2002).This procedure minimizes an objective function, given by the sum of the mean squared differences between the modelled and observed river streamflow, using the Gauss-Marquardt-Levenberg non-linear LS method.Several tests were carried out and we identified the most relevant parameters to be calibrated in our specific case study.The coefficients with a high correlation (i.e.|corr| > 0.9) or the ones that preserved almost the original values after the PEST tests have been excluded.Thus we reduced the original set of 25 tunable parameters to the 7 that are found to play a key role in the Ofanto basin.They are the surface roughness scaling factor, which controls the hydrograph shape and the timing of the peaks; the infiltration coefficient; the saturated hydraulic conductivity; and the aquifer coefficients, which control the total water volume.
As a second step, we carried out a manual calibration.Although this is a fairly rough approach, it avoids the uncertain-Table 4. Tuned coefficients of WRF-Hydro/NOAH-MP for both experiments.OVROUGHRTFAC is the overland roughness scaling factor, REFKDT is the infiltration coefficient and REFDK is the saturation of soil hydraulic conductivity.
ties arising from the tuning of highly correlated parameters when a non-linear LS method such as PEST is adopted.
The optimal values of the WRF-Hydro tuned parameters are listed in Tables 4 and 5.
Among the WRF-Hydro parameters controlling the hydrograph shape, Manning's 2-D and 1-D roughness coefficients play a crucial role as they are involved in the empirical formula used to compute the discharges of both the 2-D overland water flow (Eqs.A9-A10) and the 1-D channel streamflow (Eq.A14).Manning's 2-D roughness coefficients are indexed using the land use categories, while 1-D roughness coefficients are assigned on the basis of Strahler's stream order.In order to upgrade the computation of the 2-D roughness coefficients, we replaced the default USGS land use categories with the higher-resolution and updated Corine data (EEA dataset).We refined also the computation of Strahler's order and thus the 1-D roughness coefficients by adopting the higher-resolution and updated EU-DEM topography data (EEA dataset).
The 2-D roughness coefficients can also be calibrated using a scaling factor whose values vary between 0 and 1.0, where values equal to 1.0 mean that Manning's 2-D roughness coefficients are not changed.We found this factor needs to be adjusted and we progressively reduced it in order to obtain faster water streaming: 0.05 is the optimal value to best capture the timing of the peaks compared with the observed hydrograph.
The NOAH-MP sensitive coefficients are the infiltration and the saturated hydraulic conductivity.They affect both the surface water budget and the moisture content of the NOAH-MP soil layers through the parameterisation of the infiltration capacity (Eq.A4).They also indirectly condition the WRF-Hydro overland water flow through the source term of the Table 5. Tuned coefficients of WRF-Hydro/aquifer law for Experiment 1 and 2. Z ini is the initial value of the aquifer water depth, Z max is the maximum value of the aquifer water depth, α is the exponential law coefficient, C is the volume capacity of the aquifer.Different optimal values are set depending on sub-basin and season.
Additional sensitivity tests pointed out the seasonality of the soil physics.We found that the coefficients for the infiltration and the saturated hydraulic conductivity are seasonally dependent and thus different values are assumed in the two experiments.In winter, the soil is expected to be wetter than in autumn and the soil porosity lower, which implies that the infiltration coefficient value is fixed as lower and saturated hydraulic conductivity as higher in Experiment 1 than in Experiment 2.
The calibration of WRF-Hydro aquifer was also performed manually.The aquifer discharge directly feeds the river streamflow by affecting the river baseflow.Thus, the manual calibration of the aquifer coefficients was carried out by comparing the simulated and observed hydrograph.The tuned values of the aquifer coefficients are listed in Table 5.
Our sensitivity experiments highlighted that the aquifer discharge is dependent on soil type and season.The clay loam soil type of the Ofanto low valley is much less pervious than the upstream loam soil type.Thus, the low valley (sub-basins 2, 3 and 4 of our catchment) is characterised by a lower hydraulic conductivity, which tends to counteract the upward aquifer discharge.We set the values of the initial water depth of the aquifer, the exponential law coefficient and the volume capacity lower in the downstream sub-basins 2, 3 and 4 than in the upstream sub-basin 1, while the maximum depth of the aquifer is higher.
Experiment 1 also shows higher volume capacity coefficient and lower maximum depth than Experiment 2, as the soil column is expected to be wetter and thus the elevation of the water table is higher in winter than in autumn.
This study found that the soil infiltration and the aquifer water storage parameterisations should be seasonally dependent.This means that the present parameterisations of these processes are not capable to capture the complexity of the groundwater physical processes.

Validation of the runoff
The observed water level at the Cafiero station was used to validate the river runoff modelled by WRF-Hydro.To compare the model versus observed river level, we used the "stage-discharge" relationship for the Ofanto river at the Cafiero station, which converts the model runoff into the water level.This is because the water level as computed by WRF-Hydro is conditioned by the channel geometrical parameters such as side slope, bottom width and channel roughnesses.These parameters are prescribed as a function of the Strahler's stream order instead of being customised for the specific catchment, which makes the model water level unreliable.
Figure 9 shows the observed and modelled hydrograph of the Ofanto river in Experiment 1 by using the simulated precipitation (Fig. 9a) or the corrected precipitation (Fig. 9b).Similarly Fig. 10 refers to Experiment 2. The gap in the ob- In both experiments, working with the model precipitation, the river water level tends to be overestimated (Figs.9a  and 10a).This is reduced by the corrected precipitation (Eq. 1) and thus the runoff peaks are better captured (Figs.9b  and 10b).
For an overall assessment of the streamflow reconstruction, we calculated the CV(RMSE) and CORR at the Cafiero station.Table 6 summarises the statistical indices, considering the model grid point nearest to the Cafiero gauge station and an hourly frequency.The first significant result is that the OA+LS method is a powerful way of dealing with the precipitation uncertainties, and it has a positive impact on the hydrological reconstructions.The OA+LS procedure reduces the WRF-Hydro CV(RMSE) by 20 % in Experiment 1 and by 6 % in Experiment 2 and increases CORR by 24 % in Experiment 1 and by 19 % in Experiment 2. This is significant also in comparison with similar studies based on WRF-Hydro.For example, Yucel (2015) shows that the assimilation of the precipitation by a 3DVAR method reduces the hydrograph RMSE by 7.6 % but does not affect the correlation, which is equal to 0.90 and 0.89 without and with the assimilation of the precipitation field, respectively.
It is important to stress that by validating the Ofanto hydrograph we were able to calibrate the OA tunable coeffi-cients: the optimal values are those that ensure that the simulated hydrograph is the closest to the observed one.
Some shortcomings are evident in the representation of the Ofanto hydrograph.The first problem is the timing of the runoff peak after a severe rain event.If the observed water level peak occurs at times close (less than 24 h) to the rain peak event, the simulated water level has a low skill.A short lead time (less than 24 h) is the reason for the delay in the simulated peak centred on 23 November 2013 and for the peak underestimation starting on 2 December 2013 (Fig. 10), both of which are associated with local-scale rainfall.In contrast, the peaks observed on 19 February 2011 and on 6 March 2011 are well captured despite the rain short lead time because they are linked to weather events with large spatial scales.
It should be also noted that the Event 2 onset overlaps the start time of WRF 72 h simulation (Table 2) and this probably affects the underestimation of the runoff peak starting on 2 December 2013.
For a comprehensive analysis of the shortcomings in the representation of the Ofanto hydrograph, the limited coverage of the rain-gauge stations and the quality of the observed values should be considered.The overestimation of the hydrograph we found in November 2013 (Fig. 10), when WRF-Hydro is forced by the corrected precipitation, enabled us to identify and remove some "outliers" among the observed precipitation values on 5 and 11 November.However, the overestimation of the hydrograph persists due to the WRF overestimate of precipitation on 4 and 10 November.The daily maps of modelled and assimilated precipitation in Fig. 11 show that the WRF overestimate cannot be removed by the OA+LS corrections as there are no rain gauges in the uppermost region of the Ofanto basin where the rainfall peaks occurred.
Despite the predictability gaps in the precipitation field and the limited coverage of the rain gauges, the final configuration of our meteo-hydrological modelling chain with an appropriate calibration was found to reasonably reconstruct the Ofanto hydrograph and to correctly reproduce the water level peaks as well as the plateaus.
To conclude, we compare the water level reconstruction at Cafiero station with and without the aquifer and the same but without the hydro component of the modelling system.The blue time series in Fig. 9 shows the water level when the aquifer in not activated with respect to the red time series with the aquifer switched on.The comparison shows that the river baseflow is affected by the aquifer parameterisations and that a better representation is obtained with the aquifer.However, the aquifer parameterisations do not impact the quality of the reconstruction because of the small Ofanto catchment aquifer capacity, as shown in Fig. 9 (i.e.CV(RMSE) index reduces of only 2 % when the aquifer is switched on and the correlation is almost the same).
Figure 12 shows that coupling the land surface model NOAH-MP with the WRF-Hydro model component is cru-  cial.The "column-only" land surface model NOAH-MP parameterises the surface runoff through Eq. (A2), which is inadequate to represent the Ofanto hydrograph as shown by the blue time series in Fig. 12.

Summary, conclusions and future plans
This study investigated the ability of a new meteohydrological modelling system, based on WRF and WRF-Hydro models, to reconstruct the local water cycle of a small catchment in Southern Italy.We chose a challenging case study: the semi-perennial Ofanto river with a small-sized catchment and a porous aquifer in the downstream region.The river basin is also located in the Southern Italy, which is frequently subjected to flash flood events.We chose two time windows characterized by the occurrence of severe weather events with flooding of the Ofanto.
Our study provides the first implementation of the WRF-Hydro system in this region.
The aim was to highlight the strengths and weaknesses of the final set-up of our meteo-hydrological modelling system, as well as to develop a useful tool for reconstructing the water runoff in rivers for several months.
One of the novelties of this study lies in our OA+LS method, which we used to correct the modelled precipitation.The OA statistically interpolates the scattered observations on the WRF model regular grid and the LS method is used to merge the OA optimal estimates with the modelled precipitation.This is a powerful method to deal with precipitation uncertainties, providing a positive impact on the hydrological reconstructions.The OA+LS procedure improved the precipitation estimate by reducing RMSE by 84 % and increasing correlation by 41 %; water level connected to runoff was improved with a RMSE reduction of 20 % and correlation increase of 24 %.
The quality of the modelling system was proven by the validation of both precipitation and water level predictions.We obtained promising statistical indices for both fields, also in comparison with recent studies dealing with WRF and WRF- Hydro models (e.g.Fiori et al., 2014;Yucel and Onen, 2014;Yucel et al., 2015;Senatore et al., 2015;Givati et al., 2016).
The calibration of NOAH-MP and WRF-Hydro tunable coefficients highlighted that the infiltration and the aquifer coefficients are seasonally dependent.This means we need to account for increased complexity in the parameterisation of the groundwater physical processes.
The performance of our modelling system was affected by various error sources: (i) the predictability limit of the precipitation field using a mesoscale meteorological model due to the meso-β and meso-γ scales involved in the rainfall events; (ii) the sensitivity of the precipitation predictions in relation to the initialisation time, which cannot be too close or far from the rainfall events.The WRF model set-up thus needs a spin-up period of about 1.5 days before the start of the rain event in order to be able to realistically reproduce the local weather pattern, and 72 h was found to be an adequate reinitialisation range.
Our next step is to exploit the OA+LS method as a flexible technique to correct other WRF variables.This technique could be embedded into an operational meteo-hydrological forecasting system.We are also planning to develop a more robust WRF ensemble which consists of overlapping chains of 72 h simulations with a delayed start time.
Overall, we highlighted the two-way feedback existing between a proper reconstruction of the meteorological events and the hydrological ones.A reliable description of the river hydrograph goes through a proper description of the meteorological and soil processes, with the precipitation field playing the most relevant role.At the same time the validation of the river hydrograph works as a effective post-processing tool to calibrate the water infiltration through the soil column and the aquifer recharge/discharge and to correct the modelled precipitation with the OA+LS method.
More research is required to establish a better groundwater modelling that at the moment considers seasonally dependent, ad hoc values of the soil infiltration and the aquifer water storage.We plan to evaluate different parameterisations of the aquifer recharge/discharge.Overall a reduction of the parameterisations involved in the WRF-Hydro system could be desirable.
Data availability.Our experiments have not been conceived for collecting a data repository, and thus our output data is available on request but not directly downloadable.aquifer water z = z + Q bot dt.Tunable parameters are the initial value of the aquifer water depth z ini (units of millimetres), the maximum value of the aquifer water depth z max (units of millimetres), the exponential law coefficient α and the volume capacity of the aquifer C (units of m 3 s −1 ).

Figure 2 .
Figure 2. The Ofanto river Catchment.(a) Topography height (units of m) and location of 27 rain-gauge stations.(b) The whole basin and the 4 sub-basins (coloured zones) defined as the areas upstream of the selected monitoring points (black dots).(c) USGS Soil Type Categories in the region of the Ofanto basin.(d) The flow accumulation grid defined by the number of grid cells which drain into an individual cell along the river network grid.

Figure 4 .
Figure 4.The concatenation procedure of the simulations.

Figure 5 .
Figure 5. Maps of 24 h cumulated precipitations (in mm day −1 , colours) during the peak events on 1 March 2011 (top panels) and 1 December 2013 (bottom panels).Shaded maps of modelled (left panel) and assimilated (right panel) precipitation with overlapped observed spots over the Ofanto basin.

Figure 9 .
Figure 9. Validation of the Ofanto discharge for Experiment 1 at Cafiero station.(a) Modelled precipitation.(b) Assimilated precipitation.The blue time series refers to the additional experiment performed with the aquifer switched off.

Figure 10 .
Figure 10.Validation of the Ofanto discharge for Experiment 2 at Cafiero station.(a) Modelled precipitation.(b) Assimilated precipitation.The gaps in the black line are due to the river flood which damaged the gauge station.

Figure 11 .
Figure 11.Maps of 24 h cumulated precipitations (in mm day −1 , colours) on 4 November 2011 (top panels) and 10 November 2013 (bottom panels): shaded maps of modelled (left panel) and assimilated (right panel) precipitation with overlapping observed spots on the Ofanto basin.

Figure 12 .
Figure 12.Comparison of the Ofanto discharge for Experiment 1 at Cafiero station as provided by the best WRF-Hydro set-up and by NOAH-MP.

Figure B1 .
Figure B1.Mapping of daily precipitation (mm day −1 ) on 1 March 2011 as carried out by objective analysis (OA).(a) OA optimal estimate with overlapping observations; (b) OA error.

Table 2 .
Details on the experiments.

Table 3 .
Statistical indices for validation of modelled and assimilated precipitation by comparison with rain-gauge stations in the Ofanto basin.

Table 6 .
Statistical indices for validation of river streamflow by comparison with Cafiero gauge station.