Sensitivity Analysis of the MOHID-Land Hydrological Model: A Case Study of the Ulla River Basin

: Hydrological models are increasingly used for studying watershed behavior and its response to past and future events. The main objective of this study was to conduct a sensitivity analysis of the MOHID-Land model and identify the most relevant parameters / processes inﬂuencing river ﬂow generation. MOHID-Land is a complex, physically based, three-dimensional model used for catchment-scale applications. A reference simulation was implemented in the Ulla River watershed, northwestern Spain. The sensitivity analysis focused on sixteen parameters / processes inﬂuencing water dynamics at that scale. River ﬂow generation was inﬂuenced by the resolution of the simulation grid, soil water inﬁltration, and crop evapotranspiration. Baseﬂow was a ﬀ ected by soil hydraulic properties, the depth of the soil proﬁle, and the dimensions of the river cross-sections. Peak ﬂows were mostly constrained by Manning’s coe ﬃ cient in the river network, as well as the dimensions of the river cross-sections. The MOHID-Land model was then used to simulate daily streamﬂow during a 10-year period (2008 − 2017). Model simulations were compared against measured data at four hydrometric stations characterizing the natural ﬂow regime of the Ulla River, resulting in coe ﬃ cients of determination (R 2 ) from 0.56 to 0.85; ratios of the standard deviation of the root mean square error to observation (RSR) between 0.4 and 0.67, and Nash and Sutcli ﬀ e model e ﬃ ciency (NSE) values ranging from 0.55 to 0.84. The MOHID-Land model thus has the capacity to reproduce watershed behavior at a daily scale with reliable accuracy, constituting a powerful tool to improve water governance at the watershed scale.


Introduction
Hydrological models are, similar to any other model, a simplified representation of the "real-world" system [1]. They are usually used for two main purposes: to predict future events and to better understand different hydrological processes [2]. Hydrological models can be classified considering their complexity [3]. Empirical models are the simplest, characterized by linear and non-linear equations relating inputs and outputs without fully representing the physical processes occurring in the catchment. Conceptual models are of intermediate complexity, based on simplified equations that describe the watershed's water balance. Finally, the most complex are physical models, which are governed by laws and equations based on real hydrological responses. These models use finite difference equations as well as state variables that can be measured and are time and space dependent [2,4]. They require a large number of parameters to describe the physical characteristics of the watershed, such as the

Model Description
The MOHID-Land model [16,27] is an open source hydrological model; the code can be accessed from an online repository (github.com/Mohid-Water-Modelling-System/Mohid). The watershed has a population of about 150,000 inhabitants, mainly located in the cities of Santiago de Compostela, Estrada, Lalín, and Padrón. The Ulla River further includes three reservoirs, namely, Portodemouros, Touro, and Bandariz ( Figure 1). Portodemouros has a total capacity of 297 hm 3 , while Bandariz and Touro are smaller, storing 2.74 hm 3 and 3.78 hm 3 , respectively. The reservoirs are mainly used for energy production and flood control. The three reservoirs function together, with Bandariz and Touro being used to uniformize turbinated flows in Portodemouros during peak times, reducing its impact in downstream areas.

Model Description
The MOHID-Land model [16,27] is an open source hydrological model; the code can be accessed from an online repository (github.com/Mohid-Water-Modelling-System/Mohid). MOHID-Land simulates the water cycle considering four compartments or mediums: atmosphere, porous media, Water 2020, 12, 3258 4 of 25 soil surface, and river network. The atmosphere is not explicitly simulated but provides data necessary for imposing surface boundary conditions that may be space and time variant. The water moves between the remaining mediums based on mass and momentum conservation equations that are computed using a finite volume approach. MOHID-Land is thus a physically based, fully distributed model using an explicit algorithm with a variable time step. The time step is maximum during dry seasons and minimum when fluxes increase.
The simulated domain can be discretized by a regular grid, quadrangular or rectangular, in the surface plane, and by a cartesian coordinate system in the vertical direction. Thus, the surface land is described using a 2D horizontal grid while the porous media is a 3D domain, which includes the same horizontal grid on the surface complemented with a vertical grid with variable thickness layers. The river network is a 1D domain defined from a digital terrain model (DTM). The water lines of the drainage network are then delineated by linking surface cell centers (nodes) together.

Infiltration
The MOHID-Land model includes three options to compute soil water infiltration. The infiltration rate (i, LT −1 ) can be first estimated according to Darcy's law, as follows: where K sat is the saturated soil hydraulic conductivity (LT −1 ), h is the soil pressure head (L), and z is the vertical space coordinate (L).
The infiltration rate can also be calculated according to the Green and Ampt method [31]: where t is the time (T), D 0 is the soil water diffusivity (L 2 T −1 ), and ∆θ is the difference between the volumetric water content in the wetted region (θ 0 ) and soil initial conditions (θ i ) (∆θ = θ 0 − θ i , L 3 L −3 ). Soil water diffusivity can then be calculated as: where K 0 is the hydraulic conductivity of the wetted region (LT −1 ), and ∆h is the difference between the matric head in the wetted region (h 0 ) and at the moving front (h F ) (∆h = h 0 − h F ). Finally, the SCS runoff curve number (CN) method [32] can be further used to estimate the infiltration rate. In this method, infiltration water is the amount of water not removed by surface runoff, entering the soil at a rate computed from the ratio between the amount of water available for infiltration and the time step. The surface runoff is first calculated as follows: where Q is the runoff (L), p is the rainfall (L), S is the potential maximum retention (L) after runoff begins, and I a is the initial abstraction (L). Initial abstraction considers all losses before runoff begins, as water retained in surface depressions, water intercepted by vegetation, evaporation, and infiltration, and is calculated as: Water 2020, 12, 3258

of 25
The S parameter is related to soil properties and land use through the curve number (CN): The CN values may range between 0 and 100 (-). Higher values of CN are related with more impermeable soils and, consequently, with higher runoff values. In MOHID-Land, the user can define average curve number values (CNII). However, because runoff is affected by soil moisture before a precipitation event, these average values are adjusted for dry (CNI) and wet (CNIII) conditions. The adjustment of the CNII values considers the comparison of the accumulated precipitation during the last five days with predefined thresholds. Under dry conditions, CN values will then decrease while under wet conditions CN values will increase. Nonetheless, it is important to note that, in MOHID-Land, the amount of water in I a may not fully infiltrate if the soil is saturated, being transformed back to surface runoff and summed to Q in Equation (4).

Surface Flow
Surface flow is computed by solving the Saint-Venant equation in its conservative form, accounting for advection, pressure, and friction forces: where Q is the water flow in the river (L 3 T −1 ), A is the cross-sectional flow area (L 2 ), g is the gravitational acceleration (LT −2 ), v is the flow velocity (LT −1 ), H is the hydraulic head (L), n is the Manning coefficient (TL −1/3 ), R h is the hydraulic radius (L), and subscripts u and v denote flow directions. In the drainage network, surface flow is solved for one direction (1D domain) considering the water lines obtained from the DTM. The cross-sections in the nodes of the river network are defined by the user. Outside the drainage network, surface flow results from the amount of water that does not infiltrate or ascends by capillarity and is solved on a 2D domain considering the directions of the horizontal grid. Water exchanges between the soil surface and the drainage network are computed according to a kinematic approach, neglecting bottom friction, and using an implicit algorithm to avoid instabilities.

Porous Media
The movement of infiltrated water in the porous media is computed by the Richards equation: where θ is the volumetric water content (L 3 L −3 ), x i represents the xyz directions (-), K is the hydraulic conductivity (LT −1 ), and S is the sink term representing root water uptake (L 3 L −3 T −1 ). The soil hydraulic properties are described using the van Genuchten Mualem functional relationships [33,34]: where S e is the effective saturation (L 3 L −3 ), θ r and θ s are the residual and saturated water contents, respectively (L 3 L −3 ), K sat is the saturated hydraulic conductivity (LT −1 ), α(L −1 ) and η (-) are empirical shape parameters, m = 1 − 1/η, and l is a pore connectivity/tortuosity parameter (-). In MOHID-Land, the relation between the horizontal and vertical hydraulic conductivities is defined by a factor (f h = K hor /K ver ) that can be adjusted by the user. The model uses the Richards equation in the whole subsurface domain and simulates saturated and unsaturated flow using the same grid. A cell is considered saturated when moisture is above a threshold value (e.g., 98%) defined by the user. When a cell reaches saturation, the model uses the saturated conductivity to compute flow and the pressure becomes hydrostatic, corrected by friction. This procedure eases the implementation of the model and simplifies its use at an annual scale. The penalty is the time step that during the wetting period must be shorter to guarantee stability. The constraint is minimized using parallel computing. The water fluxes between the porous media and the drainage network are also driven by the pressure gradient in the interface of these two mediums.

Root Water Uptake
Root water uptake considers the weather conditions and soil water contents. The reference evapotranspiration rates (ET o , LT −1 ) are first computed according to the FAO Penman-Monteith method [35]. Crop evapotranspiration rates (ET c , LT −1 ) are then obtained from the product of ET o and a single crop coefficient (K c ). The K c is imposed, with the model either assuming a constant value representing the average characteristics of each vegetation type over the entire growing season as well as average effects of evaporation from the soil [27] or a crop stage-dependent value as used in Allen et al. [35]. Advantages and limitations of these approaches are discussed in Canuto et al. [27].
The ET c values are partitioned into potential soil evaporation (E p , LT −1 ) and crop transpiration (T p , LT −1 ) as a function of the simulated leaf area index (LAI, L 2 L −2 ) [36]: where λ is the extinction coefficient of radiation attenuation within the canopy (-). The LAI values are simulated using a modified version of the EPIC model [37,38], considering the heat units for the plant to reach maturity, the crop development stages, and crop stress. Additional details can be found in Ramos et al. [17]. Root water uptake reductions, i.e., T p reductions, are finally computed using the macroscopic approach proposed by Feddes et al. [39], as follows: where T a is the actual transpiration rate (T a , LT −1 ) and α is a prescribed dimensionless function of h (0 ≤ α ≤ 1) limiting T p over the root zone in the presence of depth-varying stressors, such as water stress [40,41]. According to the linear model proposed by Feddes et al. [39], root water uptake is maximum when the pressure head is between h 2 and h 3 , has a linear reduction when h > h 2 or h < h 3 , and becomes zero when h < h 4 and h > h 1 (subscripts 1−4 denote different threshold pressure heads). Finally, the actual soil evaporation (E a , LT −1 ) is calculated from E p values by imposing a pressure head threshold value [42].

Model Set-Up (Reference Simulation)
The MOHID-Land model was applied to the study area using a constant horizontally spaced grid in the eastern and northern directions (215 columns × 115 rows), with origin on 42. The drainage network was derived from the DTM. The geometry of the river cross-sections was defined according to Andreadis et al. [44]. This database relates the drained area in each node to the top width and depth of the cross-section at that node ( Corine Land Cover [29] data with a resolution of 100 m were used to identify land use. For each land use, the vegetation type and respective Manning's coefficient were defined. Vegetation type was defined according to MOHID's vegetation database, while the Manning's coefficients were defined according to Pestana et al. [45]. The K c values were set based on vegetation type. As mentioned earlier, these coefficients represent empirical average values for the growing season and are not stage dependent as in Allen et al. [35]. The resulting interpolation to the MOHID's simulation domain ( Figure 2) showed a variation for the Manning coefficient from 0.023 to 0.298 s m −1/3 while the K c values varied from 0.15 to 1.0 ( Figure 2).
The soil domain was divided into six grid layers, with a thickness from the surface to the bottom layers of 0.3, 0.3, 0.7, 0.7, 1.5, and 1.5 m (vertical grid). On the other hand, the soil profile was characterized by three horizons: 0 to 0.6 m, 0.6 to 2.0 m, and 2.0 m to the soil bottom. Thus, the first horizon included the first two grid layers, the second horizon included the two middle grid layers, and the third horizon consisted of the two bottom layers of the vertical grid domain. The soil depth in each cell was estimated considering the cell slope, with flat areas approximating a maximum predefined value while sloped areas approached a minimum. These maximum and minimum soil depths were defined as 5.0 and 0.1 m, respectively. Soil data were extracted from the multilayered European Soil Hydraulic Database [46]. This database includes the Mualem-van Genuchten hydraulic parameters for the whole of Europe at a resolution of 250 m. Information is provided at 0, 0.05, 0.15, 0.3, 0.6, 1.0, and 2.0 m depths. For the present application, the 0.3 m layer was used to characterize the first horizon, the 1.0 m layer was used for the second horizon, and the 2.0 m layer was used for the bottom horizon. Figure 2 shows the spatial distribution of soil data in the Ulla watershed while Table 2 presents the corresponding Mualem-van Genucthen parameters. The f h parameter, which relates the horizontal and vertical hydraulic conductivities, was set to 10. For the initial conditions, the soil was assumed as saturated for 95% of the profile (from the bottom to the surface), while the soil water content in the vadose zone was set to field capacity. Soil water infiltration was computed following the Darcian approach (Equation (1)). Water 2020, 12, x FOR PEER REVIEW 8 of 25 The soil domain was divided into six grid layers, with a thickness from the surface to the bottom layers of 0.3, 0.3, 0.7, 0.7, 1.5, and 1.5 m (vertical grid). On the other hand, the soil profile was characterized by three horizons: 0 to 0.6 m, 0.6 to 2.0 m, and 2.0 m to the soil bottom. Thus, the first horizon included the first two grid layers, the second horizon included the two middle grid layers, and the third horizon consisted of the two bottom layers of the vertical grid domain. The soil depth in each cell was estimated considering the cell slope, with flat areas approximating a maximum predefined value while sloped areas approached a minimum. These maximum and minimum soil depths were defined as 5.0 and 0.1 m, respectively. Soil data were extracted from the multilayered European Soil Hydraulic Database [46]. This database includes the Mualem-van Genuchten hydraulic parameters for the whole of Europe at a resolution of 250 m. Information is provided at 0, 0.05, 0.15, 0.3, 0.6, 1.0, and 2.0 m depths. For the present application, the 0.3 m layer was used to characterize the first horizon, the 1.0 m layer was used for the second horizon, and the 2.0 m layer was used for the bottom horizon. Figure 2 shows the spatial distribution of soil data in the Ulla watershed while Table 2 presents the corresponding Mualem-van Genucthen parameters. The fh parameter, which relates the horizontal and vertical hydraulic conductivities, was set to 10. For the  θ r , residual water content; θ s , saturated water content; α and µ are empirical shape parameters; l, pore connectivity/tortuosity parameter; K sat , saturated hydraulic conductivity.
Meteorological data were extracted from the ERA5-Reanalysis dataset [47]. This dataset provides several gridded meteorological parameters with an hourly timestep and with a resolution of 0.28125 • (31 km). The variables used were the u and v components of wind velocity at 10 m height, dewpoint temperature and air temperature at 2 m height, surface solar radiation downwards, surface pressure, total cloud cover, and total precipitation. The ERA5 precipitation data were first validated through comparison with measured values obtained at Melide and Santiago meteorological stations ( Figure 1). For the Melide station, a comparison was conducted from 1 January 2008 to 31 December 2012, resulting in a coefficient of determination (R 2 ) of 0.74 ( Figure 3). For the Santiago station, that period was from 1 January 2016 to 31 December 2017, with the R 2 value reaching 0.78. The ERA5 data were then interpolated to the case study grid. θr, residual water content; θs, saturated water content; α and µ are empirical shape parameters; l, pore connectivity/tortuosity parameter; Ksat, saturated hydraulic conductivity.
Meteorological data were extracted from the ERA5-Reanalysis dataset [47]. This dataset provides several gridded meteorological parameters with an hourly timestep and with a resolution of 0.28125° (31 km). The variables used were the u and v components of wind velocity at 10 m height, dewpoint temperature and air temperature at 2 m height, surface solar radiation downwards, surface pressure, total cloud cover, and total precipitation. The ERA5 precipitation data were first validated through comparison with measured values obtained at Melide and Santiago meteorological stations ( Figure 1). For the Melide station, a comparison was conducted from 1 January 2008 to 31 December 2012, resulting in a coefficient of determination (R 2 ) of 0.74 ( Figure 3). For the Santiago station, that period was from 1 January 2016 to 31 December 2017, with the R 2 value reaching 0.78. The ERA5 data were then interpolated to the case study grid.

Sensitivity Analysis
The computational cost of one single run with the MOHID-Land model made it impossible to carry out thousands of runs to perform a global and detailed sensitivity analysis. Thus, a local sensitivity analysis was performed by modifying selected MOHID-Land parameters and processes

Sensitivity Analysis
The computational cost of one single run with the MOHID-Land model made it impossible to carry out thousands of runs to perform a global and detailed sensitivity analysis. Thus, a local sensitivity analysis was performed by modifying selected MOHID-Land parameters and processes one at the time and evaluating their impact on river flow. The respective flow-duration curves, which express the exceedance probability of a certain flow [48], were then compared against the reference simulation. A sensitivity index was finally computed following the methodology proposed by Ranatunga et al. [7].
The simulated river flows at a specific location were arranged in a descending order and ranked from 1 to N with the exceedance probability (frequency of occurrence) being given as follows: where F is the exceedance probability (expressed as a percentage of time a flow value is equaled or exceeded), R is the rank, and N is the total number of flow values resulting from the simulation. The flow duration curve is thus the graphical representation of flow values and the corresponding F values. The flow duration curves were then divided into five zones based on the percentage exceedance as shown in Figure 4. Flows with an F value from 0 to 10% were considered high flows; 10 to 40% corresponded to moist conditions; 40 to 60% corresponded to mid-range flows; 60 to 90% corresponded to dry flows; and 90 to 100% corresponded to low flows [7]. exceeded), R is the rank, and N is the total number of flow values resulting from the simulation. The flow duration curve is thus the graphical representation of flow values and the corresponding F values. The flow duration curves were then divided into five zones based on the percentage exceedance as shown in Figure 4. Flows with an F value from 0 to 10% were considered high flows; 10 to 40% corresponded to moist conditions; 40 to 60% corresponded to mid-range flows; 60 to 90% corresponded to dry flows; and 90 to 100% corresponded to low flows [7]. The sensitivity index (SI) aimed at measuring the relative influence of the analyzed parameters/processes on river flow. This index results from the normalization of the root mean square error (RMSE) by dividing the error of the estimate of each simulated scenario by the range of flow values of the reference simulation:  The sensitivity index (SI) aimed at measuring the relative influence of the analyzed parameters/processes on river flow. This index results from the normalization of the root mean square error (RMSE) by dividing the error of the estimate of each simulated scenario by the range of flow values of the reference simulation: A larger influence of a parameter/process on watershed hydrology is represented by higher values of SI, while lower values of SI mean a lower influence of the analyzed parameter/process. Model simulations were performed from 1 January 2008 to 31 December 2012 (five years). The first four months of simulations were considered as the warm-up period and were not included in the analysis of the results. Flows were analyzed at a daily scale. Due to MOHID-Land's high computational requirements and having as the ultimate goal the calibration of Ulla River watershed, the sensitivity analysis focused on a single variation of the selected parameters/processes, illustrating the impact of those variables on model performance. The analyzed parameters/processes as well as the variations imposed on model inputs were defined based on previous calibrations of the MOHID-Land model for different watersheds [22,[25][26][27]. Additionally, the calibration processes performed on similar physically based models were also considered [49][50][51][52][53][54]. Thus, sixteen parameters/processes influencing water dynamics at the catchment scale were analyzed: • The resolution of the simulation grid was modified from 0.005 • (≈500 m) used in the reference simulation to 0.01 • (≈1000 m; 140 columns × 100 rows) (simulation 1, S1).

•
The DTM was changed from the EU-DEM (30 m) to the one provided by the National Geographic Institute of Spain, with a resolution of 5 m [55]. The new DTM was interpolated to the simulation grid, with the model also delineating a new catchment area and drainage network based on that input (S2).
• The effect of cross-section geometry on river flow was assessed in two simulations. In one (S3), the top and bottom widths were increased by 25% (i.e., larger river) while the depth and the shape of the cross-section remained the same. In the other (S4), the river depth was increased by 100%, while the top and bottom widths and the shape of the cross-section were maintained as in the default simulation. Table 1 shows the variations introduced to those parameters per drainage area.

•
The K sat value of each cell was multiplied by a factor of 10 while f h was maintained (S5). As a result, the horizontal hydraulic conductivity was also modified since f h = K hor /K vert .

•
The f h value was analyzed in a separate test by changing this parameter from 10, in the reference simulation, to 20 (S6). The K sat,vert values were the same as in the default scenario, meaning that a change in f h led to an increase in the K hor . • The number of layers in the vertical grid increased from 6 to 12 as defined in Table 3 (S7), thus decreasing the thickness of the layers when compared with the reference simulation.

•
The soil depth also increased from 5 to 10 m (S8), with the number of layers in the vertical grid increasing from six to nine (Table 3).

•
The surface Manning coefficients increased by 50% when compared with the reference simulation (S9).

•
The SCS curve number method was used to compute runoff and soil water infiltration as an alternative to Equation (1) (S11). The CN values were defined for each grid cell according to the soil type and land cover. The hydrologic soil groups (HSGs) were extracted from the HYSOGs250 m dataset [56], which derived that information from the soil texture classes and depth to bedrock available in the SoilGrids250 m product [57]. That information was then merged with Corine Land Cover [29] data following the United States Department of Agriculture [58] guidelines to derive the CN values. Figure 5 presents the CN values adopted in this study. Additionally, changes in the CN values were also assessed by decreasing the values set in S11 by 25% (S12).

•
The Green and Ampt infiltration method was now used as an alternative to Equation (1) (S13).
The MOHID-Land model needed the values of K sat,ver , suction head, porosity, and wilting point in each cell as inputs ( Figure 6). These inputs were obtained by combining the information available in the LUCAS database [59] with data from Rossman [60], who related the soil texture classes with the soil hydraulic characteristics.

•
The importance of the porous media and vegetation growth processes for river flow results were investigated in three separate simulations. Firstly, vegetation growth processes were deactivated (S14), meaning that no evapotranspiration occurred in the catchment area. Secondly, both the porous media and vegetation growth processes were deactivated (S15). In the absence of porous media, the SCS CN method was used to compute the partitioning of rainfall data between surface runoff and infiltration. Infiltration water was then lost to the system since soil porous processes were not considered. The CN values presented in Figure 5 were adopted for this analysis. Lastly, both porous media and vegetation growth processes remained deactivated, but the CN values were reduced by 25% (S16) as in S12. this study. Additionally, changes in the CN values were also assessed by decreasing the values set in S11 by 25% (S12). • The Green and Ampt infiltration method was now used as an alternative to Equation (1) (S13).
The MOHID-Land model needed the values of Ksat,ver, suction head, porosity, and wilting point in each cell as inputs ( Figure 6). These inputs were obtained by combining the information available in the LUCAS database [59] with data from Rossman [60], who related the soil texture classes with the soil hydraulic characteristics.  this study. Additionally, changes in the CN values were also assessed by decreasing the values set in S11 by 25% (S12). • The Green and Ampt infiltration method was now used as an alternative to Equation (1) (S13).
The MOHID-Land model needed the values of Ksat,ver, suction head, porosity, and wilting point in each cell as inputs ( Figure 6). These inputs were obtained by combining the information available in the LUCAS database [59] with data from Rossman [60], who related the soil texture classes with the soil hydraulic characteristics.  The time required for MOHID-Land to compute each simulation was also quantified. The average, maximum, and minimum time (in seconds) required to complete each day of simulation were compared for all scenarios. This information was critical to understand which parameters/processes would improve computational speed during model calibration/validation.

Model Calibration/Validation
Based on results from the sensitivity analysis, the following parameters/processes were modified from those adopted in the reference simulation during model calibration/validation: the vertical hydraulic saturated conductivity (K sat,ver ); the relation factor between the horizontal and vertical hydraulic conductivities (f h ); and the dimensions of the cross-sections in the river network. Model calibration consisted of modifying these parameters one at a time, and adjusting them by trial-and error, until deviations between model simulations and flow measurements at Sar, Ulla, Arnego-Ulla, and Deza hydrometric stations (Figure 1) were minimized. Validation was then performed by simply comparing model simulations using the calibrated parameters with flow measurements at the same hydrometric stations. The Sar, Ulla, Arnego-Ulla, and Deza hydrometric stations were selected because their data describe the natural flow regime of the Ulla watershed, which was obtained from Augas de Galicia [61]. The Ulla-Teo station (Figure 1) was under the influence of reservoir management and therefore its data were not considered for evaluating model performance. Model simulations were carried out from 1 January 2008, to 31 December 2017 (10 years). The first four months were considered as the warm-up period, while the period between 1 May 2008 and 31 December 2012 was considered for calibration. The validation period was defined between 1 January 2013 and 31 December 2017.
Different goodness-of-fit tests were considered for assessing model performance, including the Nash and Sutcliffe model efficiency (NSE), the percent bias (PBIAS), the RMSE-observation standard deviation ratio (RSR), and the coefficient of determination (R 2 ). The NSE [62] was computed as: where Q i obs is the observed flow on day i, Q i sim is the simulated flow for day i, Q mean obs is the observed mean flow for the period under consideration, and p is the total number of days in that same period. The NSE is used to assess the relative magnitude of residual variance compared to the measured data variance and it ranges between −∞ and 1.0, being 1.0 the optimal value. Values between 0.0 and 1.0 are classified as acceptable levels of performance, and values ≤ 0.0 indicate that the mean observed value is a better predictor than the simulated value. The PBIAS is a statistical parameter that measures the average tendency of the simulated data to be larger or smaller than their observed counterparts, and was computed as follows: The optimal value of PBIAS is 0.0 and low-magnitude values indicate accurate model simulation. Positive values demonstrate model underestimation while negative values represent model overestimation.
The RSR results from the ratio between the RMSE and the standard deviation of observed values were obtained as follows: RSR incorporates the benefits of error index statistics and includes a scaling/normalization factor. RSR is equal to 0.0 when RMSE is 0.0, indicating that the variation is residual, and the model is perfect. Thus, low RSR values correspond to low RMSE values and a good model simulation performance. Finally, R 2 describes the degree of collinearity between simulated and measured data and ranges from 0 to 1, with higher values indicating less error variance. This statistical parameter was computed as follows: According to Moriasi et al. [63], model performance for streamflow can be classified as satisfactory when NSE > 0.50, RSR ≤ 0.70, PBIAS ± 25%, and R 2 > 0.5.

Impact of Model Parameters/Processes on River Flow
Although flow-duration curves were compared for all flow stations, only results for the Ulla-Teo station (Figure 1) are presented graphically, to limit the number of figures and to maintain consistency. This station was not considered for evaluating model performance because of the influence of reservoir management on river flow. However, its larger drainage area produced the most contrasting differences between the simulation scenarios and the reference simulation, which helped demonstrate model behavior during the sensitivity analysis. Figure 7 shows the flow-duration curves for each simulation included in the sensitivity analysis, while Tables 4 and 5 present the mean flow and respective SI values for each exceedance probability class.
Decreasing the resolution of the base grid from 500 m (reference simulation) to 1000 m (S1) had a substantial impact on river flow, with mean values decreasing between 71% and 97% in all ranges of the flow-duration curve ( Table 4). That was also visible in the SI values, which ranged from 0.42 to 0.93 (Table 5). These results showed that setting up the base grid, which basically defined how detailed the study area would be represented, was an important step for accurately simulating river flow at the catchment scale. A more detailed grid led to a more dynamic watershed than when using a coarser one. These results contrast with Sreedevi and Eldho [54], who, after testing three grid sizes (4000, 2000, and 1000 m), found no scale dependency on streamflow generation values using the SHETRAN model. Nevertheless, the use of a coarser resolution grid in the Ulla River watershed had a contrasting result at the hydrometric stations located in the river heads above 600 m (Sar, Arnego-Ulla, and Deza). In those locations with a drainage area smaller than at Ulla-Teo, the coarser steeper slopes in the DTM ended up promoting surface runoff and higher flow peaks, which eventually disappeared by the time the iflow reached the Ulla-Teo station.
On the other hand, for the same resolution grid (500 m), results with a more detailed DTM (S2) did not differ much from the reference simulation (Tables 4 and 5). The use of high-resolution data in coarse scale applications seems thus to be irrelevant when not accompanied by a base grid also with greater resolution to be able to consider such detailed information. Yet, improvements n streamflow simulations using more detailed grids and DTMs can only go so far, as shown by Pieri et al. [50]. These authors found no statistically significant differences in the accuracy of DTMs varying between 10 and 2 m resolution when simulating streamflow and sediment yield using the WEPP model in a considerably smaller basin than the Ulla River catchment (the Centonara catchment in northern Italy, 1.92 km 2 ). Also, Zhang et al. [51] and Nazari-Sharabian et al. [53] showed the influence of the DTM resolution on the calibration of streamflow simulations using different physically based models.
Increasing the width (S3) and depth (S4) of the cross-sections (i.e., the river network) promoted higher river flow in all the exceedance intervals except for the moist conditions (Q 10−40 ), where it slightly decreased (Figure 7). The largest increase was in the higher flow class (Q 0−10 ), i.e., the flow peaks, with S3 and S4 leading to an increase of the mean flow by 11% and 39%, respectively (Table 4). This was explained by the fact that increasing the dimensions of the cross-sections meant that the boundary between the riverbed and the porous media would also increase, promoting water exchanges between these two mediums, mainly from the porous media to the riverbed. However, only S4 resulted in higher SI values in the extremes of the flow duration curve, reaching 0.29 for the Q 90−100 class and 0.26 for the Q 0−10 class. Water 2020, 12, x FOR PEER REVIEW 15 of 25 Figure 7. Flow-duration curves for the simulations considered during sensitivity analysis. Impact of: (a) grid resolution (S1) and elevation data (S2); (b) cross-section widths and heights (S3 and S4); (c) vertical and horizontal saturated hydraulic conductivities (S5 and S6); (d) vertical soil discretization and soil depth (S7 and S8); (e) surface and channel Manning coefficients (S9 and S10); (f) infiltration methods (S11, S12, and S13); (g) porous media and vegetation processes (S14, S15, and S16) (see Section 2.4 for details). (a) grid resolution (S1) and elevation data (S2); (b) cross-section widths and heights (S3 and S4); (c) vertical and horizontal saturated hydraulic conductivities (S5 and S6); (d) vertical soil discretization and soil depth (S7 and S8); (e) surface and channel Manning coefficients (S9 and S10); (f) infiltration methods (S11, S12, and S13); (g) porous media and vegetation processes (S14, S15, and S16) (see Section 2.4 for details).
The increase in K sat,ver (S5) and f h (S6) led to a rise of river flow values in the mid-range (Q 40−60 ), low (Q 60−90 ), and dry (Q 90−100 ) classes of the flow duration curve, which varied from 48% to 188% when compared with the reference simulation. On the other hand, flow values in the Q 0−10 range (high flows) showed a decreasing trend while the Q 10−40 class (moist condition) remained basically unchanged.
Thus, the increase in those parameters promoted the infiltration process (exchange between the surface and the porous media), the faster movement of soil water (subsurface flow), and the exchange between the porous media and river sections. Higher subsurface flows caused the water to allocate to the baseflow instead of generating flow peaks, increasing the mean values in the lower classes of the flow-duration curve. As a result, the SI values were notoriously higher in the lower classes of the flow duration curve (Q 40−60 to Q 90−100 ), reaching values from 1.20 to 1.52 in S5 and from 0.48 to 1.52 in S6 (Table 5) due to the increasing baseflow. Note that a scaling factor of 10 was considered for varying K sat,ver , in line with previous calibrations of this parameter performed by Canuto et al. [27] in the Guadiana catchment, Iberian Peninsula. However, this parameter is known to be one of the most variable in nature, being affected by soil physical, chemical, and biological properties. There is also the fact that Tóth et al. [46] maps were developed from a database containing measurements of this parameter carried out in the laboratory on samples of limited size, with issues related to their representativeness being often raised when describing actual flow conditions, transport, and reaction processes occurring at the field/catchment scales due to limitations in the porous medium continuum [64,65]. Similarly, a scaling factor of 20 was considered in this application for f h . The literature is far less rich on variations of this parameter since it is specific to the MOHID-Land model. Yet, f h has been found to vary between 3.0 [27] and 25.0 [25]. The latter was used to fit simulated streamflow to measurement values in the Alegria River watershed, northern Spain. Thus, different scaling factors could have been here considered for K sat,ver and f h , but results would not differ much since the main point was to identify which hydrological processes would be affected by these parameters.
The influence of vertical discretization (S7) on river flow was practically null. This explains the relatively coarse discretization of the vertical grid in existing applications of the MOHID-Land model at the watershed scale [22,[25][26][27]. These studies reported vertical grids varying from 7 to 12 layers when describing soil profiles reaching 8 to 30 m depth. More detailed grids would have increased the computational burden with no practical benefit for model performance. Nonetheless, those numbers contrast with the one-dimensional application presented in Ramos et al. [17], who used 100 grid cells with 0.02 m thickness each to describe soil water dynamics in a soil profile with 2 m depth and respective interactions between the vadose zone and a shallow groundwater table. On the other hand, a deeper soil profile (S8) led to an increase in the lower exceedance probability classes of the flow-duration curve (Q 40−60 , Q 60−90 and Q 90−100 ) from 53% to 289% (Table 4), resulting in higher SI values in those classes (0.52 to 2.71). The increase in soil depth represented a larger volume of water stored below the surface, at depths (bottom profile) that were not affected by evapotranspiration. This water stored at deeper depths became thus an additional resource for exchanges between the porous media and the river network, leading to higher baseflow.
The impact of modifying the surface's Manning coefficients was small (S9). This agreed with Canuto et al. [27], who considered this parameter in the calibration of streamflow simulations in the Guadiana basin, but changes to the model's default values ended up being only minor. Changing the channel's Manning coefficient (S10) resulted in a decrease of the mean flow in the Q 0−10 class by 23% and an increase in the Q 90−100 class by 10% (Figure 7, Table 4). Increasing the roughness and, consequently, the friction between the water and the surface led to a decrease in flow velocity. The infiltration process was then promoted, causing an increase in baseflow and a reduction of the flow peaks. Nonetheless, the SI values were always smaller than 0.14, which show the reduced influence of this parameter in the generation of streamflow.
The use of the SCS CN method (S11) promoted the increase in mean flow values for the moist (Q 10−40 ), mid-range (Q 40−60 ), and dry flow classes (Q 60−90 ) by 8%, 19%, and 14%, respectively, when compared with the reference simulation (Table 4). Flow in the Q 40−60 class also corresponded to the highest SI value (0.20; Table 5). Reducing the CN values of all grid cells by 25% (S12) naturally led to less runoff, with the mean flow values for the same classes referred to above being 3%, 9%, and 6% higher than in the reference simulation while in the Q 90−100 class they were lower by 6%. The use of the Green and Ampt method as an alternative to Equation (1) had no real impact on river flow. Obviously, future works need to analyze the sensitivity of the MOHID-Land model to inputs used for computing soil water infiltration with this method. Notwithstanding, the MOHID-Land model offers analytical, semi-analytical, but also empirical solutions for modeling a key process in the hydrological cycle, which can be selected depending on the complexity of users' applications.
Finally, the deactivation of vegetation (S14) and porous media (S15) modules produced a strong modification of the flow-duration curve at the Ulla-Teo station (Figure 7), which was expected since most processes included in the reference simulation were disregarded. Ignoring the evapotranspiration process in the catchment (S14) led to an increase in the mean flow in all intervals of the flow duration curve, particularly in the mid-range (Q 40−60 ), dry (Q 60−90 ), and low flow (Q 90−100 ) classes (Table 4). Flow in the Q 90−100 class also returned the highest SI value of 14.32 (Table 5). Without the main driver for soil water dynamics, soil water contents as well as exchanges between the porous media and the river network simply increased, promoting mainly baseflow. The porous media module was responsible for the occurrence of subsurface flow and baseflow. This constituted a major component in the water budget of the Ulla River catchment. Disregarding this module naturally led to a decrease in the mean river flow values in all classes of the flow-duration curve because the infiltrated water simply disappeared from the system. With the reduction of the CN values by 25% (S16), that drop was even higher (Table 4).  Table 6 summarizes those results by presenting the minima, mean, and maxima values. The fastest computation time was naturally obtained by running the simplest applications, i.e., S2 with the coarser grid resolution of 1000 m, and S15 and S16 without considering the porous media processes. The former can be explained by a simulation grid four times smaller than the one in the reference simulation, substantially decreasing the calculations needed to run the simulation, and the latter by disregarding subsurface flow, with MOHID-Land distancing from its physically based nature and relying on a more empirical basis. On the other hand, the longest simulation performance was achieved with a more detailed vertical grid domain (S7). The higher number of vertical layers increased the calculations related to the porous media processes, with the model spending more time to perform the simulation. However, as showed earlier, this had no impact on river flow. All other simulations presented a similar computation time, with no real influence on the number of days (seven days) the reference simulation took to run. computation time, with no real influence on the number of days (seven days) the reference simulation took to run.   238  402  1764  S1  29  83  550  S2  190  319  2011  S3  228  389  1513  S4  255  399  1269  S5  213  422  1947  S6  237  407  1702  S7  354  528  2829  S8  303  447  2155  S9  235  404  1752  S10  234  399  1723  S11  216  360  1711  S12  221  359  1600  S13  231  334  1599  S14  209  309  1437  S15  6  65  475  S16  5  52  448 Simulations: 1-grid resolution; 2-elevation data; 3 and 4-cross-section widths and heights; 5 and 6-vertical and horizontal saturated hydraulic conductivities; 7 and 8-vertical soil discretization and soil depth; 9 and 10-surface and channel Manning coefficients; 11, 12, and 13-infiltration methods; 14, 15, and 16-porous media and vegetation processes (see Section 2.4 for details)  Table 6. Minimum, mean, and maximum computation times for each simulation test.

Minimum Mean Maximum
Reference Simulation  238  402  1764  S1  29  83  550  S2  190  319  2011  S3  228  389  1513  S4  255  399  1269  S5  213  422  1947  S6  237  407  1702  S7  354  528  2829  S8  303  447  2155  S9  235  404  1752  S10  234  399  1723  S11  216  360  1711  S12  221  359  1600  S13  231  334  1599  S14  209  309  1437  S15  6  65  475  S16  5  52  448 Simulations: 1-grid resolution; 2-elevation data; 3 and 4-cross-section widths and heights; 5 and 6-vertical and horizontal saturated hydraulic conductivities; 7 and 8-vertical soil discretization and soil depth; 9 and 10-surface and channel Manning coefficients; 11, 12, and 13-infiltration methods; 14, 15, and 16-porous media and vegetation processes (see Section 2.4 for details)  Table 7. These results were obtained after modifying K sat,ver , f h , and the dimensions of the cross-sections in the river network defined for the reference simulation. The K sat,ver values for each soil horizon in Table 2 were multiplied by a scaling factor of 10, similarly as performed in the sensitivity analysis. The f h value was then automatically updated since this parameter represents the relation between horizontal and vertical conductivities in each cell. The dimensions of the cross-sections were defined as shown in Table 1 and consisted of increasing river depth and interactions between the porous media and the river network. Thus, by modifying these parameters, the calibration process mainly focused on increasing baseflow in the Ulla River watershed.  Table 7. These results were obtained after modifying Ksat,ver, fh, and the dimensions of the cross-sections in the river network defined for the reference simulation. The Ksat,ver values for each soil horizon in Table 2 were multiplied by a scaling factor of 10, similarly as performed in the sensitivity analysis. The fh value was then automatically updated since this parameter represents the relation between horizontal and vertical conductivities in each cell. The dimensions of the cross-sections were defined as shown in Table 1 and consisted of increasing river depth and interactions between the porous media and the river network. Thus, by modifying these parameters, the calibration process mainly focused on increasing baseflow in the Ulla River watershed.    The R 2 values varied from 0.56 to 0.75 during the calibration period and from 0.76 to 0.85 during the validation period ( Table 7), showing that the model could explain most of the variability of the measured data in all hydrometric stations. The errors of the estimates were relatively small as shown by the low RSR values obtained at the different stations during both simulation periods (RSR ≤ 0.67). On the other hand, the PBIAS values showed some underestimation of measured values at the Sar hydrometric station (PBIAS ≤ 16.09%), and some overestimation of measured data for the remaining stations (−18.54% ≤ PBIAS ≤ −4.35%). Finally, the NSE values were relatively high in most locations, ranging from 0.55 to 0.84 during both periods, indicating that the residual variance was much smaller than the measured data variance. These indicators are comparable with the best flow estimates under natural regimes in Canuto et al. [27], with model predictions for the Ulla River watershed being considered extremely satisfactory if Moriasi et al. [63] guidelines are considered.

Prediction of River Flow in the Ulla River Watershed
Despite the good statistical results, Figure 9 further showed the difficulty of the MOHID-Land model in predicting the highest flow peaks in the Ulla River watershed. This was attributed to the precipitation inputs provided by the ERA5 dataset, which also failed to represent higher precipitation values as shown in Figure 3 for the Santiago meteorological station. This tendency was also verified by Hénin et al. [66], who concluded that reanalysis products such as ERA5 present an underestimation of heavy precipitation events. Hence, the amount of water entering the basin during heavy rain events in that dataset was clearly below measured rainfall and proved to be insufficient for the model to reach higher flow peaks. However, the replacement of the ERA5 data by measured values from the different weather stations was not considered as a viable solution to overcome this problem. Rainfall measurements were only available at a daily time step while the ERA5 data have an hourly time step. In the MOHID-Land model, daily measured data are distributed evenly during the day, reducing rainfall intensity rates, which would further reduce or even miss flow peaks. Besides that, one should consider that the model was calibrated for the entire period without special distinction between dry and wet seasons, which could also explain the difficulty for the model to reach peak flows in the Ulla River watershed. One way to overcome the problem from the model point of view would have been to define distinct geometries of the cross-sections per subbasin, and not per drainage area. This would increase model accuracy of river flow predictions at the watershed scale.

Conclusions
The MOHID-Land model is a complex, physically based, three-dimensional model used for catchment scale applications. The sensitivity analysis helped to identify the most relevant parameters/processes influencing river flow generation and for accurately modeling baseflow and peak flows. For the application in the Ulla River catchment, the resolution of the simulation grid, the choice of the infiltration method, and the evapotranspiration process were the main factors influencing river flow generation. The soil hydraulic properties, the depth of the soil profile, and the dimensions of the river cross-sections, which basically control the interactions between the porous media and the river, influenced baseflow. On the other hand, peak flows were mostly constrained by the channel's Manning coefficient, as well as the dimensions of the river cross-sections.
The sensitivity analysis further showed that the use of a too coarse resolution grid as well as the deactivation of the porous media and vegetation processes can compromise the quality of results, which should then be subjected to careful revision. Also, model simulations of soil infiltration considering the Darcian and the Green and Ampt approaches produced very similar outputs, both in terms of river flow values and computational time, meaning that users can choose between the two solutions depending on available data. Finally, the number of layers in the vertical simulation grid can lead to a substantial increase in the time needed to compute model simulations with no effect on river flow predictions. It is also important to note that the sensitivity analysis focused on just a few input parameters/processes used in MOHID-Land for simulating river flow at the watershed scale. Others should also be considered in future analysis, namely, the remaining soil hydraulic parameters, the crop coefficients, as well as the parameters used for water quality modeling.
Nevertheless, the MOHID-Land model is a powerful tool for simulating river flow at a daily scale in areas under natural flow regimes. This was demonstrated in simulations of the Ulla River flow at four locations, with comparisons between model predictions and measured values returning R 2 ≥ 0.56, RSR ≤ 0.67, and NSE ≥ 0.55. These same simulations showed a clear underestimation of river flow peaks, which was attributed to the quality of the ERA5 dataset and the misrepresentation of higher rainfall events.