Implementation of WRF-Hydro at two drainage basins in the region of Attica, Greece

An integrated modeling approach for simulating flood events is presented in the current study. An advanced flood forecasting model, which is based on the coupling of hydrological and atmospheric components, was used for a twofold objective: first to investigate the potential of a coupled hydrometeorological model to be used for flood forecasting at two drainage basins in the area of Attica (Greece) and second to investigate the influence of the use of the coupled hydrometeorological model on the improvement of the precipitation forecast skill. For this reason, we used precipitation and hydrometric in-situ data for 7 events at two selected drainage regions of Attica. The simulations were carried out with WRF-Hydro model, which is an enhanced version of the Weather Research and Forecasting (WRF) model complemented with the feedback of terrestrial hydrology on the atmosphere, where surface and subsurface runoff were computed at a fine resolution grid of 95 m. Results showed that WRF-Hydro is capable to produce the observed discharge after the adequate calibration method at the studied basins. Besides, the WRF-Hydro has the tendency to slightly improve the simulated precipitation in comparison to the simulated precipitation produced the atmospheric only version of the model. These outcomes provide confidence that the model configuration is robust and, thus, can be used for flood research and operational forecasting purposes in the area of Attica.


Advanced Research WRF
The Advanced Research Weather Research and Forecasting model Version 3.9.1.1 was used in this study (Skamarock et al., 2008) for the atmospheric only simulations. Numerical simulations were carried out with four nested grids ( The parametrization schemes used for the simulations are given in Table 2. For the cloud microphysics processes, the WRF Single-Moment 6-Class Microphysics scheme (WSM6) is used for the coarser domain (d01; Hong and Lim, 2006) and WRF Double-Moment 6-Class Microphysics scheme (WSM6) is used for the higher resolutions domain (d02, d03 and d04;Hong et al, 2010). The shortwave and longwave radiation fluxes were parameterized with the Dudhia (Dudhia, 1989) and the Rapid Radiative Transfer Model (RRTM;Mlawer et al., 1997) schemes. For the surface layer parameterization the Eta geophysical fluid dynamics laboratory (GFDL) scheme (Schwarzkopf and Fels, 1991) was adopted. The Noah land surface model scheme (Chen and Dudhia, 2001) and Mellor-Yamada-Janjic (MYJ) parameterization (Janjic, 2002) were chosen as land surface and the planetary boundary layer schemes. Cumulus parameterization, namely the Kain-Fritch scheme (Kain et al., 1992), was activated only for d01 and d02.
The simulations were initialized and forced at its lateral boundaries by meteorological data derived from ERA5 reanalysis data (Hersbach and Dee, 2016) of European Center for Medium-Range Weather Forecasts (ECMWF). The reanalysis data have a spatial resolution of 0.5° × 0.5°, 37 pressure levels in the vertical direction and are provided at 6 h intervals.
Using the aforementioned setup, a series of sensitivity tests were performed in order to explore the best spin-up time for each event. Precisely, four numerical simulations were conducted for each event, starting at 24h, 18h, 12h and 6h before the initiation of the rainfall of the event .The choice of the best spin-up time for each simulation was made by comparing the temporal evolution of precipitation reproduced by WRF model to the observed precipitation by the rain gauge station at Vilia for the basin of Sarantapotamos and the rain gauge station at N. Makri for the basin of Rafina. An example of the temporal evolution of the rainfall in Vilia for event #2 is given in Fig. 2. The simulation periods for each event are presented in Table 1.

WRF -Hydro
The WRF-Hydro modeling system version 3.0 was used for this study. WRF-Hydro is a distributed hydrological modeling system which couples with WRF providing multiple physics options for surface overland flow, saturated subsurface flow, channel routing, and base-flow processes (Gochis et al., 2015). The main advantage of WRF-Hydro is 125 130 135 140 145 150 the ability to assimilate the specialized components of water cycle such as soil moisture and ground water, through the routing processes of the infiltration capacity excess and the saturated subsurface water.
In the present study, the WRF-Hydro is configured for the d04 domain, in a coupled manner with physics options of surface flow, sub-surface flow and channel routing activated. The surface and subsurface runoff are computed in a fine resolution grid (95m) that permits to better resolve the local topography. The soil infiltration and redistribution is computed in 4 layers (0-10, 10-40, 40-100, and 100-200 cm) in the fine resolution grid and then is aggregated in the coarser grid of d04.
Subsurface lateral flow of soil is calculated by applying the methodology proposed by Wigmosta et al. (1994) and Wigmosta and Lettenmaier (1999) and is computed at the high resolution grid prior to the routing of overland flow, allowing the exfiltration from fully saturated grid cells to be added to the surface flow of the coarser grid. The effects of topography and the saturation depth of soil are included in the calculation of subsurface flow. Thus, when the depth of ponded water on a grid cell exceeds a threshold, the overland flow is solved with a diffusive wave formulation adapted from Julien et al. (1995) and Ogden (1997).

Calibration method
WRF-hydro is capable to simulate soil moisture and predict the streamflow and the associated flood events after application of a suitable calibration method (e.g Givati et al., 2012;Yucel et al., 2015). The calibration process is essential in order to predict with reasonable accuracy the runoff in a sub-basin, while it is required to repeat the procedure for each sub-basin separately. The aim of the calibration is to improve the spatial resolution of parameters that control the total water volume and the shape of the hydrograph. There are two categories of calibration processes for WRF-Hydro: the stepwise (e.g. Li et al., 2017) and the automate calibration process (e.g. Cuntz et al., 2016). The stepwise approach of calibration is recommended in order to minimize the high number of model runs and the measured data which are required for the automate calibration approach.
WRF-Hydro has numerous tabulated parameters that influence the streamflow and the associated discharge, while Yucel et al. (2015) showed that four parameters are the most critical. Thus, in this study, calibration procedure was based on the stepwise method suggested by Yucel et al. (2015), which has been implemented by other authors also (e.g. Li et al., 2017, Naabil, 2017. The stepwise calibration was performed in two basic steps: firstly, we defined the parameters that influence the total water volume and then we calibrated the parameters controlling the shape of the hydrograph. The parameters that control the total water volume, are the runoff infiltration factor (REFKDT) and the surface retention depth (RETDEPRTFAC). The REFKDT parameter controls the amount of water that flows into the channel network, while the RETDEPRTFAC influences the surface slope and thus the accumulation of the water. The parameters that control the shape of the hydrograph, are related to the surface (OVROUGHRT) and channel roughness (Manning's roughness, MannN). The parameters are abbreviated following the nomenclature of WRF-Hydro namelist. The calibrated values for each parameter are shown in Table 3. In the stepwise calibration method, sensitivity tests were performed for each parameter and when a parameter is calibrated its optimum value was kept constant when the sensitivity tests for the next parameter were performed. Further details on the calibration of the aforementioned parameters for each basin (Sarantapotamos & Rafina) are given in the following section.

. Calibration of Sarantapotamos basin
Due to limited availability of streamflow data, the calibration process was performed only for event #2 at the sub-basin of Sarantapotamos, while the rest of the events were used to evaluate the performance of the calibration process. Fig. 3 shows the evolution of the discharge (observed and simulated) for each calibrated parameter. The choice of the optimum value for each parameter was based on the selected objective criteria, the Nash-Sutcliffe efficiency and the correlation coefficient (R), between simulated and observed discharges. Fig. 3a shows the results for the first parameter of the step-wise calibration method (REFKDT). As possible values for the REFKDT parameter range from 0.5 to 5, we firstly performed several simulations for possible REFKDT's values of 1, 2, 3, 4 and 5 (not shown) in order to find the appropriate range of the scaling factor. Thus, the appropriate range of REFKDT was found to be from 0.5 to 1.5 and then additional simulations were performed within this range with increment of 0.1. Fig. 3a shows that the discharge decreases as the REFKDT's values increase. For the selection of the optimum value of each parameter we implemented two basics steps. Firstly, a visual comparison of the simulated and observed discharge was performed. Secondly, we applied statistical analysis tests. More precisely, the statistical analysis included the computation of correlation coefficient and the Nash-Sutcliffe coefficient between the observed and predicted discharge calculated per 15 min for each possible value of REFKDT (Fig. 4). Thus, the value which has the best correlation/ Nash-Sutcliffe coefficient was chosen as the optimum value, after the visual comparison of the simulated and observed discharge. Namely, the value of 0.5 for REFKDT parameter was selected. Table 4 has the correlation and the Nash-Sutcliffe coefficient for the optimum value for each parameter It is noted that there is a lag at the time of maximum discharge between the observations and the model results. This discrepancy is attributed to the time lag between the simulated and observed temporal evolution of precipitation at Vilia station ( Fig. 2). After the implementation of cross correlation analysis, it was found that the maximum correlation between the simulated and the observed temporal evolution of precipitation is achieved with a delay of 5 hours. It must be noted that the results of the statistical analysis presented in Table 4 are computed after the displacement of the temporal evolution of the simulated discharge. This displacement of 5-h was necessary in order to derive the optimum value of each parameter. For instance, if we do not take into account the 5-h gap, the correlation and the Nash-Sutcliffe coefficients are not in the acceptable limits, thus the choice of the optimum value for each parameter cannot be determined. Figs 3c and 3d show the temporal evolution of discharge for the parameters which control the hydrograph shape (OVROUGHRTAC and Manning's roughness). The OVROUGHRTAC parameter is related to the surface roughness of the channel and was calibrated for values between 0.1 and 1.0 with 0.1 increments (Fig. 3c). Finally, a scaling factor value of 0.4 for OVROUGHRTAC parameter was selected.

Validation of the calibration of Sarantapotamos basin
After the calibration of WRF-Hydro over Sarantapotamos basin based on the event #2, the four parameters defined above were validated for the events #1 and #3 of Sarantapotamos basin. Figures 5a and 6a show the comparison of the temporal distribution of the observed and simulated discharges for the events #1 and #3, respectively. For the event #1, the simulated temporal distribution of the discharge shows similarity to the observed one (Fig. 5a), as the time of maximum occurrence almost coincide while the two temporal distributions don't show similar maximum values of discharge (the observed discharge is 12.8 m 3 /s and the simulated is 5.7 m 3 /s).
The correlation coefficient of the two temporal distributions is 0.83. For the event #3, the simulated and observed temporal distribution of the discharges show similarity in the time of the maximum values but the simulated discharge underestimates the observed one throughout the duration of the event (Fig. 6a), as the maximum value of the simulated discharge is 10.6 m 3 /s while the observed one is 7 m 3 /s. This is due to the underestimation of the simulated rainfall at the station of Vilia compared to the observed one (Fig. 6b). The correlation coefficient between the simulated and observed discharges is 0.75.

Calibration of Rafina basin
The stepwise calibration method suggested above, was implemented for the calibration of Rafina basin using event #5. Fig. 7 shows the temporal distribution of the precipitation as observed at the station of N. Makri and simulated using WRF atmospheric only simulations and WRF-Hydro coupled simulations, while Fig. 8 shows the temporal evolution of the observed and simulated discharges for the possible values of each calibrated parameter. The observed and simulated precipitation (provided by WRF-Hydro) are highly correlated (correlation coefficient: 0.83) while quantitatively they also compare very well (Fig. 7). The choice of the optimum values for each parameter was based on the visual comparison of the simulated and observed discharge (Fig. 8) and statistical analysis (Table 5), as it was explained for Sarantapotamos basin. The stepwise calibration method suggested above, was implemented for the calibration of Rafina basin using event #5. Fig. 7 shows the temporal distribution of the precipitation as observed at the station of N. Makri and simulated using WRF atmospheric only simulations and WRF-Hydro coupled simulations, while Fig. 8 shows the temporal evolution of the observed and simulated discharges for the possible values of each calibrated parameter. The observed and simulated precipitation (provided by WRF-Hydro) are highly correlated (correlation coefficient: 0.83) while they also compare very well quantitatively (Fig. 7). The choice of the optimum values for each parameter was based on the visual comparison of the simulated and observed discharge (Fig. 8) and statistical analysis (computation of correlation coefficient and the Nash-Sutcliffe coefficient between the observed and predicted discharge calculated per 15 min; Table 5).
In consistency to the calibration of Sarantapotamos, we firstly performed several simulations for possible REFKDT's values between 1 to 5 and we also found that the appropriate range of the scaling factor from 0.5 to 1.5. Thus, the additional simulations were performed within this range with increment of 0.1 and it was selected the value of 0.5 as optimum value for REFKDT parameter. The simulations for RETDEPRTFAC were performed within the range from 0 to 10, with increment of 1. As in the case of Sarantapotamos, the simulated discharge is decreasing with increasing values of RETDEPRTFAC (Fig. 7b). After the comparison of the selected statistical criteria, the optimum value for the RETDEPRTFAC parameter was found the value 6.
Regarding the parameters controlling the shape of the hydrograph, 10 (from 0.1 to 1.0 with increment of 0.1) and 16 (from 0.6 to 2.1 with increment of 0.1) simulations performed for the parameters related to the surface and channel roughness, respectively. After the computation of correlation coefficient and Nash-Sutcliffe parameter for each simulation, the optimum values of 0.3 and 1.2 for OVROUGHRTAC and MannN parameters were selected. At the end of the calibration procedure, the two temporal distributions (observed and discharge) have correlation coefficient 0.62, while the Nash-Sutcliffe is close to 0.5 (Table 5).

Validation of the calibration of Rafina Basin
The validation of the calibration process of Rafina basin was held by comparing the temporal distributions of simulated and observed discharges of the events #4, #6 and #7 (Figs 9b, 10b and 11b), using the optimum values of the calibration's parameters. The correlation coefficients between the simulated and observed discharges are 0.77, 0.86 and 0.62, respectively. Therefore, it is obvious that WRF-Hydro is capable to forecast the discharge after the calibration process. The simulated discharge is dependent on the simulated precipitation, thus a possible underestimation of the simulated discharge is influenced by a possible underestimation of the precipitation. For instance, at event #4, the maximum simulated discharge is 5.0 m 3 /s while the observed one is 8.0 m 3 /s (Fig. 9b). This is attributed to the underestimation of the total precipitation, as the total simulated precipitation is 27.6 mm while the observed is 37.0 mm.
Besides, the lag between the observed and simulated discharge is attributed to the lag of the observed and simulated precipitation (Fig. 9a).

Precipitation
In this section the influence of the use of the coupled model (WRF-Hydro) on the improvement of the precipitation forecast skill as compared to the atmosphere-only simulations performed with WRF model will be investigated.
Namely, WRF-Hydro contributes to a better simulation of the soil moisture content, due to the computation of the lateral redistribution and re-infiltration of the water (Gochis et al., 2013). The improved simulation of the soil moisture affects the computation of the sensible and latent heat fluxes, which influence humidity and temperature in the lower atmosphere and consequently precipitation (Seneviratne et al., 2010). Therefore, the physical process of the coupling of land-atmosphere is expected to improve the forecast skill of precipitation.
Figs 2, 5a, 6a, 7, 9a, 10a and 11a show the temporal distribution of the precipitation observed and simulated by WRF only and WRF-Hydro for each studied event observed in Sarantapotamos and Rafina basins. In all cases, the precipitation reproduced by WRF-Hydro has differences compared to WRF (atmospheric only) simulations. The temporal distribution of WRF-Hydro and WRF follow the same pattern as this is reflected in the same calculated correlation coefficients shown in Table 6. WRF-Hydro performs better than the WRF in terms of quantitative precipitation forecasting and this is reflected to the lower calculated Root Mean Square Errors and the lower Mean Absolute Error (MAE), which have been computed on hourly values of precipitation (Table 6). It must be noted that for events #4 and #7 despite the fact that the correlation coefficient is low, due to the lag between simulated and observed discharge (Figs. 9b and 11b), the values of total amount of the simulated and observed precipitation are similar. Also, the low correlation coefficient and the high MAE at event #2 is attributed to the time lag between the simulated and observed temporal evolution of precipitation (Fig. 2). Fig. 12 shows the difference between the total amount of precipitation observed minus the total amount of precipitation simulated by a) WRF-Hydro and b) WRF-only for each event. Therefore, values close to zero mean that the total amount of precipitation simulated is close to the observed one. For each case, the difference between the total amount of