Comparison of Hydrological Platforms in Assessing Rainfall-Runoff Behavior in a Mediterranean Watershed of Northern Morocco

Abstract

This research evaluates the applicability of different types of hydrological models to simulate discharge behavior scenarios in a northern Moroccan watershed, Oued Laou Watershed (OLW). In this context, an improved understanding of the runoff mechanisms through hydrological modeling of the OLW can assist in the hazard risk management and facilitate the effective planning of water resources. For that end, a multitude of hydrological models were used to perform a very efficient modelling, and a comparative approach was adopted. Comparison of the models allowed the determination of potential sources of uncertainty in hydrological modelling of a subhumid watershed. Three models (ATelier Hydrologique Spatialisé (ATHYS), Hydrologic Modeling System (HEC-HMS), and Soil and Water Assessment Tool (SWAT)) with different characteristics were employed for a continuous modelling approach. The models were calibrated and validated using observed daily rainfall and streamflow data for 4 years (2004–2008) and 3 years (2009–2011), respectively. The multi-criteria model comparison (R², NSE, RSR, and PBIAS) showed that all three models are capable of reproducing the observed flows. The SWAT model performed well over both periods (NSE = 0.76 for calibration), with an improvement in validation (NSE = 0.84). A good agreement was also observed in the HEC-HMS model outputs, with an approximately stable NSE of 0.77 and 0.78 for calibration and validation phases, respectively. The ATHYS model showed a NSE value of 0.67 during the calibration, with a decrease of 0.06 towards the validation period. The other performance criteria confirmed these findings. Additionally, results suggest that semi-distributed and conceptual hydrological models are particularly suitable for the OLW given their physical heterogeneity. Generally, the integration of these models may be suitable for water resources assessment in OLW.

1. Introduction

Water is a major component of life and a universally recognized regulator of the development of countries. For that, it must be managed and quantified by the implementation of tools that could help decision-support and management processes. Among these tools, hydrological models have emerged over the last thirty years [1,2,3,4,5,6,7], particularly models of rainfall–runoff transformation, recharge, and evaporation [8,9,10]. The primary focus of hydrological modeling is to reproduce the flows with the smallest possible error. A successful model is expected to be robust to changes in watershed conditions. A hydrological model is one of the tools used to obtain “trustworthy” hydrological data for the conception of essential water infrastructure for social and economic well-being [11,12].

The increase in computing capacity and the improvement of computer tools make the simulation of the partitioning of precipitation into runoff by mathematical models a very reactive and globally relevant subject [13,14]. A multitude of hydrological models are available to users, where distributed, semi-distributed, and global conceptual models are the most used in hydrological studies due to their adaptability in different climatic settings. This applicability generally differs from one climate setting to another and within the same climate zone. In addition, the models’ performance depends significantly on the simulation period [15,16,17,18,19,20,21,22,23,24,25]. Conceptual models have several advantages in terms of applicability, especially in watersheds in developing countries, where input and calibration data are not always available and/or accessible to hydrological model users [26,27]. Among these advantages is the reduced number of usable input data and calibration parameters. Hence, the applicability of these conceptual models to a wide range of watersheds around the world can be easily understood given their great simplicity. However, the hydrologic modeling process used by these models is subject to a number of uncertainties that prevent robust parameter estimation.

During the last decades, several conceptual hydrological models platforms have been developed and used worldwide, such as MIKE-SHE [28], HBV [29], HYPE [30], PRMS-IV [31], SWAT [32], HEC-HMS [33], and ATHYS [34]. Despite the different descriptions of the integrated processes, conceptual hydrologic model applications often face uncertainties in model structure, input data accuracy, and calibration parameter [35]. In the watersheds of Morocco, HEC-HMS [36,37], SWAT [38,39], and recently ATHYS [40,41,42] are among the models most commonly used by hydrologists. These three model use different methods to represent each component of the rainfall-runoff process, including loss (runoff volume), transformation (direct runoff), routing (channel flow), and basic flow. HEC-HMS is a physically and conceptually based semi-distributed model for simulating rainfall-runoff processes in various geographic regions. SWAT is a continuous-time, semi-distributed [43] model, and it is among the most internationally used eco-hydrological models, as demonstrated by previous review studies and a large collection of papers now exceeding 5000 studies [44]. ATHYS is a spatially distributed hydrological modelling platform that includes several production and transfer models [41]. The principal working hypotheses of these three models are detailed in the method Section 2.3.

The novelty of this study is the testing, applying, and comparison of these three hydrological models to quantify the water balance in OLW. It should also be noted that the comparison of the model performance is widely used to detect the weaknesses and advantages of hydrological models when moving from one area to another or from one climate zone to another [45,46,47,48,49]. The large Mediterranean watersheds of Morocco, including the OLW, have a major constraint the flooding and inundation [50,51]. The OLW witnessed a large number of flood events that significantly affected the economic activities of this watershed and caused considerable damage to basic infrastructure and agricultural production, which is the main source of income for most of the inhabitants in this region [52]. A major challenge for the OLW’s decision makers is to predict runoff responses to rainfall events in order to anticipate inundations and to better understand the water balance in this region. Hydrological models significantly contribute to the proper functioning and management of water resources in OLW and similar watersheds. Results of these models (e.g., runoff rates and locations) can effectively be used to predict and mitigate risk [53,54,55].

The objective of this study is to apply a continuous modelling approach to better understand uncertainty in rainfall runoff in the Mediterranean watershed (OLW) and to discuss the best hydrologic modelling approaches by comparing simulation results based on three different model platforms, namely ATHYS, HEC-HMS, and SWAT, using more comprehensive performance criteria approaches (NSE, R², RSR, and PBIAS).

2. Materials and Methods

2.1. Description of OLW Area

The OLW (area: 940 km²) is located in the northern part of Morocco (Figure 1). The OLW is surrounded by the summits of Jebel Soukna (elevation: 1800 m above sea level) and Tissouka (2180 m) to the southeast, Jebel Kelti (1928 m) to the west, and Jebel Tazoute (1800 m) to the northeast and the Mediterranean Sea to the north (Figure 1). The OLW is marked by its rugged slopes [36]. In terms of climate, the OLW is stretched over a Mediterranean sub-humid climate zone with high rainfall in winter and hot, dry summers [50].

Precipitation is the highest of all in Morocco [56], reaching up to 1602 mm/year in the southern upstream parts of the OLW and decreasing to 300 mm/year in the northern parts of the watershed (Figure 2).

Geologically, the OLW is characterized by a significant geological diversity given its location in the Rifan Chain; it is composed of formations of the internal domain (limestone ridge and Paleozoic sheet) and external domain and also of flysch sheet or shale (Figure 3A). In terms of permeability, the OLW is covered by impermeable or low-permeability facies, and only the limestone chain, the plains, and the alluvial valleys witness high infiltration [56]. The OLW is occupied by forest and moderately cultivated land (60%) and bare soil (33%), while urban, rural, and water areas occupy 2.46%, 1.78%, and 2.76% of the OLW, respectively (Figure 3B).

2.2. Data

Daily rainfall and streamflow data from 2004 to 2012, which are equivalent to 2920 days, were selected to assess the performance of three hydrological models platforms to continuously simulate streamflow patterns in the OLW. The hydroclimatic data used in this modelling approach mainly consist of precipitation, solar radiation, temperature (max, min, average daily), relative humidity, wind speed, and discharge time series recorded at two rain gauges (Kodiet Kouriren and Bab Taza; Figure 1) and at a streamflow gauge (Kodiet Kouriren; Figure 1). The data used in this project were obtained from the Loukkos Hydraulic Basin Agency in Tetouan. Furthermore, a 12.5 m spatial resolution digital elevation model (DEM) downloaded over the OLW area from the Alaska Satellite Facility was used for the geographic analysis. Land-use data from (www.globallandcover.com, accessed on 1 January 2019) with a 30 m spatial resolution were used for the daily time step as the watershed coverage. The World Soil Database (www.fao.org, accessed on 31 December 2018) was used as the input data for the soil maps.

2.3. Description of Selected Rainfall-Runoff Models

The availability of a wide variety of rainfall-runoff models provides unprecedented opportunities to model the hydrogeologic behaviors of different hydrologic units. However, the model selection is a challenging task. In this context, Melsen et al. [57] have indicated that for an appropriate selection of a hydrological model, the primary need is to clarify the purpose of the study (e.g., water resources management or flood modelling). In the same line of thought, Paul et al. [58] underlined that the selection of an adequate model should be based on its suitability for the research questions comprising the study objective, the spatiotemporal scales, as well as the landscape of the region. In this study, three structurally and physically different hydrological models, namely ATHYS [59], HEC-HMS [60], and SWAT [43], were selected to quantify runoff rates and locations in the OLW and determine the optimal model for assessment of the rainfall-runoff behavior in the OLW and similar watersheds. These models are different in their physics, structures, input/output parameters, and spatial discretization.

2.3.1. HEC-HMS

The HEC-HMS model was developed by the United States (U.S.) Army Corps of Engineers and designed to simulate continuous hydrological processes and flood events [60]. This semi-distributed model has been successfully used worldwide to simulate hydrological processes by different hydrologists [61,62,63,64,65,66].

In this study, the basin model files and meteorological model files were prepared using the HEC-GeoHMS extension tool of Arc-GIS software. The watershed is given as an interconnected system with hydrologic and hydraulic components that are combined to simulate the basin processes. Each component is representative of the factors to convert precipitation to runoff within a part of the basin, which is usually considered as the sub-basin. The SCS-CN (soil conservation service-curve number) production function is used for the estimation of the infiltration and surface drainage process in the Mediterranean watershed [67,68].

$$P_e=\frac{{\left(P-I_a\right)}^2}{P-I_a+S}$$

(1)

where P_e is the excess precipitation, P is the total precipitation, I_a is the initial losses, and S is the maximum retention potential. In the SCS-CN method, the initial losses are given by the relationship I_a = 0.2S. The retention potential S is related to the CN, which in turn can be estimated as a function of hydrologic soil group HSG, land-cover classification, and hydrologic conditions according to USDA and NEH-4 standard lookup tables [69] or by calibration with observed data [70]:

$$S=\frac{\mathrm{254,400}-254 CN}{CN}$$

(2)

The Muskingum–Cunge method was adopted for the representation of the flood propagation process in the water bodies, while the rainfall transformation in runoff was estimated by the SCS Unit Hydrograph method.

In this model, the adopted modelling approach for simulating the hydrological behavior of OLW is a lumped one without spatialization of rainfall, land-use, or soil properties. The SCS-CN function, which is known to be the most sensitive parameter, was selected to be calibrated, while the concentration time parameter of the transfer function was computed and fixed based on an empirical formula. These methods are the most used and appropriate in the Moroccan watersheds [36,37,71]. Details of the underlying principles of the HEC-HMS model can be found in [33].

2.3.2. ATHYS

The ATHYS modelling platform (http://www.athys-soft.org, accessed on 15 January 2019) [72] has been used successfully to reproduce the hydrological behavior of several Mediterranean watersheds [40,73,74]. The ATHYS platform includes a set of spatialized hydrological models associated with hydroclimatic and geographic data processing. Such advantages allowed selecting the coupling model that combines the SCS production function and the lag and route (LR) transfer function. The SCS-LR model is a distributed model with reservoirs based on a discretization of the considered watershed into regular square cells in the Mediterranean watershed. The SCS production function is used to estimate the runoff production of the cell based on the following equation:

$$R\left(t\right)=C\left(t\right).Pb\left(t\right)$$

(3)

where R(t) is the runoff of the cell (mm.h⁻¹), C(t) is a time variable runoff coefficient (%), and Pb(t) is the instantaneous precipitation (mm.h⁻¹).

The cell production is routed to the downstream based on the LR transfer function as follows:

$$q_m\left(t\right)=0 \mathrm{if} \mathrm{t} <t_0+T_m$$

(4)

$$q_m\left(t\right)=\frac{p_e\left(t_0\right)}{K_m}exp\left(-\frac{t-\left(t_0-T_m\right)}{K_m}\right).A \mathrm{if} \mathrm{t}<t_0+T_m$$

(5)

where P_e(t₀) is the effective rainfall of cell (mm) at time to, Tm is propagation time at the outlet (s), Km is a diffusion time (s), A is the cell size (m²), and q(t) is the elementary discharge (m³/s). The complete flood hydrograph is produced by adding up all the contributions of the cells at each time.

In this model, the OLW divided into regular cells with the rainfall information is spatially interpolated based on the Thiessen polygons method. The SCS-LR combined model does not take into account the physical properties of the considered watershed. The S parameter of SCS production function and flow velocity V0 of the LR transfer function were calibrated, while the remaining parameters were set based on the literature values. An extended background on the ATHYS model can be viewed at (http://www.athys-soft.org, accessed on 15 January 2019) [72].

2.3.3. SWAT

SWAT was developed by the United States Department of Agriculture (USDA). It is a physical-based spatially distributed model, which requires specific information related to the watershed’s weather conditions, soil properties, land management, cultural practices, and flows. SWAT, as a semi-distributed model, has been largely used to model runoff, water quality, and sediment in many different watersheds [75,76,77,78,79].

In this model, the simulations are performed on a continuous time scale, and the basin is divided into sub-basins, which are then divided into hydrological response units (HRU). HRU are defined as homogeneous spatial units that have similar hydrological and geomorphological properties [80]. The hydrologic cycle in SWAT is simulated based on the water balance:

$$S{W}_t=S{W}_0+{{\displaystyle \sum }}_{n=1}^t\left(R_{day}-Q_{surf}-E_a-W_{seep}-Q_{gw}\right)$$

(6)

where SW_t refers to the final soil water content (mm), SW₀ refers to the initial soil water content (mm), Rday refers to the amount of precipitation (mm), t is time (days), Q_surf represents the amount of surface runoff (mm), E_a represents the amount of evapotranspiration (mm), Wseep is the amount of water entering the vadose zone from the soil profile (mm), and Q_gw is the amount of return flow (mm).

In SWAT model, the surface runoff is simulated based on the SCS-CN function. The surface flow is routed based on kinematic wave transfer function, while the return flow is predicted assuming a shallow aquifer [81,82,83].

In this model, the OLW is divided into several HRUs. To simulate the hydrological behavior of the studied watershed, SWAT provides us with 28 parameters to set up. The large number of parameters included in the SWAT model makes its parameterizations particularly challenging. A sensitivity analysis was carried out to select sensitive parameters to be included for calibration. After 250 runs, the assessment of sensitivity of 28 different SWAT parameters revealed that six parameters were the most influential: CN (curve number, moisture condition II (dimensionless)); ESCO (soil evaporation compensation factor (dimensionless)); SOL_AWC (available soil moisture capacity, mm h⁻¹); ALPH_BF (base flow recession constant, days); CH_N2 (Manning’s “n” value for the main channel); and CH_K2 (effective hydraulic conductivity in the channel alluvium), which were all identified to be calibrated. A full description of the SWAT model can be found in the literature [38,84].

2.4. Calibration and Validation Procedures and Performance Criteria

The process of model calibration and validation is performed manually by applying the classical split sample test. In this method, the streamflow time series is divided into two sets. The first set (2004 to 2009) is used for the calibration phase and the remaining set (2009 to 2011) for model validation. The median values of calibrated parameters for best fit between observed and estimated streamflow were used for validation.

In all three models, the simulations were performed in continuous mode with a daily time step.

2.5. Performance Criteria

In order to test the performance of the parameterized models and their ability to reproduce the flows at the outlet of the OLW, the modeled flows were assessed using four indices, and their equations are described below:

Nash and Sutcliffe (NSE)

$$NSE=1-\left[\frac{{{\displaystyle \sum }}_{i=1}^n\left(Q_i^{obs}-Q_i^{sim}\right)²}{{{\displaystyle \sum }}_{i=1}^n\left( Q_i^{obs}-Q^{mean}\right)²}\right]$$

(7)

It is estimated that the simulation is of poor quality when the Nash criterion is low (<0.5), it is acceptable when it is greater than (>0.7), and it is perfect when it is equal to (1) [85].

Linear correlation coefficient (R²)

$$R²=\frac{{{\displaystyle \sum }}_{i=1}^n(Q_{obs}-\overline{Q_{obs}})\left(Q_{sim }-\overline{Q_{sim }} \right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n(Q_{obs}-\overline{Q_{obs }} )²}\sqrt{{{\displaystyle \sum }}_{i=1}^n(Q_{sim }-\overline{Q_{sim }} )²}}$$

(8)

The closer R² is to 1, the closer the result of simulation is to the observation [86].

Root mean square error observations standard deviation ratio (RSR)

To quantify the prediction error in terms of units of the variable computed by the model, we selected the (RSR). It uses observations’ standard deviation (SD) to standardize root mean square error (RMSE). RSR is the ratio of the RMSE and SD of observed data. This is a frequently used indicator and its definition is as follows:

$$RSR=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{(Q_i^{obs}-Q_i^{sim})}^2}{{{\displaystyle \sum }}_{i=1}^N{(Q_i^{obs}-Q_{ave}^{obs})}^2}}$$

(9)

where $Q_i^{sim}$ and $Q_i^{obs}$ are the sample (of size N) containing the model estimates and the observations, respectively. It ranges from 0 to +∞, where RSR = 0 indicates a perfect fit [87].

Percent bias coefficient (PBIAS)

$$PBIAS=\left[\frac{{{\displaystyle \sum }}_{i=1}^n\left(Q_i^{obs}-Q_i^{sim}\right)\times \left(100\right)}{{{\displaystyle \sum }}_{i=1}^n\left(Q_i^{obs}\right)}\right]$$

(10)

The optimum value is 0.0, with values of low magnitude denoting accurate simulation of the model. A positive value denotes model overestimation bias, and negative values denote model underestimation bias. It is useful for continuous long-term simulations and can help identify average model simulation bias (over prediction vs. under prediction) and can incorporate measurement uncertainty [88].

Here, $Q_i^{obs} \mathrm{and} Q_i^{sim}$ are the observed and estimated flow, respectively; n is the number of data points; and ${\mathrm{Q}}^{\mathrm{mean}}$ is the mean of the observed flow peak over the investigated period.

It is significant to consider that the performance of the evaluation criteria differs depending on the targeted portion (peak, medium, and low flows) of the hydrograph [60].

3. Results and Discussion

In this study, the overall performance and robustness of the three models were assessed by analyzing the physical meaning of calibrated parameters and evaluating models performances in the calibration and validation periods.

3.1. Optimized Values of Calibrated Parameters

The optimized values of calibrated parameters for each model are shown in

Table 1. Calibrated parameter median values for the three models.

Model	Calibrated Parameters	Parameters Names	Values	Parameter Range
HEC-HMS	CN (dimensionless)	Curve number, moisture condition II	75	(69.8–93.0)
	Lag (min)	Lag time	4.5
ATHYS	CN (dimensionless)	Curve number, moisture condition II	73.57	(69.8–93.0)
	ds (1/day)	Drainage coefficient f(t) infiltration	0.12	-
	VO (m/s)	Flow velocity	1.80	-
	S (mm/mm)	Soil reservoir capacity (the potential maximum retention)	91.25	(0–+∞)
SWAT	CN2 (dimensionless)	Curve number, moisture condition II	74	(69.8–93.0)
	ESCO (dimensionless)	Soil evaporation compensation factor	0.83	(0.01–1)
	SOL_AWC (mm/mm)	Available soil moisture capacity, mm h −1	0.4	(−0.6 to 0.6)
	ALPHA_BF (1/day)	Base flow recession constant, days	0.43	(0–0.048)
	CH_N2 (s/m1/3)	Manning’s “n” value for the main channel	0.15	(0.01–0.25)
	CH_k2 (mm/h)	Effective hydraulic conductivity in the channel alluvium	58	(0–150)

. The number of sensitive parameters increases with increasing model complexity (2 for HEC-HMS, 4 for ATHYS, and 6 for SWAT). In this sense, the CN parameter showed a significant impact on the simulated flow in all three models, which confirms that the land-use and soil parameters have a significant influence on the hydrological simulations. The CN values are approximately equal in the SWAT, HEC-HMS, and ATHYS models, with calibrated values of 74, 75 and 73, 57, respectively (Table 1).

The values of the calibrated parameters in the HEC-HMS model may be related to the fact that the model was designed to consider fixed values for the physical parameters during the simulation period (CN, Lag, and Ia); however, land-use changes resulting from human activities sometimes lead to uncertainty in the CN. In addition, the CN method was originally designed for a temperate regime.

The sensitivity analysis revealed that the efficiency of the ATHYS model depends mainly on three parameters: S, ds, and V0. The S variable relies on several factures of the impermeable layer of the roads and urban regions, and changing these variables leads to the increase of the flows at the outlets [41,79].

The SWAT model employs certain parameters that the HEC-HMS and ATHYS models do not take into account (Table 1). However, the structure of the SWAT model is more complex than that of ATHYS, which is a bit complex when compared with the HEC-HMS model. This difference in complexity may raise the uncertainty of the modeled streamflow outputs. Among the six SWAT model parameters shown in (Table 1), the initial SCS-CN was found to be the most sensitive. Despite SWAT being a physical-based model, certain processes are represented by empirical functions, such as the runoff, and for this reason, the CN was quite important during simulation. The sensitivity analysis allowed decreasing the number of parameters to be calibrated for the SWAT model, which enabled avoiding the problem of over parameterization and improvement of the results of the simulations.

3.2. Analysis of the Calibration and Validation Performances

The NSE of the SWAT-derived runoff during the validation phase is higher than that of HEC-HMS and ATHYS (Figure 4). During the calibration phase, the Nash–Sutcliff coefficients of HEC-HMS and SWAT are approximately the same, while during the validation, a good improvement is detected in the Nash–Sutcliff coefficient of SWAT, which reaches 0.84. A minor Nash–Sutcliff coefficient improvement is detected in the validation phase of the HEC-HMS model, while the ATHYS model shows a decrease in performance from 0.67 during the calibration period to 0.61 during the validation period, which might be related to the structure and concept of the model. The inter-comparison of the three models shows a clear improvement in terms of discharge restitution when moving from calibration to validation for the SWAT and HEC HMS models (Figure 4).

The SWAT model performs well (NSE = 0.76) over both periods, with an improvement of 0.08 in validation. A good agreement can be observed in the HEC-HMS model outputs, with an approximately stable NSE of 0.77 and 0.78 for the calibration and validation phases, respectively. The coefficient of determination (R²) and scatterplots for all models are shown in (Figure 5). The results of the linear correlation of simulated discharge in the ATHYS model show an R² of 0.68 during the calibration with a decrease of 0.02 towards the validation. The observed drop in the ATHYS model performance based on R² and NSE criteria in the validation phase reveals that it is the least adapted to the particularities of the study area.

The results of the other performances criteria, namely RSR and PBIAS (

Table 2. RSR and PBIAS values for the calibration and validation periods.

Model	RSR		PBIAS
	Calibration	Validation	Calibration	Validation
HEC-HMS	0.53	0.51	15%	14%
ATHYS	0.51	0.53	14%	16%
SWAT	0.52	0.41	12%	9%

), confirm the findings mentioned above. The performance criteria show good results; however, they differ significantly. The SWAT model underestimates the discharge. This could be due to the fact that SWAT-derived lateral discharge relies on the relative slope, and the OLW possesses a relatively steep slope. Therefore, we may argue that the SWAT model underestimates the flow in the steep slope area in OLW. This deficiency leads to a surface runoff fraction of nearly 35% in OLW, as surface runoff and base flow compensate for the missing interflow. Soil characteristics, particularly hydraulic conductivity data, create uncertainty and cause problems in representing infiltration characteristics correctly.

Examination of Table 2 indicates that the ATHYS does not manage to improve the simulation in the validation phase, contrary to HEC-HMS and SWAT, which show a certain improvement in this phase by reduction in residue error, especially in the SWAT model. (RSR decreases from 0.52 in calibration to 0.41 in validation.)

The ATHYS model exhibits some prediction uncertainties in the simulations and in the discharge estimate. It can be observed that simply introducing more realistic soil characteristics improves the simulation, increasing the high flow peaks, especially in the driest seasons where low flow decreases in a more realistic simulation. In the ATHYS model, the continuous mode simulation is very time-consuming compared to other models (e.g., SWAT and HEC-HMS). The performance of the HEC-HMS model for both simulation periods indicates a positive trend with a small improvement in validation. In the HEC-HMS model, the watershed is presented as an interconnected system with hydrologic and hydraulic components. Several components are combined to simulate processes in the watershed; each component is representative of the factors that convert precipitation to runoff in parts of the watershed that are usually considered sub-basins. The graphical results show that the HEC-HMS model gives better results for the OLW.

The graphic comparison of the model-derived discharge allows the identification of possible sources of uncertainty in the hydrological modelling of the OLW. Figure 6 shows the difference between observed and measured streamflow for the three models (e.g., SWAT, ATHYS, and HEC-HMS) during calibration and validation phases.

The HEC-HMS model overestimates the discharge at the beginning of the calibration period, and it underestimates the discharge towards the end of the simulation during the calibration and during the validation period (Figure 6).

Despite the fact that the SWAT model underestimated the base flow, this model perfectly reproduces the observed hydrograph in term of peak flows on the OLW (Figure 6).

The ATHYS model gave satisfactory results for the two simulation periods in the studied watershed. The model underestimates the observed flow at the beginning of the calibration and overestimates the flow in the middle and towards the end of the calibration; from the graphical analysis, we can notice that the model oscillates between underestimating and overestimating of the discharge during the calibration and validation period (Figure 6).

The maximum simulated flows are often underestimated by the three models. The majority of hydrologists agree that underestimation of rainfall measurements is the main driver of uncertainty [89,90,91]. The results obtained for the three model platforms are relatively good. These results reinforced the multi-functionality and reliability of these three models as relevant tools for the management and control of water resources in Morocco.

3.3. Detailed Discussion on Models’ Performance

An overview of the performance of the three models was provided to verify which of them is most applicable to the study area. Figure 4 shows that the SWAT model performs best in terms of the NSE criterion, especially during the validation period, with the other performance criteria (R², PBIAS, and RSR) showing the same results. According to the results presented by Shekar and Vinay [92] in a sub-humid tropical watershed, SWAT performed better then HEC-HMS, while a study in the Ethiopian Rift Valley lake basin comparing HEC-HMS with the SWAT model found contradictory results [93]: they found an increase in the performance of the HEC-HMS model in the validation period, and they found that the results provided by the HEC-HMS model are more satisfactory than those provided by the SWAT model. From our results and these studies, we could see that the performance of these models relies on the soil type, elevation, and other characteristics of the watershed. The SWAT model uses characteristics such as soil type, land-use and cover, and slope, which influence the characteristics of the flow of a watershed, minimum temperature, solar radiation, precipitation, relative temperature, humidity, and wind speed for modelling various physical processes. Although the models of the SWAT platform give us good results in simulations, SWAT proves to be more adapted to the simulation of continuous scenarios of large scales of time (more than 30 years) because SWAT requires a great deal of input data. This model can be used in the OLW to study the phenomena of water erosion caused by deforestation and flash floods. SWAT can also be used for futuristic prediction of climate change, land use and land cover, and other watershed management scenarios [38].

The HEC-HMS model presents good result in the two simulations periods, similar to SWAT, while the HEC-HMS model fails to correctly estimate peak flows although the overall performance is better. The HEC-HMS model is less demanding than the SWAT model in terms of input data for runoff calibration, and the simulation results are very good; after taking into account all types of uncertainties during the model input and daily parameterization process, better correlation and good agreement between the measured data and the estimated daily flow data were shown. During the application, it is observed that the modelling procedure by HEC-HMS platform is flexible; however, some appropriate ground measurements and analysis can be used to improve the simulation results.

The most noticeable distinctions between the HEC-HMS and SWAT and ATHYS models are several properties described as follows: The HEC-HMS and SWAT models are treated as semi-distributed models that spatially interpolate rainfall inputs, and ATHYS is processed as a spatially distributed model. The simulated runoff differs between them probably due to the complexity of the model and the rationality of the parameterization. The reason for the greater uncertainty in model calibration and validation observed for ATHYS compared to SWAT and HEC-HMS can be explained by the complexity of the model structure.

The results show that the simulated peak flows are underestimated for the majority of the simulation periods, and the origin of this underestimation is probably related to the underestimation of rainfall, which is marked by significant spatial heterogeneity, which therefore represents the main difficulties in quantifying rainfall. The heterogeneity of rainfall was demonstrated by the observations of [94], who found that the more pronounced the rainfall intensity is, the smaller the affected area will be [95], confirming that the runoff production and hydrological process in a watershed are strongly influenced by the spatial variability of rainfall. Bouizrou et al. [41] demonstrated that the performance of the ATHYS model was higher when precipitation information was obtained from multiple rain gauges, suggesting that additional knowledge of the spatial variability of precipitation can improve model performance.

The difference between the coefficient of determination and the Nash–Sutcliffe coefficient for the three models is due to the underestimation of peak flows by SWAT model, underestimates in the major phases by HEC-HMS, and fluctuations between underestimation and overestimation of flow by the ATHYS model.

The advantages and limitations of these three models remain to be explored in more comprehensive applications. As discussed here, they hold considerable promise for meeting important contemporary challenges, for characterizing OLW diversity, and more widely for helping to move toward a more unified watershed-scale hydrologic theory. We emphasize that the generation and testing of multiple model alternatives must be carefully controlled. It is equally critical to avoid “over-fitting” a model setup.

4. Conclusions

In this study, we examined three rainfall-runoff patterns in the OLW located in northern Morocco. Daily rainfall data for 4 years (2004–2008) were used for calibration and 3 years (2009–2011) to validate the three structurally and physically different hydrological models. The multi-criteria model comparison shows that all three models are capable of reproducing the observed flows at calibration and validation phases. The SWAT model performs well over both periods, with an improvement in validation, which confirms its applicability in the sub-humid context. Furthermore, a good agreement can be observed in the HEC-HMS model outputs, with an approximately stable R² = 0.82 during the two simulation periods, while the ATHYS model shows a decrease in performance from 0.67 in the calibration period to 0.61 in the validation period, which might be related to the structure and concept of the model.

The precision of the three hydrological models used in simulation depends on the accuracy of the precipitation, land-use, and soil parameters. Therefore, a detailed spatial distribution and good quality of data measurements are necessary to obtain more accurate hydrological model results. However, rain gauges are scarce and unevenly distributed in the OLW due to financial and other technological constraints. However, we propose to implement more accurate estimates of surface precipitation and rain gauge stations in the near future in order to obtain specific hydrological models for the OLW. As a result, in this study, SWAT and HEC-HMS are the most appropriate models for OLW rainfall-runoff simulation because they have the best compromise between parameterization efforts and physical representativeness, followed by ATHYS.