Evaluating the Forecast Skill of a Hydrometeorological Modelling System in Greece

: A hydrometeorological forecasting system has been operating at the Institute of Marine Biological Resources and Inland Waters (IMBRIW) of the Hellenic Centre for Marine Research (HCMR) since September 2015. The system consists of the Advanced Weather Research and Forecasting (WRF-ARW) model, the WRF-Hydro hydrological model, and the HEC-RAS hydraulic–hydrodynamic model. The system provides daily 120 h weather forecasts focusing on Greece (4 km horizontal resolution) and hydrological forecasts for the Spercheios and Evrotas rivers in Greece (100 m horizontal resolution), also providing ﬂash ﬂood inundation forecasts when needed (5 m horizontal resolution). The main aim of this study is to evaluate precipitation forecasts produced in a 4-year period (September 2015–August 2019) using measurements from meteorological stations across Greece. Water level forecasts for the Evrotas and Spercheios rivers were also evaluated using measurements from hydrological stations operated by the IMBRIW. Moreover, the forecast skill of the chained meteorological–hydrological–hydraulic operation of the system was investigated during a catastrophic ﬂash ﬂood in the Evrotas river. The results indicated that the system provided skillful precipitation and water level forecasts. The best evaluation results were yielded during rainy periods. They also demonstrated that timely ﬂash ﬂood forecasting products could beneﬁt ﬂood warning and emergency responses due to their efﬁciency and increased lead time.


Introduction
The interconnected physical processes between the atmosphere and hydrosphere affect the water cycle of the planet and sometimes trigger severe hydrometeorological phenomena like floods [1]. Flooding is one of the most frequent natural hazards often posing a threat to human lives and structural utilities. As floods present an overall increasing trend worldwide during the three last decades [2], there is the necessity to increase the preparedness and the protection of socio-economic activities and human lives. Given the fact that no preventative measures or defense structures can be completely effective, flood-risk management systems should provide as much sufficient lead time as feasible. Therefore, the implementation and operation of flood forecasting and warning systems are important parts of flood management services [3][4][5].
In this context, several operational flood forecasting systems have been implemented, especially during the last two decades. These can be divided into two main categories, the global-scale and the continental-scale flood forecasting systems [6]. As regards the globalscale flood forecasting systems, the Global Flood Awareness System (GloFAS [7]) is one of the most popular. Since July 2011, GloFAS has been producing global ensemble streamflow forecasts and flood early warnings. In addition, GloFAS was recently upgraded to produce seasonal hydrological forecasts [8]. As regards the continental-scale flood forecasting systems, the National Water Model (NWM) and the European Flood Alert System (EFAS) are two of the most popular in the United States and Europe, respectively. The NWM system is a hydrological modelling framework providing forecasts from 18 h to 30 days in advance. The NWM has been operated by the National Oceanic and Atmospheric Administration (NOAA) over the entire continental United States since August 2016 [9]. The core of the NWM system is the Weather Research and Forecasting Hydrological model (WRF-Hydro [10]) set up on a 250 m × 250 m grid. WRF-Hydro has been developed by the National Center for Atmospheric Research (NCAR), also including community-based development processes. Regarding input data, it assimilates streamflow information from hydrological stations while ingesting meteorological forecasting data from a variety of sources such as the Global Forecast System (GFS) of the National Centers for Environmental Prediction (NCEP). A detailed description of the NWM system is online available at: https://water.noaa.gov/about/nwm (accessed on 12 July 2021). The EFAS system has been performing probabilistic hydrological forecasts at the Joint Research Centre (JRC) of the European Commission since 2005 [11]. It aims at providing river flood forecasts to hydrometeorological services in European countries so as to increase their preparedness. EFAS is set up on a 5 km × 5 km grid and provides forecasts from 4 h to 8 weeks in advance. EFAS uses historical and real-time hydrometeorological measurements to approximate flood thresholds and to define initial hydrological conditions, respectively. As regards meteorological information, EFAS uses weather forecasts from the German Weather Service (DWD) and the European Centre for Medium-Range Weather Forecasts (ECMWF) as well as 51 probabilistic forecasts from the Ensemble Prediction System (EPS) of the ECMWF [12].
In the Mediterranean area, catastrophic flash floods are usually characterized by their small spatiotemporal scales [13][14][15][16]. Greece also suffers from local flash floods causing severe damages, economic losses, and even fatalities [17]. Noteworthy flash flood events in Greece occurred at Mandra town on 15 November 2017 [17][18][19][20][21][22][23], at Volos city on 9 October 2006 [24], and at Rafina catchment on 11 December 2009 and 22 February 2013 [25]. The hydrometeorological forecasting and warning services in small-scale catchments demand a high resolution to capture the specific spatiotemporal characteristics of the flash flood events [26]. However, operational small-scale flash flood forecasting systems still require a considerable effort in research and development. Recently, flash flood modelling in smallscale basins of Greece has usually been based on high-resolution hydrological models forced by precipitation estimates using nowcasting approaches (e.g., [22]) regarding short-term (e.g., 1-3 h in advance) forecasts or by numerical weather prediction results (e.g., [17,18,24]) regarding forecasts exceeding even 1-2 days in advance.
In this context, the Institute of Marine Biological Resources and Inland Waters (IM-BRIW) of the Hellenic Centre for Marine Research (HCMR) has been operating and continuously upgrading a high-resolution hydrometeorological forecasting system since September 2015. The forecasts are online available from the website: https://meteo.hcmr.gr (accessed on 12 July 2021). This study aims at assessing the forecast skill of the IMBRIW's hydrometeorological system over Greece for the 4-year period from September 2015 to August 2019. Precipitation forecasts were evaluated using measurements from land surface meteorological stations across Greece included in the network of the World Meteorological Organization (WMO). Moreover, water level forecasts in the Evrotas and Spercheios rivers were evaluated using water level measurements from hydrological stations operated by the IMBRIW. The hydrological forecast skill was also investigated during a flash flood in the Evrotas river that occurred on 7 September 2016 and caused one fatality, severe damages, and the overflow of a bridge at Skala town.
The rest of the manuscript is organized as follows: In Section 2, an overview of the IMBRIW's hydrometeorological modelling system is presented describing also the setup of each model. Section 3 presents the forecast skill evaluation methodology. Section 4 presents the 4-year evaluation results regarding precipitation and river water level as well as the investigation of forecast skill in the flash flood event at the Evrotas river. The results are discussed and interpreted in Section 4. Finally, Section 5 contains the conclusions of this work.

Overview of the Hydrometeorological System
The IMBRIW's hydrometeorological system consists of atmospheric, hydrological, and hydraulic-hydrodynamic models. It is noteworthy that the integration of interdisciplinary modelling approaches better facilitates the representation and prediction of water and energy cycles, increasing the physical content of forecasting products. In particular, the system consists of the Advanced Weather Research and Forecasting (WRF-ARW) model [27], the WRF-Hydro hydrological model [10], and the HEC-RAS hydraulic-hydrodynamic model [28]. The WRF-ARW model is set up to daily produce 120-h weather forecasts for three nested domains at different horizontal resolutions (36,12, and 4 km) across the Mediterranean basin, focusing on Greece (4 km). Respectively, WRF-Hydro is set up to daily produce 120-h streamflow and water level forecasts for two river basins in Greece, i.e., the Spercheios and Evrotas river basins. It is important to note that, the capability to provide inundation forecasting using the HEC-RAS model forced by WRF-Hydro results is a recent enhancement of the system based on knowledge obtained during previous flash flood studies in Volos city [24] and Mandra town [18]. Therefore, this study only presents the application of HEC-RAS in the flash flood event that occurred on 7 September 2016 at the Evrotas river overflowing a bridge at Skala town. Figure 1 presents a flowchart which summarizes the main data used and model setup as described in the following sub-sections.

Atmospheric Model Setup
As regards the atmospheric component of the system, the WRF-ARW model has been installed and appropriately configured on the parallel computing infrastructure of the Laboratory of Hydrometeorology (LHM) of the IMBRIW-HCMR. In particular, version 3.7 of the WRF-ARW model [27] was used during the evaluation period from September 2015 to August 2019. As shown in Figure 2a For the lower boundary conditions over the sea, sea surface temperature (SST) was refreshed daily preserving it constant throughout the 5-day simulations. During the evaluation period, initial SST fields were based on the high-resolution (0.083 • × 0.083 • ) real time global (RTG) SST analyses produced by the NCEP. Furthermore, land use and topographic data were obtained from the United States Geological Survey (USGS) database. As far as the parameterization schemes are concerned, many tests were performed before September 2015 to find which of them resulted in realistic simulations of weather conditions over the area covered. Hence, the revised Monin-Obukhov scheme [29] was used during the evaluation period to represent the processes in the atmospheric surface layer. For the simulation of planetary boundary layer processes, the Yonsei University scheme [30] was used. The land surface and soil processes were simulated using the Unified Noah scheme [31]. In order to resolve the longwave and shortwave radiation processes the RRTMG scheme [32] was used, while the Thompson scheme [33] was used for cloud microphysics. Moreover, the Grell-Freitas ensemble scheme [34] was used to parameterize the convective processes in the D1 and D2 domains while explicit resolve was applied in the D3 domain. where IMBRIW operates hydrological stations are also shown. A bridge located at Skala town, which is crossed by the Evrotas river, is also depicted. Moreover, only stream orders from 2 to 6 are illustrated.

Hydrological Model Setup
In order to estimate streamflow and water level in the Spercheios and Evrotas rivers during the evaluation period, LHM applied the version 3.0 of the WRF-Hydro hydrological model [10]. The WRF-Hydro model was installed on the parallel computing infrastructure of the LHM and configured adopting an offline coupling approach using forcing fields provided by the WRF-ARW model. Forcing meteorological data used in the WRF-Hydro simulations were provided by the D3 (4 km × 4 km) domain of the WRF-ARW model. The forcing fields used are the following: liquid water precipitation rate, air temperature at 2 m, specific humidity at 2 m, incoming shortwave and longwave radiation, u-and vcomponents of wind at 10 m, and surface pressure [10]. As regards the areas covered by the hydrological simulations, the WRF-Hydro model was configured on two nested domains of WRF-ARW including Evrotas (73 × 73 Arakawa-C grid points) and Spercheios (73 × 37 Arakawa-C grid points) river basins, respectively. The nested domains used in the hydrological simulations were designed with horizontal resolution of 1 km × 1 km (Figure 2b,c). WRF-Hydro used the Noah land surface model (LSM) [35] in the 1 km × 1 km horizontal resolution to represent land processes, but routing processes were represented in a 10times finer horizontal resolution (100 m × 100 m), applying aggregation/disaggregation procedures [10]. The WRF-Hydro model estimates infiltration and exfiltration using a diffusive wave overland routing scheme [36,37] while it calculates channel routing using the Muskingum-Cunge approach [38].
Additionally, the Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) [39] data in the resolution of 90 m × 90 m distributed by the National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS) was used to construct high-resolution topographic datasets used in WRF-Hydro simulations. More specifically, the void-filled version [40] of this data obtained from the Hydrological Data and Maps Based on Shuttle Elevation Derivatives at Multiple Scales (HydroSHEDS; https://hydrosheds.cr.usgs.gov/index.php (accessed on 12 July 2021) was used, regridding it to the 100 m × 100 m Evrotas and Spercheios domains (Figure 2b,c). These topographic datasets were also used to estimate stream order classification [41] for the two river basins, which was necessary for the hydrological simulations as it implies flow accumulation information. For the channel routing setup in the two rivers, the Manning roughness coefficient, the channel bottom width, and the slide slope were set for each stream order value as demonstrated in Table 1. These channel parameter values were selected performing many simulations before 1 September 2015 in adverse and mild weather conditions. Moreover, two important parameters adjusting the overland flow before water reaching in channels were set up. In more detail, the overland flow roughness scaling factor and the initial retention height scaling factor were set to 0.3 after extensive tests in the two rivers. Given the lack of continuous streamflow measurements in the two rivers, we exploited useful information derived from previous studies using the WRF-Hydro model in Greek catchments [17,18,22,24] as well as in catchments located at adjacent countries, such as Cyprus [42], Italy [43,44], and Turkey [45,46], to set the abovementioned calibration parameters so as to give reasonable results.

Hydraulic-Hydrodynamic Model Setup
The Hydrologic Engineering Center's (CEIWR-HEC) River Analysis System (HEC-RAS) was used for river flood modelling and mapping. The free software HEC-RAS 6.0 version is capable of 2D unsteady flow simulation [47]. HEC-RAS 2D modelling and mapping had already been implemented successfully in several recent studies [18,24,[48][49][50][51][52][53][54]. Moreover, some studies tested the recent capabilities of the updated HEC-RAS [55] and the efficiency of the two-dimensional (2D) HEC-RAS model [56]. Therefore, HEC-RAS 2D was selected as the most appropriate hydraulic-hydrodynamic model to be incorporated into the IMBRIW's hydrometeorological system and, in this study, its performance is assessed for the catastrophic flash flood event of 7 September 2016 at the Evrotas river.
However, no conventional flood-related data (such as streamflow or water depth data from gauging stations close to the floodplain, and high-water marks of the most affected areas) exists for the study event. Thus, the hydraulic modelling setup and calibration were based on limited available non-conventional flood data (photographs and videos) collected from various sources (e.g., local mass-media reports). Therefore, the hydraulic modelling application was performed at the bridge's location at Skala town (Figures 2b and 3).
In 2D fluvial flood modelling, the Digital Elevation Model (DEM) accuracy is very important [57][58][59][60][61]; therefore, special care has been undertaken for the development of the area's Digital Surface Model (DSM) and Digital Terrain Model (DTM). In particular, this procedure was based on high-resolution river geometry data created by processing Unmanned Aerial Vehicle (UAV) imagery (Figure 3a,b). The spatial resolution of the area's DSM/DTM used was 11.521 cm.
The upstream boundary condition was defined using flood hydrograph (i.e., streamflow) forecast provided by the WRF-Hydro hydrological model, while the downstream boundary condition was determined as normal water depth or energy slope [28,47,62,63]. The bridge representation was based on data retrieved from the topographical survey contacted for the flood management risk plans of east Peloponnese (Figure 3c,e) [64]. The technical details of the bridge piers and abutments are presented in Figure 3d.  [64], and (f) Photographs of Skala's bridge during the flash flood event studied [65,66].
Based on the collected non-conventional flood data, i.e., photographs of Skala's bridge taken during the flash flood event (Figure 3f) [65,66], we estimated the flood water depth to be approximately 6 m at the bridge location. According to this estimation, and following a typical trial and error optimization technique, the Manning roughness coefficient was set to 0.09.

Evaluation Methodology
The precipitation forecasts were evaluated for the 4-year period from 1 September 2015 to 31 August 2019 using precipitation data measured from 24 land surface meteorological stations over Greece. These stations are included in the network of the WMO and the measurements were available from the ECMWF. Except for precipitation, these stations measure additional conventional meteorological parameters such as temperature and dew point at 2 m, atmospheric pressure, wind speed, and direction at 10 m, etc. The locations of the 24 meteorological stations are illustrated in Figure 4. Additionally, Table 2 presents information for the 24 stations used in the evaluation of precipitation forecasts. Table 2 includes information about WMO identification code (WMO ID), the latitude, the longitude, and the altitude of each meteorological station over Greece.  The evaluation methodology was based on a point-to-point comparison between simulated and measured values [67]. In more detail, the simulated 6-h accumulated precipitation values were spatially co-located with the measured ones, using the four nearest neighboring points of the D3 model grid. Similar methodology was used by [68]. It is important to note here that the measurements of precipitation are available at 00:00 and 12:00 UTC (6-h accumulated values) and at 06:00 and 18:00 UTC (12-h accumulated values). In this study, we used 6-h accumulated precipitation values that resulted from a processing of these 6-h and 12-h values. To combine model results at the four points giving the nearest points higher impact, a form of inverse distance weighting interpolation was used. Considering that p 1 , p 2 , p 3, and p 4 are the simulated 6-h accumulated precipitation values at the four nearest neighboring points around the station location, as well as r 1 , r 2 , r 3, and r 4 are the respective distances between the points and the station location, the simulated 6-h accumulated precipitation value (p) used in the evaluation is estimated as follows The evaluation scores used for the 6-h accumulated precipitation values were based on the contingency table approach [69]. As shown in Table 3, the contingency table is a   Considering the table elements, the precipitation forecast skill was evaluated, estimating the bias score (BIAS) and the equitable threat score (ETS), using pairs of simulated and measured 6-h accumulated precipitation values while considering 11 precipitation thresholds (1,3,6,9,12,15,18,21,24,27, and 30 mm).
Thus, the BIAS is defined as while the ETS is defined as where E, a term to remove the number of random correct forecasts, is defined by with N representing the total number of measurements (i.e., N = A + B + C + D) [67]. At given precipitation thresholds, the BIAS may yield a systematic overestimation (when BIAS > 1) or underestimation (when BIAS < 1), and the ETS may present poor forecasts (when ETS ≈ 0) or perfect forecasts (when ETS = 1). More details about the BIAS and ETS scores are available in [67,70]. The evaluation procedure for water level was also based on a point-to-point comparison between daily average simulated values and daily average measured values from the Sentenikos and Anthili hydrological stations of IMBRIW in the Evrotas and Spercheios rivers, respectively (Figure 2b,c). The measurements used in this study covered a period from 8 April 2016 to 31 August 2019 regarding the Sentenikos station ( Figure 2b) and a period from 4 September 2015 to 31 August 2019 regarding the Anthili station (Figure 2c). It is important to mention here that IMBRIW operates a network of hydrological stations covering various Greek catchments measuring water level among other parameters [71], while it also archives measurements in a database. The simulated water level values were averaged daily and spatially co-located with the daily average measured values. For the spatial co-location, the nearest to the stations neighboring river points of the 100 m × 100 m hydrological model grids were used. Consequently, the water level forecast skill was evaluated estimating the mean bias error (MBE [69]), the root mean square error (RMSE [69]), and the Nash-Sutcliffe efficiency (NSE [72]). The formulas of these evaluation scores are shown below where S i and M i are the daily average simulated and measured values, respectively, M is the mean of measured values, and N is the total number of measurements.

Results and Discussion
This section presents the results of forecast skill evaluation regarding precipitation over Greece and water level in the Evrotas and Spercheios rivers. Specific forecast skill characteristics of the atmospheric and hydrological models are revealed. Afterwards, this section presents the tri-model (i.e., meteorological, hydrological, and hydraulic-hydrodynamic) forecast of the Evrotas river flash flood event that occurred during a severe storm in the morning of 7 September 2016. The results of the sequential simulations are analyzed to investigate the operational forecast skill of the integrated modelling system even in such adverse weather conditions. Figure 5 demonstrates the forecast skill of the meteorological model regarding 6-h accumulated precipitation. The forecast skill is expressed by the dependence of BIAS and ETS evaluation scores to the forecast horizon up to 120 forecast hours (FH) ahead. It is important to note that only 6-h accumulated precipitation values exceeding the threshold of 1 mm were considered to show the relation between the forecast skill and the forecast horizon in weak, moderate, and heavy precipitation conditions. The BIAS score ranges from 0.73 at the 6th FH to 1.24 at the 108th FH and it is characterized by variations around one as the forecast horizon increases, implying alternations of overestimation (higher than one) and underestimation (lower than one). BIAS alternations seem to follow an afternoonevening-morning-noon pattern during the 5-day forecast horizon implying a dependence of precipitation forecast on semi-diurnal cycles. The worst BIAS (i.e., 0.73) is presented in the first 6 FH and it is related with the reduced predictability of the model during its spinup period. Overall, the BIAS score does not present systematic trends with the forecast horizon. On the other hand, the ETS presents an increase in the first 18 FH, reaching a value of 0.38, and then it gradually reduces to 0.15. Comparable ETS values and similar linear ETS reduction with the forecast horizon have also been documented by previous studies (e.g., [67,73,74]).  Figure 6 shows the dependence of BIAS and ETS on the 11 thresholds of 6-h accumulated precipitation used in the evaluation procedure. The diagrams refer to all the 24 meteorological stations and the forecasts for the first day (i.e., from the 12th to the 36th FH) were considered, excluding the first 12 FHs to avoid spinup issues. Figure 6a shows very good BIAS equal to 0.98 regarding 6-h accumulated precipitation that exceeds 1 mm, thus including almost all of the precipitation intensities (weak, moderate, and heavy). However, considering higher precipitation thresholds, BIAS ranges from 0.83 to 0.65 implying underestimation is more pronounced for heavy precipitation. Respectively, ETS ranges from 0.34 to 0.03 with the ETS values worsening for the most intense precipitation thresholds (>27 mm/6-h) (Figure 6b). Figure 6. Diagrams of (a) BIAS and (b) ETS for 6-h accumulated precipitation (mm) over all stations for the first forecast day (from the 12th to the 36th forecast hour). Y-axis represents BIAS and ETS, respectively, while X-axis represents 6-h accumulated precipitation (mm) thresholds considered in the estimation of the contingency table.

Evaluation of Precipitation and River Water Level Forecast Skill for the 4-Year Period
A seasonal analysis of precipitation forecast skill is presented in Figure 7. In more detail, Figure 7a-d shows BIAS for winter, spring, summer, and autumn, respectively, while Figure 7e-h demonstrates ETS for the same seasons. Similarly to the overall BIAS and ETS of Figure 6, 6-h accumulated precipitation values considering the 11 thresholds were used in the evaluation procedure. Precipitation underestimation and reduction of ETS with forecast horizon are evidenced in all the seasons. Nevertheless, there are seasonal variabilities of both BIAS and ETS with the worst forecast skill in summer and the best in autumn and winter. In spring, BIAS is better compared to the other seasons but ETS is worse than the others in winter and autumn. It is noteworthy that an improvement in both scores is noted in summer regarding high thresholds. This may be explained by the explicit resolve of convection used by the model, thus resulting in sufficient representation of physical processes during the formation of summer thunderstorms, which usually produce heavy precipitation. Figure 8a,b illustrates all the 24 stations included in the evaluation of the model colored to demonstrate precipitation BIAS and ETS, respectively. Similarly to Figure 6, the 6-h accumulated precipitation (mm) amounts over 1 mm were considered in the estimation of the contingency table for each station. There are spatial variabilities of BIAS and ETS scores that are attributed to the local climate of each area studied. Even though it is difficult to draw robust conclusions for the forecast skill due to the blended BIAS and ETS values presented across Greece, it is interesting to note the systematic underestimation (BIAS lower than one) illustrated at stations in western Greece. This may be attributed to the weakness of the model to represent large amounts of precipitation during the passage of barometric lows and fronts that originated from western sea areas. This weakness is partly explained by the lack of sea variabilities representation in the model forecasts. This could be feasible by implementing an air-sea coupled forecasting system "online" representing air-sea processes and feedbacks in the surface layer calculations of the meteorological model [75][76][77][78][79][80]. Such modelling systems not only improve the representation of surface fluxes but also simulate more realistically water cycle processes resulting in improvements of precipitation forecasts [17,18].   Regarding the simulated values, forecasts of daily average water level referring to the first day (i.e., from the 12th to the 36th FH) were used in the evaluation and are presented in Figure 9a,b. As shown in Figure 9a,b, the hydrological model sufficiently captures the water level variations at both river locations following the rainy-period peaks and yielding reasonable statistical scores. Regarding the Sentenikos station, MBE, RMSE, and NSE are equal to −3.52, 16.59, and 0.36, respectively. The scores are better at the Anthili station, reaching 0.15, 13.09, and 0.41, respectively. This may be attributed to the downstream location of the Anthili station compared to the upstream location of the Sentenikos station, facilitating a more complete hydrological representation of the river basin.

Evaluation of Flash Flood Forecast Skill: The Case of 7 September 2016, Evrotas River
In an effort to investigate the flood forecast skill of the integrated system in its operational application, a flash flood in Evrotas river basin was selected to be studied. The flash flood occurred on 7 September 2016 and resulted in one fatality, severe economic losses, and extended damages (transportation networks, buildings, and agricultural areas). The flash flood was triggered by a severe storm over the Evrotas river basin. As shown in Figure 10a, a cut-off trough over the central Mediterranean Sea supported southern atmospheric circulation over the study area on 7 September at 06:00 UTC. The southern circulation transferred warm and moist air from sea areas to the mountainous area around Evrotas. At the same time, the colder air in the center of the cut-off trough supported the formation of a cold front over western Greece (Figure 10b). This cold front swept through the Evrotas basin causing torrential rainfall, especially over the mountainous areas of the basin, which subsequently triggered the flash flood. The results demonstrate that the atmospheric model was capable of reasonably simulating the storm causing the flash flood. As demonstrated in Figure 11, the meteorological model produced heavy 24-h precipitation (mm) from 6 September at 21:00 to 7 September at 21:00 UTC, a period that corresponds to daily precipitation on 7 September (i.e., from 7 September at 00:00 to 8 September at 00:00) if considering local time (UTC + 3 h) to be consistent with measurements. The model seems to underestimate the daily precipitation compared to measurements from meteorological stations across the Evrotas basin or at adjacent areas ( Figure 11). In more detail, despite the fact that the model overestimated the 24-h precipitation at the Sparta (measurement: 67 mm and forecast: 85.6 mm) and the Molaoi (measurement: 38.2 mm and forecast: 61.3 mm) stations, it yielded noticeable underestimation at the Geraki (measurement: 200.2 mm and forecast: 49.3 mm) and the Kardamyli (measurement: 130 mm and forecast: 55.5 mm) stations ( Figure 11). Moreover, the model resulted in a maximum 24-h accumulated precipitation of 160 mm at the west of the river basin, however, the maximum measured amount was 200.2 mm at the Geraki station located at the east of the basin. Figure 11. Spatial distribution of 24-h accumulated precipitation (mm) resulting from the meteorological forecast for the period from 6 September at 21:00 to 7 September at 21:00 UTC (corresponding to the period from 7 September at 00:00 to 8 September at 00:00 local time, i.e., UTC + 3 h). The map also shows precipitation measurements from land surface meteorological stations at Geraki, Sparta, Molaoi, and Kardamyli for the same period. These stations are included in the network of the National Observatory of Athens (NOA) and their measurements are available at the NOA's archive database (https://meteosearch.meteo.gr/ (accessed on 12 July 2021)).
Despite the quantitative and spatial inaccuracies of the forecast, the simulated precipitation over mountainous areas, especially for a 4-h period from 01:00 to 05:00 UTC, was enough to trigger the flash flood (Figure 12a-d). Subsequently, the WRF-Hydro hydrological model captured the river flooding characteristics during the period from 05:00 to 10:00 UTC (Figure 12d-i). The overall maximum simulated streamflow occurred at 07:00 UTC reaching 767 m 3 s −1 (Figure 12f). The streamflow peak upstream from Skala's bridge (location: 36.8602 • N, 22.6677 • E) occurred at 08:00 and 09:00 UTC reaching~705 m 3 s −1 (Table 4). Streamflow values resulting from the hydrological simulation at this location were used as input hydrograph in the hydraulic-hydrodynamic simulation in order to produce the inundation forecast at Skala's bridge (Table 4).  Figure 13 presents the maximum water depth and the flood extent simulated by the 2D flood inundation model using the flood hydrograph generated from the WRF-Hydro hydrological model (Table 4). In addition, Figure 13b,c shows the maximum water depth profile for the Profile line (indicated by the green color in Figure 13a) that is upstream of the bridge and the timeseries of water depth at the deepest part of the profile line, respectively. According to the results presented in Figure 13, the simulated water depth converges to a maximum value of 6 m and there is overbank flow from both sides of the bridge. Indeed, qualitative analysis of the collected flood-related data (photographs and videos) at the Skala's bridge location, in combination with the topographical survey sketch, indicated that the maximum water depth at the specific location was around or even higher than 6 m, creating flooding on both sides of the bridge.

Conclusions
This study assesses a hydrometeorological system operating since 2015 at the Institute of Marine Biological Resources and Inland Waters (IMBRIW) of the Hellenic Centre for Marine Research (HCMR). Precipitation forecasts produced for a 4-year period (September 2015-August 2019) were evaluated using measurements from 24 meteorological stations across Greece. Water level forecasts were also evaluated using measurements from two Greek hydrological stations of IMBRIW on the Evrotas and Spercheios rivers. Furthermore, the forecast skill of the system set up in a coupled meteorological-hydrological-hydraulic operation was investigated during a flash flood in Evrotas river (7 September 2016), causing one fatality, extensive damages, and overflow of a bridge at Skala town.
The results of the statistical evaluation indicated that the use of the hydrometeorological modelling system provides skillful precipitation and water level forecasts. The forecast skill evaluation for the 4-year period indicated that the modelling system is reliable during all of the seasons, showing better scores during the rainy periods of the year (i.e., autumn and winter). However, this finding shows that the model exhibits a slightly different performance per season when it is applied with the same configuration throughout the year. Hence, if the model is configured according to the season, then it is likely to achieve better performance during its operational implementation. In several case studies, the model has been applied using different configurations from that used in the operational application, giving impressive results, consequently, this knowledge could be useful in this effort. Nevertheless, this procedure has not been applied elsewhere and needs special care, time, and huge computational resources to be accomplished. This effort will also be facilitated by newer model versions including additional features and capabilities. The spatial distribution of statistical scores used was complex and, thus, only a trend of underestimation for western Greece could be partially attributed to the lack of full coupling with dynamically-varying sea processes. As regards the Evrotas and Spercheios rivers, it is important to note that future discharge measurements or estimates could better support a more adequate evaluation of the hydrological forecasts showcasing more specific forecast skill characteristics. Nevertheless, the hydrological response of the system in respect of the water level measurements seems to be sufficient.
The system was also capable of reasonably simulating the severe storm that caused the flash flood of Evrotas river on 7 September 2016, while also simulating high streamflow values reaching~705 m 3 s −1 , which caused overflow of a bridge at Skala town. Due to the absence of official flood datasets allowing for detailed validation of critical flood characteristics for the event of 7 September 2016 at the Evrotas river, estimated flood water depth for a specific location was the main criterion used for the verification of the hydraulic-hydrodynamic simulations presented in this work. Thus, in order to examine the model performance, relevant photos and press documentary were used to assess the observed water depth in the location of Skala's bridge. The river flood modelling results demonstrated in this study show acceptable correspondence between the highest values of the simulated water depth with the empirical evidence obtained after processing information from local authorities and press at Skala's bridge location. The integrated modelling system has good potential to be used not only as operational forecasting system but also as a useful research tool in interdisciplinary flood simulation studies. The physically-based modelling system can increase the forecast lead time of weather-driven floods; therefore, it can be exploited as a forecasting guidance for flash flooding management and mitigation design.
Author Contributions: Conceptualization, formal analysis, investigation, methodology, resources, validation, visualization, writing-original draft, writing-review and editing, G.V. and G.P.; conceptualization, formal analysis, investigation, methodology, resources, validation, visualization, vriting-original draft, supervision, writing-review and editing, A.P.; conceptualization, resources, validation, supervision, writing-review and editing, E.D. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.