Stochastic Hybrid Event Based and Continuous Approach to Derive Flood Frequency Curve

: This study proposes a methodology that combines the advantages of the event-based and continuous models, for the derivation of the maximum ﬂow and maximum hydrograph volume frequency curves, by combining a stochastic continuous weather generator (the advanced weather generator, abbreviated as AWE-GEN) with a fully distributed physically based hydrological model (the TIN-based real-time integrated basin simulator, abbreviated as tRIBS) that runs both event-based and continuous simulation. The methodology is applied to Peacheater Creek, a 64 km 2 basin located in Oklahoma, United States. First, a continuous set of 5000 years’ hourly weather forcing series is generated using the stochastic weather generator AWE-GEN. Second, a hydrological continuous simulation of 50 years of the climate series is generated with the hydrological model tRIBS. Simultaneously, the separation of storm events is performed by applying the exponential method to the 5000- and 50-years climate series. From the continuous simulation of 50 years, the mean soil moisture in the top 10 cm (MSM10) of the soil layer of the basin at an hourly time step is extracted. Afterwards, from the times series of hourly MSM10, the values associated to all the storm events within the 50 years of hourly weather series are extracted. Therefore, each storm event has an initial soil moisture value associated (MSM10Event). Thus, the probability distribution of MSM10Event for each month of the year is obtained. Third, the ﬁve major events of each of the 5000 years in terms of total depth are simulated in an event-based framework in tRIBS, assigning an initial moisture state value for the basin using a Monte Carlo framework. Finally, the maximum annual hydrographs are obtained in terms of maximum peak-ﬂow and volume, and the associated frequency curves are derived. To validate the method, the results obtained by the hybrid method are compared to those obtained by deriving the ﬂood frequency curves from the continuous simulation of 5000 years, analyzing the maximum annual peak-ﬂow and maximum annual volume, and the dependence between the peak-ﬂow and volume. Independence between rainfall events and prior hydrological soil moisture conditions has been proved. The proposed hybrid method can reproduce the univariate ﬂood frequency curves with a good agreement to those obtained by the continuous simulation. The maximum annual peak-ﬂow frequency curve is obtained with a Nash–Sutcliffe coefﬁcient of 0.98, whereas the maximum annual volume frequency curve is obtained with a Nash– Sutcliffe value of 0.97. The proposed hybrid method permits to generate hydrological forcing by using a fully distributed physically based model but reducing the computation times on the order from months to hours.


Introduction
The estimation of floods for return periods higher than the length of the observed data series and the associated peak-flow frequency curves is often needed for the design of hydraulic infrastructures, analysis of hydrological safety, urban and rural planning, estimation of flood areas, among others. This estimation is a concern in many countries due to the high economic and social consequences associated to it. Even though the risk assumed is low, the associated risk should be recalculated, as the legal regulations, climate conditions, basin, infrastructure and region development, among others, may vary along time [1]. In the standard engineering approach, the design flood for a given return period (characterized by its peak-flow and volume) is analysed with deterministic methods. The widespread use of this methodology is due to its ease of application and moderate information requirements. However, scientific advances in the field of hydrology combined with an increase in computational capacity, the development of geographic information systems and spatial modelling processes have enabled a new framework based on the application of stochastic models for the extreme characterization of the hydrological behaviour of basins. Many of the involved variables in hydrology have a stochastic nature [2,3]. In recent decades, numerous authors proposed probabilistic approaches accounting for the randomness associated with different variables [2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Several authors obtained accurate maximum peak-flow frequency curves within a Monte Carlo framework [5,[15][16][17][18][19]. Focusing on hydrologic design, a distinction should be made between the basic evaluation criteria for the selection of floods and the calculation methods for developing these criteria. The former is generally defined by the risk that a given community is willing to accept in the face of possible failure or collapse of a hydraulic structure, or possible rural and urban areas flooded with a certain recurrence associated and depends not only on technical and economic considerations but also on social, political, cultural and environmental ones. As far as flood estimation methods are concerned, there are numerous classification criteria based on: the origin of the information to be considered (based on rainfall or flow data), calculation methodology (deterministic or probabilistic), expected results (peak-flows or hydrographs), and temporal scope (continuous or eventbased methods), among others. However, they do not specify whether the return period is associated with the peak-flow, the volume of the hydrograph or the entire hydrograph. The return periods corresponding to the design flood hydrographs are commonly associated to the return period of the rainfall events that generate hydrographs with peak-flows with the same return period. Many authors questioned this hypothesis [18,[20][21][22][23]. It is also well-known that return periods associated to peak-flows are not the same as those associated to hydrograph volume.
To characterize the floods, several methods can be applied and are mainly classified according to two main groups: (a) statistical flood frequency analysis: statistical methods are generally based on processing existing local and regional flow through statistical procedures such as distribution fitting [24,25], making appropriate use of historical references where available, and (b) derived flood frequency simulation methods: by the use of a hydrological model, a set of hydrographs is obtained from which a flood frequency distribution can be derived. Statistical methods need large flow records that are hardly available [26,27], have the drawback of the uncertainty associated to the distribution fitting for large return periods [28], and provide a value of peak-flow, volume, or duration but not the hydrograph shape. Furthermore, when they are applied in ungauged basins, the physical processes that occur in the watershed are not usually considered and the uncertainty on the flood quantiles estimation increases [29]. Due to these reasons, derived flood frequency methods are generally preferred over statistical ones [30].
Derived flood frequency methods can be divided into two approaches: continuous simulation and event-based methods [30]. Using a continuous distributed hydrological model provides the user the possibility of deriving flood frequency distributions from the continuous hourly streamflow obtained at any desired point in the drainage network of the basin by forcing the distributed hydrological model with an hourly stochastic weather generator. This approach has the advantage of estimating the variables for the entire period of simulation. However, continuous models tend to be more complex than event-based models, with computational efforts that could be very intensive (mainly when the time series simulated are extended to hundreds or thousands of years and the size of the basin modelled increases to hundreds or thousands km 2 ), even when using high performance computing and parallelization processes [31][32][33]. On the other hand, event-based simulations require much shorter simulation times. However, the main disadvantage of the event-based models is the uncertainty in the estimation of the existing soil moisture within the basin, prior to the storm event to be simulated, which is crucial in the derivation of the corresponding hydrograph. In addition, it is usually assumed that the flood hydrograph has the same return period as the storm event, which could be not realistic [20][21][22][23]. Furthermore, different properties of storm events could affect the derivation of flood frequency curves: rainfall temporal distribution, event duration, maximum intensity, and total storm depth, among others. Examples of both derived flood frequency approaches can be found in the existing literature. In the case of event-based approaches, some studies have combined non-complex stochastic storm generators [3,13] or complex rainfall generators [34]; with semidistributed [35] or distributed models [36]. Some authors highlighted the importance of considering simultaneously both, the antecedent soil moisture and extreme rainfall in flood frequency analysis [37,38]. For example, Wright et al. [37] suggested that these analyses should focus on flood estimation techniques that incorporate a realistic site-specific representation of the spatial and temporal structure of storm rainfall. Moreover, Cea and Fraga [38] proposed a methodology that accounts for the observed intra-event spatial and temporal variability of rainfall, as well as for site-specific antecedent soil moisture conditions. In the case of continuous simulations, in order to reduce the computational cost, most authors have worked with lumped or semidistributed models [39,40], and some of them have worked with distributed models by using high performance computers [31,33,41].
To overcome the limitations in the two approaches, some authors have proposed combining event-based models with continuous models. Paquet et al. [14] proposed the SCHADEX method, which consists of replacing stochastic rainfall events within a short continuous simulation with observed data and a Monte Carlo framework. Li et al. [30] developed the hybrid CE approach, which consists of combining continuous long-term simulations of rainfall and short continuous hydrological simulations to probabilistically characterize the rainfall and initial soil moisture (respectively) to force an event-based model. However, both approaches used lumped models.
Potential future water scarcity may lead in the medium and long term to the need of planning new dams in some countries, while in others, scheduling proper maintenance plans for the existing ones, accounting for a framework for assessing the potential risks of failure [42,43].
The objective of this study is to propose a methodology that combines the advantages of the event-based and continuous models, for the derivation of the maximum flow and maximum volume frequency curves, by combining a stochastic continuous weather generator (the advanced weather generator, abbreviated as AWE-GEN) with a fully distributed physically based hydrological model (the TIN-based real-time integrated basin simulator, abbreviated as tRIBS) that allows both event simulation and continuous simulation. The proposed hybrid method permits to generate hydrological forcing for extreme return periods by using a fully distributed physically based model reducing the computation times on the order from months to hours.

Case Study
The hybrid method was applied to Peacheater Creek (herein, PCH). This basin is located within Baron Fork at Eldon River basin (Oklahoma, USA) and it was selected since it was originally used to test, calibrate, and validate the tRIBS model [44,45]. PCH basin has a drainage area of 64 km 2 , 16.6 km long stream, with elevations ranging between 432 and 248 m.a.s.l. (Figure 1). The basin is characterized by gentle (slopes between 2% and 5%) to steep (slopes between 15% and 40%) hillslopes. The vegetation is a mixture of oak-hickory-pine forest grasslands (mainly located in the southern part of the basin) and cropland-urban areas (mainly in the northern region). The impervious fraction of the basin is estimated on 1.6%. The predominant soil consists of gravelly silt loams [44,45]. The mean annual rainfall is 1142 mm, being spring and summer the rainiest seasons of the year. Additional details on the climate, hydrology, and basin characteristics have been presented in Ivanov et al. [44,45]; and Vivoni et al., [46]. it was originally used to test, calibrate, and validate the tRIBS model [44,45]. PCH basin has a drainage area of 64 km 2 , 16.6 km long stream, with elevations ranging between 432 and 248 m.a.s.l. (Figure 1). The basin is characterized by gentle (slopes between 2% and 5%) to steep (slopes between 15% and 40%) hillslopes. The vegetation is a mixture of oakhickory-pine forest grasslands (mainly located in the southern part of the basin) and cropland-urban areas (mainly in the northern region). The impervious fraction of the basin is estimated on 1.6%. The predominant soil consists of gravelly silt loams [44,45]. The mean annual rainfall is 1142 mm, being spring and summer the rainiest seasons of the year. Additional details on the climate, hydrology, and basin characteristics have been presented in Ivanov et al. [44,45]; and Vivoni et al., [46].

Set Up of Modelling Experiments
Two different modelling experiments were carried out: stochastic weather generation and hydrological simulations. First, we generated stochastically 5000 years of punctual hourly weather by applying the model AWE-GEN [47], previously calibrated by Gabriel-Martín et al. [48] with the climate data at Westville, a weather station located within the basin. Second, we carried out hydrological simulations by applying the tRIBS model. tRIBS is a physically-based, distributed hydrologic model developed at the Ralph M. Parsons Laboratory, Massachusetts Institute of Technology. The model emphasizes the dynamic relationship between a partially saturated vadose zone and the land surface response to the continuous storm and inter storm cycle. One-dimensional infiltration in the surface normal direction is redistributed by both the lateral fluxes in the vadose zone and in the phreatic aquifer during storm and inter storm periods. By modelling these detailed processes, tRIBS produces runoff through a number of generation mechanisms including infiltration-excess, saturation from below and perched saturation-excess runoff, among others. In combination, these processes allow tRIBS to model the development of partial contributing areas within the basin. The runoff produced from these areas is routed to the outlet using a simplified scheme that allows non-linearity in the routing through the channel network. A detailed description of the model can be found in Ivanov et al. [44] and in http://vivoni.asu.edu/tribs.html (last accessed on 12 July 2021).

Set Up of Modelling Experiments
Two different modelling experiments were carried out: stochastic weather generation and hydrological simulations. First, we generated stochastically 5000 years of punctual hourly weather by applying the model AWE-GEN [47], previously calibrated by Gabriel-Martín et al. [48] with the climate data at Westville, a weather station located within the basin. Second, we carried out hydrological simulations by applying the tRIBS model. tRIBS is a physically-based, distributed hydrologic model developed at the Ralph M. Parsons Laboratory, Massachusetts Institute of Technology. The model emphasizes the dynamic relationship between a partially saturated vadose zone and the land surface response to the continuous storm and inter storm cycle. One-dimensional infiltration in the surface normal direction is redistributed by both the lateral fluxes in the vadose zone and in the phreatic aquifer during storm and inter storm periods. By modelling these detailed processes, tRIBS produces runoff through a number of generation mechanisms including infiltration-excess, saturation from below and perched saturation-excess runoff, among others. In combination, these processes allow tRIBS to model the development of partial contributing areas within the basin. The runoff produced from these areas is routed to the outlet using a simplified scheme that allows non-linearity in the routing through the channel network. A detailed description of the model can be found in Ivanov et al. [44] and in http://vivoni.asu.edu/tribs.html (last accessed on 12 July 2021).
We generated a triangular irregular network (TIN) of 6680 nodes derived from a US Geological Survey (USGS) 30-m DEM using the procedure described in Vivoni et al. [49] ( Figure 2a). The simulations were based on a model calibration conducted by Ivanov et al. [44,45]. To make the experiment approachable, the continuous simulation (for validation purposes) was carried out by using parallel computing techniques. The basin was partitioned into eight different parts using a surface-flow partitioning script [33] with the graph-partitioning software METIS [50], which balances the number of TIN nodes across processors minimizing the dissections that occur in the channel network ( Figure 2b). Each core of the processor calculates tRIBS variables within each part of the basin and, every time step each core of the processor sends messages to the others to combine all the calculations done within the basin. The experiment was carried out by using Magerit, a high-performance computer owned by the Technical University of Madrid. An Intel family processors architecture was used, consisting of 41 nodes with two intel processors Intel XEON E5-2670 of eight cores each, and 64 GB RAM/core. The event-based model was carried out by using a standard computer (Intel processor i7-8650U, 8 cores, 16 GB RAM/core). We generated a triangular irregular network (TIN) of 6680 nodes derived from a US Geological Survey (USGS) 30-m DEM using the procedure described in Vivoni et al. [49] ( Figure 2a). The simulations were based on a model calibration conducted by Ivanov et al. [44,45]. To make the experiment approachable, the continuous simulation (for validation purposes) was carried out by using parallel computing techniques. The basin was partitioned into eight different parts using a surface-flow partitioning script [33] with the graph-partitioning software METIS [50], which balances the number of TIN nodes across processors minimizing the dissections that occur in the channel network ( Figure 2b). Each core of the processor calculates tRIBS variables within each part of the basin and, every time step each core of the processor sends messages to the others to combine all the calculations done within the basin. The experiment was carried out by using Magerit, a highperformance computer owned by the Technical University of Madrid. An Intel family processors architecture was used, consisting of 41 nodes with two intel processors Intel XEON E5-2670 of eight cores each, and 64 GB RAM/core. The event-based model was carried out by using a standard computer (Intel processor i7-8650U, 8 cores, 16 GB RAM/core).  Figure 3 shows a general scheme of the modelling framework and methodology of the hybrid method, based on the integration and application of the mentioned models AWE-GEN and tRIBS. First, two continuous set of hourly weather forcing series were generated by using AWE-GEN: 5000-and 50-years series.   Second, we identified the independent rainfall events from the 5000-and 50-years generated weather climate series by applying the exponential method [51]. Rainfall events were identified by fixing the duration of the minimum inter-event time (MIT) that follows, or precedes, a rainfall event. Two rainfall events are independent whether there is no rainfall (considering in this study a minimum rainfall threshold of 0.25 mm/h) between them in a period equal or higher than the MIT. The exponential method assumes that the event arrival time and the time between two rainfall events follows an exponential distribution. Consequently, the duration of dry periods can be approximated by an exponential distribution with a mean that is equal to the standard deviation, and therefore the coefficient of Second, we identified the independent rainfall events from the 5000-and 50-years generated weather climate series by applying the exponential method [51]. Rainfall events were identified by fixing the duration of the minimum inter-event time (MIT) that follows, or precedes, a rainfall event. Two rainfall events are independent whether there is no rainfall (considering in this study a minimum rainfall threshold of 0.25 mm/h) between them in a period equal or higher than the MIT. The exponential method assumes that the event arrival time and the time between two rainfall events follows an exponential distribution. Consequently, the duration of dry periods can be approximated by an exponential distribution with a mean that is equal to the standard deviation, and therefore the coefficient of variation is equal to one. Once the MIT is obtained, all the rainfall events within the rainfall time-series are identified. A detailed explanation of the methodology applied to identify independent rainfall events and the exponential method can be found in Restrepo-Posada and Eagleson [51] and Gabriel-Martin et al. [48].

The Hybrid Method
Third, a continuous hydrological simulation of 50 years (forced with the 50 years series of generated weather) was performed with tRIBS [44]. From the simulation, the mean soil moisture in the top 10 cm (MSM10) of the soil layer of the basin at an hourly time step was extracted. Afterwards, we selected the MSM10 values at the beginning of each rainfall event previously identified. Thus, each rainfall event has an initial soil moisture value associated (MSM10Event). It is important to point out that we assumed the independence between rainfall events and prior hydrological soil moisture conditions. This hypothesis was previously assumed by other authors with similar purposes [14] and verified in this study. Based on all MSM10Event identified, we determine 12 probability distributions of MSM10Event, one for each month of the year.
Fourth, based on previous results of Gabriel-Martin et al. [48] the five major events of each of the 5000 years in terms of total depth were identified and simulated in tRIBS within an event-based framework, and by assigning an initial moisture state value from the 12 probability distributions of MSM10Event. The initial moisture value for each event was assigned randomly by following a probability distribution previously determined and corresponding to the month of occurrence of the event. A total of 25,000 events were simulated.
Fifth, maximum annual hydrographs were identified in terms of maximum peak-flow and volume, and the associated frequency curves were derived.

Validation of the Method
To validate the method, the results obtained by the hybrid method were compared to those obtained by deriving the flood frequency curves from the continuous simulation of 5000 years, analyzing (a) the maximum annual peak-flow and maximum annual volume, and (b) the dependence between the peak-flow and volume.
As stated, the five major events of each of the 5000 years in terms of total depth were identified and simulated in tRIBS within an event-based framework. The maximum annual hydrographs were identified in terms of maximum peak-flow and volume. The return periods (Tr) associated to the maximum annual peak-flows and to the maximum annual volumes were obtained independently, by applying the Gringorten plotting position formula [52]. Finally, an analysis on how the different storm ranks and the hybrid method affected the dependence between the peak-flow and volume was performed.

Limitations of the Methodology
The applied methodology has some limitations that should be noted. The rainfall was considered spatially uniform within the basin, as a result of using a punctual stochastic weather generator. As pointed out by Liuzzo et al. [53], due to the small size of PCH, there is not a big difference between considering the spatial and punctual rainfall. The study is focused on return periods up to 500 years. For higher return periods, the length of the generated weather forcing series should be extended to ensure that it is representative enough of the return period analyzed. The procedure did not account for climate-change effects or land-use changes but considered stationary climate conditions. However, interannual variations of the mean annual rainfall were accounted for due to the nature of AWE-GEN. We assumed that the characteristics of the storm events are independent of the prior initial soil moisture conditions of the basin. In the case of modelling frameworks that couple hydrological with atmospheric models (i.e., WRF-Hydro [54]), this assumption might not be valid due to the relationship between terrestrial hydrology and weather. Finally, the methodology is applied to one basin, which limits the extrapolation of the results obtained.

Calibration and Validation of Models
As stated, in this study we used a stochastic weather generator (AWE-GEN) and a fully distributed physically-based model (tRIBS). The calibration and validation of AWE-GEN, applied to the basin PCH, was carried out in Gabriel-Martín et al. [48]. They analyzed the ability of the AWE-GEN to reproduce rainfall extremes by comparing the stochastic rainfall frequency curves with those observed at different aggregation periods (1, 6, 12, and 24 h). They used the hourly climate data of Westville (period from 1997 to 2016), a weather station located within the basin. The variables considered were rainfall, temperature, wind speed, radiation, cloudiness, relative humidity, and atmospheric pressure. Raw data had a time resolution of 5 min. The behavior of the observed rainfall frequency curves was mostly within the range of the simulated ones. Furthermore, they validated the seasonality of the storm events. The reader is referred to Fatichi et al. [47] and Ivanov et al. [55] for further details on AWE-GEN. In this study, by applying the exponential method [51] to the 5000-years generated weather climate series, we obtained a value of MIT equal to 15 h. Afterwards, we identified the independent rainfall events from the time series of 5000and 50-years. Gabriel-Martín et al. [48] also demonstrated that, for the case of PCH, by considering the five biggest storms per year (named as Storm Rank 5) classified by their total depth can result in an accurate derivation of the flood frequency curve when using continuous simulations.
Notwithstanding the above, in this study we verified these results by applying the same methodology to the generated weather series (5000 years). Figure 4 shows that considering the five biggest storms per year (Storm Rank 5) ranked by their total depth, the maximum annual peak-flow is identified in more than 95% of the 5000 analyzed years and the maximum annual volume of the hydrographs are identified in more than 98% of the analyzed years. However, when the maximum annual intensity or maximum annual storm duration are considered as variables to determine the biggest storms in a year, 86% and 26% of the maximum annual peak-flows are identified (for Storm Rank 5) respectively (Figure 4a), and 77% and 42% of the maximum annual volumes of the hydrographs are identified (for Storm Rank 5), respectively (Figure 4b). The total duration criterion does not significantly catch the effect on the maximum values of peak-flows and hydrograph volumes. Storms that ranked the highest according to total duration criterion do not account for events with high intensity and total depth which may occur in a short-time period, deriving sometimes in higher peak-flows and volumes than storms with bigger duration but with a more stable rainfall rate. Maximum intensity criterion shows some effect on hydrograph volumes and, specially, on peak-flows. However, this effect is less significant than the one achieved with the total depth criterion. Thus, the five major events of each year in terms of total depth were selected from all the storm events within the 5000 years to perform the hybrid method simulations (total of 25,000 events analyzed).
As stated, the calibration and validation of tRIBS, applied to the basin PCH, was carried out by Ivanov et al. [44,45] as part of the distributed model intercomparison project (DMIP [56]). The period analyzed was hourly streamflow from April of 1994 to July of 2000. They obtained a correlation coefficient of 0.763 and a Nash-Sutcliffe coefficient of 0.565, which was considered as satisfactory following the guidelines exposed in Sirisena et al. [57]. count for events with high intensity and total depth which may occur in a short-time period, deriving sometimes in higher peak-flows and volumes than storms with bigger duration but with a more stable rainfall rate. Maximum intensity criterion shows some effect on hydrograph volumes and, specially, on peak-flows. However, this effect is less significant than the one achieved with the total depth criterion. Thus, the five major events of each year in terms of total depth were selected from all the storm events within the 5000 years to perform the hybrid method simulations (total of 25,000 events analyzed).

Initial Soil Moisture Analysis
A continuous simulation of 50 years was performed, forcing tRIBS with 50 years of continuous hourly weather. First, the relationship between MSM10Event and the volume of the rainfall event (VEvent) for all the rainfall events in the 50 years of simulation was analyzed ( Figure 5). If all the rainfall events (3797 events in the 50 years, MIT 15 h) are considered, the value of Kendall's tau for the pair of values MSM10Event-VEvent is 0.034, whereas if considering only the five largest values (250 events) the value is 0.063. Thus, in accordance with Paquet et al. [14] we verified that the rainfall events can be considered as independent (Kendall's tau values tend to zero) of the prior initial soil moisture conditions of the basin. The average duration and total depth of the 3797 rainfall events identified was 10.2 h and 15.2 mm. respectively. Moreover, the average duration and the total depth of the 25,000 storm events selected was 25.8 h and 78.5 mm. respectively.  As stated, the calibration and validation of tRIBS, applied to the basin PCH, was carried out by Ivanov et al. [44,45] as part of the distributed model intercomparison project (DMIP [56]). The period analyzed was hourly streamflow from April of 1994 to July of 2000. They obtained a correlation coefficient of 0.763 and a Nash-Sutcliffe coefficient of 0.565, which was considered as satisfactory following the guidelines exposed in Sirisena et al. [57].

Initial Soil Moisture Analysis
A continuous simulation of 50 years was performed, forcing tRIBS with 50 years of continuous hourly weather. First, the relationship between MSM10Event and the volume of the rainfall event (VEvent) for all the rainfall events in the 50 years of simulation was analyzed ( Figure 5). If all the rainfall events (3797 events in the 50 years, MIT 15 h) are considered, the value of Kendall's tau for the pair of values MSM10Event-VEvent is 0.034, whereas if considering only the five largest values (250 events) the value is 0.063. Thus, in accordance with Paquet et al. [14] we verified that the rainfall events can be considered as independent (Kendall's tau values tend to zero) of the prior initial soil moisture conditions of the basin. The average duration and total depth of the 3797 rainfall events identified was 10.2 h and 15.2 mm. respectively. Moreover, the average duration and the total depth of the 25,000 storm events selected was 25.8 h and 78.5 mm. respectively. Figure 5. Comparison of the initial soil moisture value associated to each rainfall event (MSM10Event) and the total depth associated to the rainfall event (VEvent) during the 50 years of continuous simulation for: all the rainfall events (blue dots) and the five biggest total depth events per year (red dots).
Second, the cumulated frequency distributions of MSM10Event values of the 50 years' continuous simulation for each month of the year were obtained ( Figure 6). These frequency distributions were used as input for the hybrid method (event-based simulations).
V Event (mm) Figure 5. Comparison of the initial soil moisture value associated to each rainfall event (MSM10Event) and the total depth associated to the rainfall event (VEvent) during the 50 years of continuous simulation for: all the rainfall events (blue dots) and the five biggest total depth events per year (red dots).
Second, the cumulated frequency distributions of MSM10Event values of the 50 years' continuous simulation for each month of the year were obtained ( Figure 6). These frequency distributions were used as input for the hybrid method (event-based simulations). Finally, within a Monte Carlo framework, to each storm event, a value of MSM10Event was associated according to the frequency distribution corresponding to the month of occurrence of the storm event. Once MSM10Event was obtained, a synthetic spatial variability state of soil moisture was generated, with a MSM10 value in the basin equal to the MSM10Event. Figure 7 shows, as an example, three synthetic initial basin states for different MSM10Event values. It should be note that MSM10Event are mean values of a fully distributed basin state, that is, the initial soil moisture condition considered as input to the hybrid method (for the event-based model mode) is spatially distributed whereas the mean value is the MSM10Event obtained from the monthly frequency distributions is a basin averaged value.  Finally, within a Monte Carlo framework, to each storm event, a value of MSM10Event was associated according to the frequency distribution corresponding to the month of occurrence of the storm event. Once MSM10Event was obtained, a synthetic spatial variability state of soil moisture was generated, with a MSM10 value in the basin equal to the MSM10Event. Figure 7 shows, as an example, three synthetic initial basin states for different MSM10Event values. It should be note that MSM10Event are mean values of a fully distributed basin state, that is, the initial soil moisture condition considered as input to the hybrid method (for the event-based model mode) is spatially distributed whereas the mean value is the MSM10Event obtained from the monthly frequency distributions is a basin averaged value. Finally, within a Monte Carlo framework, to each storm event, a value of MSM10Event was associated according to the frequency distribution corresponding to the month of occurrence of the storm event. Once MSM10Event was obtained, a synthetic spatial variability state of soil moisture was generated, with a MSM10 value in the basin equal to the MSM10Event. Figure 7 shows, as an example, three synthetic initial basin states for different MSM10Event values. It should be note that MSM10Event are mean values of a fully distributed basin state, that is, the initial soil moisture condition considered as input to the hybrid method (for the event-based model mode) is spatially distributed whereas the mean value is the MSM10Event obtained from the monthly frequency distributions is a basin averaged value.

Hybrid Method Performance
Once the initial basin states were obtained, they were coupled to each of the 25,000 storm events (5 × 5000 years), serving both as forcing for the hydrological model tRIBS simulations (event-based model mode). Thus, five flood hydrographs per year were obtained, amounting 25,000 flood hydrographs. Flood frequency curves (peak-flows and maximum hydrographs volume) and dependence (between the peak-flows and hydrographs volume) were compared to those obtained from the 5000 years of continuous simulation.

Flood Frequency Distributions
By identifying the highest annual value of peak-flow and hydrograph volume, the frequency distributions curves were estimated (both, for the hybrid method and continuous simulation). Figure 8a shows, in the case of peak-flows, there is an almost perfect agreement between the frequency curves of both procedures. The Nash-Sutcliffe coefficient has a value of 0.98 by comparing the hybrid-based values to those obtained by the continuous simulation. When analyzing the maximum annual volume frequency curve (Figure 8b), it can be seen how the hybrid method underestimates the maximum annual volumes for return periods over 100 years. However, the Nash-Sutcliffe coefficient is 0.97. Maximum differences between volume flood frequency curves for the continuous and hybrid method achieve values up to 15%.

Flood Frequency Distributions
By identifying the highest annual value of peak-flow and hydrograph volume, the frequency distributions curves were estimated (both, for the hybrid method and continuous simulation). Figure 8a shows, in the case of peak-flows, there is an almost perfect agreement between the frequency curves of both procedures. The Nash-Sutcliffe coefficient has a value of 0.98 by comparing the hybrid-based values to those obtained by the continuous simulation. When analyzing the maximum annual volume frequency curve (Figure 8b), it can be seen how the hybrid method underestimates the maximum annual volumes for return periods over 100 years. However, the Nash-Sutcliffe coefficient is 0.97. Maximum differences between volume flood frequency curves for the continuous and hybrid method achieve values up to 15%.
The hybrid method diminishes the computational effort significantly. By assuming that the same hydrological model is carried out following both approaches, and a single computing node is used with similar computational characteristics, the continuous approach over 5000 years would imply a modelling time of around 420 days. On the other hand, the hybrid approach would take 4.2 days to model 50 years in continuous mode. Moreover, each event simulation would imply on average nine seconds, that is, the simulation of 25,000 events distributed in 5000 years would take 2.6 days. Thus, the modelling time of the hybrid model would be 6.8 days, a 1.6% of the 420 days needed within the continuous approach. In this study, since the time spent in the continuous hydrological simulations (5000 years, hourly time step) was approximately 49 days (by using parallel computing techniques and the computer Magerit), the hybrid method simulations were performed within 27 h, 18 for the continuous simulation of 50 years (by using parallel computing techniques and a standard 8 cores computer) and nine hours for the eventbased simulations (by using single-mode computing and a standard 8 cores computer). The weather simulations were carried out within one day of computation.  The hybrid method diminishes the computational effort significantly. By assuming that the same hydrological model is carried out following both approaches, and a single computing node is used with similar computational characteristics, the continuous approach over 5000 years would imply a modelling time of around 420 days. On the other hand, the hybrid approach would take 4.2 days to model 50 years in continuous mode. Moreover, each event simulation would imply on average nine seconds, that is, the simulation of 25,000 events distributed in 5000 years would take 2.6 days. Thus, the modelling time of the hybrid model would be 6.8 days, a 1.6% of the 420 days needed within the continuous approach. In this study, since the time spent in the continuous hydrological simulations (5000 years, hourly time step) was approximately 49 days (by using parallel computing techniques and the computer Magerit), the hybrid method simulations were performed within 27 h, 18 for the continuous simulation of 50 years (by using parallel computing techniques and a standard 8 cores computer) and nine hours for the event-based simulations (by using single-mode computing and a standard 8 cores computer). The weather simulations were carried out within one day of computation. Moreover, Figure 9 shows the relation between the peak annual flow and the maximum annual volume of the hydrograph for both, the hybrid and the continuous method. It can be seen that each method covers similar ranges of values. Both methods also highlight that, although there is positive trend, where higher peak-annual flows derive in higher maximum annual volumes, there are combinations with different peak-flows for similar volumes, and vice versa. This is mainly observed for low to medium range values (peakflow < 400 m 3 /s and volume < 10 hm 3 ). For peak-flows over 400 m 3 /s (approximately, Tr 100), the pairs (peak-flow, volume) tend to group along a line.
Water 2021, 13, x FOR PEER REVIEW 12 of 18 Moreover, Figure 9 shows the relation between the peak annual flow and the maximum annual volume of the hydrograph for both, the hybrid and the continuous method. It can be seen that each method covers similar ranges of values. Both methods also highlight that, although there is positive trend, where higher peak-annual flows derive in higher maximum annual volumes, there are combinations with different peak-flows for similar volumes, and vice versa. This is mainly observed for low to medium range values (peak-flow < 400 m 3 /s and volume < 10 hm 3 ). For peak-flows over 400 m 3 /s (approximately, Tr 100), the pairs (peak-flow, volume) tend to group along a line. Finally, we performed a sensibility analysis of the effect of considering different storm ranks on the identification of the maximum annual peak-flow and hydrograph volume and classified by return period (Tr) ranges. Figure 10 shows that, the lower Tr of the events analyzed, the higher storm ranks are needed (number of biggest storms per year considered, total depth criteria) to ensure the inclusion of both the maximum annual peakflow and the maximum annual hydrograph volume. Agreeing with previous studies [20][21][22][23]48], it is shown that it is not necessarily true the common assumption that the same probability of occurrence adopted for the precipitation is also associated to the peak-flow and maximum volume, basically in years where events with low return period occurs. For example, for Tr ranging between 1 and 10 years, if rank 1 is considered, 47% of maximum peak-flows and 58% of maximum volume are achieved for the hybrid method. Analogously, for the continuous model, 51% and 66% of maximum peak-flows and volumes, respectively, are obtained. However, by focusing on extreme return period values, this criterion, usually applied in the field of professional practice, is reasonably fulfilled. Results also highlight that the behavior of the hybrid method is similar to the continuous one by considering Storm Rank 2 or higher. For both methods, hybrid and continuous, Rank 5 ensures to achieve, at least, 95% of maximum peak-flows or maximum volumes for low range of return periods (1 < Tr < 10). For higher Tr, this percentage is obtained for Rank 3. Finally, we performed a sensibility analysis of the effect of considering different storm ranks on the identification of the maximum annual peak-flow and hydrograph volume and classified by return period (Tr) ranges. Figure 10 shows that, the lower Tr of the events analyzed, the higher storm ranks are needed (number of biggest storms per year considered, total depth criteria) to ensure the inclusion of both the maximum annual peak-flow and the maximum annual hydrograph volume. Agreeing with previous studies [20][21][22][23]48], it is shown that it is not necessarily true the common assumption that the same probability of occurrence adopted for the precipitation is also associated to the peak-flow and maximum volume, basically in years where events with low return period occurs. For example, for Tr ranging between 1 and 10 years, if rank 1 is considered, 47% of maximum peak-flows and 58% of maximum volume are achieved for the hybrid method. Analogously, for the continuous model, 51% and 66% of maximum peak-flows and volumes, respectively, are obtained. However, by focusing on extreme return period values, this criterion, usually applied in the field of professional practice, is reasonably fulfilled. Results also highlight that the behavior of the hybrid method is similar to the continuous one by considering Storm Rank 2 or higher. For both methods, hybrid and continuous, Rank 5 ensures to achieve, at least, 95% of maximum peak-flows or maximum volumes for low range of return periods (1 < Tr < 10). For higher Tr, this percentage is obtained for Rank 3.

Dependence between the Peak-Flow and Hydrograph Volume
An analysis of the dependence between the maximum annual peak-flow and the maximum annual volume of the hydrographs was performed. In addition, the ability of the hybrid method to preserve that dependence was analyzed (in comparison with the continuous model). First, for each year (from a total of 5000 years), and each method (hybrid and continuous) we identified the storm ranks which the peak-flow and maximum volume (independently) are achieved (Table 1). Second, for both models and each storm rank, we identified the number of years which both the peak-flow and the maximum volume were identified (Table 1). For example, by considering only the biggest storm per year (Rank 1, total depth criteria), in 2481 and 2396 years the peak-flow and maximum volume occurred simultaneously for the hybrid and continuous method respectively, that is, the hybrid method identified the 96.5% of the total of cases identified by the continuous model (3.5% of error). If we consider the two biggest storms per year (Rank 2, total depth criteria), in 3653 and 3576 years the peak-flow and maximum volume occurred in the biggest storm or in the second one for the hybrid and continuous model respectively, that is, the hybrid method identified the 97.8% of the total of cases identified by the continuous model (2.2% of error). It should be noted that, for this study, we considered the five biggest storm per year (total depth criteria) to apply the hybrid method. Accordingly, in 234 years the peak annual flow and in 95 years the maximum annual volume occur in storms with rank higher than 5, and thus, they were not considered (from a total of 5000 years). Analyzing the continuous model, peak-flows and maximum volumes were identified in storms with rank up to 20.

Dependence between the Peak-Flow and Hydrograph Volume
An analysis of the dependence between the maximum annual peak-flow and the maximum annual volume of the hydrographs was performed. In addition, the ability of the hybrid method to preserve that dependence was analyzed (in comparison with the continuous model). First, for each year (from a total of 5000 years), and each method (hybrid and continuous) we identified the storm ranks which the peak-flow and maximum volume (independently) are achieved (Table 1). Second, for both models and each storm rank, we identified the number of years which both the peak-flow and the maximum volume were identified (Table 1). For example, by considering only the biggest storm per year (Rank 1, total depth criteria), in 2481 and 2396 years the peak-flow and maximum volume occurred simultaneously for the hybrid and continuous method respectively, that is, the hybrid method identified the 96.5% of the total of cases identified by the continuous model (3.5% of error). If we consider the two biggest storms per year (Rank 2, total depth criteria), in 3653 and 3576 years the peak-flow and maximum volume occurred in the biggest storm or in the second one for the hybrid and continuous model respectively, that is, the hybrid method identified the 97.8% of the total of cases identified by the continuous model (2.2% of error). It should be noted that, for this study, we considered the five biggest storm per year (total depth criteria) to apply the hybrid method. Accordingly, in 234 years the peak annual flow and in 95 years the maximum annual volume occur in storms with rank higher than 5, and thus, they were not considered (from a total of 5000 years). Analyzing the continuous model, peak-flows and maximum volumes were identified in storms with rank up to 20. As this study focuses on the extreme values (peak-flow and hydrograph volume), we analyzed the dependence according to different ranges of Tr. By applying the hybrid method Figure 11 shows the percentage of years (from 5000 years analyzed) that both the peak-flow and the maximum annual volume are identified for different ranges of return periods and storm ranks considered. For example, by considering only the biggest storm per year (Rank 1, total depth criteria) and events with return period ranging between 1 and 10 years, only in the 30% of the years the peak-flow and maximum volume occurred simultaneously, that is, in the same event. The percentage increases up to 93% for events with return period ranging between 100 and 1000 years. If we consider the two biggest storms per year (Rank 2, total depth criteria) the percentage increases up to 59% (1 < Tr < 10) and 100% (100 < Tr < 1000) where the peak-flow and maximum volume occurred in the biggest storm or in the second one. Results also highlights that for extreme value analysis less than the five biggest storms per year are needed to be considered, to identify both the peak-flow and the maximum volume. In this study, for return periods ranging between 100 and 1000 years only Storm Rank 2 is needed to identify the maximums in the 5000 years analyzed. Figure 11 also highlights that the behavior of the hybrid method is similar to the continuous one by considering no less than Storm Rank 2.  As this study focuses on the extreme values (peak-flow and hydrograph volume), we analyzed the dependence according to different ranges of Tr. By applying the hybrid method Figure 11 shows the percentage of years (from 5000 years analyzed) that both the peak-flow and the maximum annual volume are identified for different ranges of return periods and storm ranks considered. For example, by considering only the biggest storm per year (Rank 1, total depth criteria) and events with return period ranging between 1 and 10 years, only in the 30% of the years the peak-flow and maximum volume occurred simultaneously, that is, in the same event. The percentage increases up to 93% for events with return period ranging between 100 and 1000 years. If we consider the two biggest storms per year (Rank 2, total depth criteria) the percentage increases up to 59% (1 < Tr < 10) and 100% (100 < Tr < 1000) where the peak-flow and maximum volume occurred in the biggest storm or in the second one. Results also highlights that for extreme value analysis less than the five biggest storms per year are needed to be considered, to identify both the peak-flow and the maximum volume. In this study, for return periods ranging between 100 and 1000 years only Storm Rank 2 is needed to identify the maximums in the 5000 years analyzed. Figure 11 also highlights that the behavior of the hybrid method is similar to the continuous one by considering no less than Storm Rank 2. Figure 11. Percentage of years (from a total of 5000 years) where the annual peak-flow and maximum annual hydrograph volume are identified by considering different storm ranks and classified according to different Tr ranges. Left graph corresponds to the hybrid model and the right one to the continuous one. Black bars show the events ranging between 1 and 10 years return period (Tr), grey bars between 10 and 100 years and the light grey bars between 100 and 1000 years. Figure 11. Percentage of years (from a total of 5000 years) where the annual peak-flow and maximum annual hydrograph volume are identified by considering different storm ranks and classified according to different Tr ranges. Left graph corresponds to the hybrid model and the right one to the continuous one. Black bars show the events ranging between 1 and 10 years return period (Tr), grey bars between 10 and 100 years and the light grey bars between 100 and 1000 years.

Conclusions
A hybrid methodology has been presented. By applying an hourly weather stochastic generator and a physically based distributed hydrological model, the hybrid method accurately estimate (compared with the results of a continuous hydrological model) the frequency curves of peak-flows and maximum annual volumes. The following summarizes the key derived conclusions:

1.
Independence between rainfall events and prior hydrological soil moisture conditions has been proved within the framework of the hybrid method developed.

2.
The hybrid method reproduces the univariate flood frequency curves with a good agreement to those obtained by the continuous simulation. The maximum annual peak-flow frequency curve is obtained with a Nash-Sutcliffe coefficient of 0.98, whereas the maximum annual volume frequency curve is obtained with a Nash-Sutcliffe value of 0.97. 3.
The correlation between the occurrence of maximum annual peak-flows and volumes is preserved by the hybrid method similarly to the continuous modelling approach. An error lower than 6 % was estimated (hybrid method and Storm Rank 5) compared to the results obtained through continuous modelling.

4.
The higher return period of the events analyzed, the lower storm ranks are needed (number of biggest storm per year considered, total depth criteria) to ensure the inclusion, either the maximum annual peak-flow and/or the maximum annual hydrograph volume. 5.
The proposed hybrid method reduces the time computation of continuous simulations from 49 days (5000 years, hourly time step, high performance computer) to 27 h of computation (standard computer) for Peacheater Creek basin (64 km 2 ), 18 h to carry out the continuous simulation of 50 years and 9 h to perform the event-based simulations (25,000 events).
The promising results obtained in this paper encourage us to continue studying these issues. Following, we highlight the future research to be carried out:

•
To deepen on the analysis of uncertainty related to the initial moisture determination, analysing the correlation with other variables (basin characteristics, rainfall, other soil moisture variables, etc.).

•
To compare the combined flood frequency law (simultaneously considering peak flow and flood volume) obtained with the continuous model and the hybrid approach. To analyse the effect of the proposed approach performance on the derived frequency law of maximum levels achieved in reservoirs, as a measure of the hydrological safety of dams.

•
To analyse how flood seasonality is preserved by the proposed approach compared to continuous modelling.

•
Accounting for the capabilities of AWE-GEN to be perturbed in order to assess the impact of climate change in flood events. Funding: This research was funded by the Spanish Ministry of Science and Innovation, grant number PID2019-105852RA-I00: "Simulation of Climate Scenarios and Adaptation in Water Resources Systems (SECA-SRH)", the funds from Universidad Politécnica de Madrid in the framework of their Program "Ayudas para contratos predoctorales para la realización del doctorado en sus escuelas, facultad, centro e institutos de I+D+i", "ayudas a proyectos de I+D de investigadores posdoctorales" (Ref. number: VJIDOCUPM19AFSW) and "XVI Convocatoria de ayudas del consejo social para el fomento de la formación e internacionalización de doctorandos".
Informed Consent Statement: Not applicable.

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