Modeling Streamﬂow at the Iberian Peninsula Scale Using MOHID-Land: Challenges from a Coarse Scale Approach

: Hydrological modeling is nowadays critical for evaluating the status, past trends, and future perspectives of water availability at the global, regional, and local scales. The Iberian Peninsula is registering more frequent and severe droughts and water scarcity caused not only by extreme meteorological events, but also by increased demand for water for urban, industrial, and agricultural supplies. Better simulation models are thus needed for accurately quantifying the availability of local water resources. In this study, the natural ﬂow regime in different watersheds of the Iberian Peninsula was simulated using the process-based, fully distributed, MOHID-Land model from 1979 to 2013. Streamﬂow results were compared with measurements at 73 hydrometric stations not inﬂuenced by reservoirs, and with the data available in the management plans of each hydrographic region. The results showed a high dispersion of the goodness-of-ﬁt indicators, with the coefﬁcient of determination (R 2 ) ranging between 0 and 0.91, and the modeling efﬁciency (NSE) being higher than 0.35 at only 22 (calibration) and 28 (validation) hydrometric stations. Considering the scale of application, results were acceptable but evidenced the difﬁculties in simulating streamﬂow in watersheds using a coarse resolution. As such, this paper further deals with the difﬁculties and challenges of the adopted modeling approach. Nevertheless, this study constitutes a further step towards the more accurate assessment of water resources availability at the Iberian Peninsula scale using process-based modeling.


Introduction
According to the Intergovernmental Panel on Climate Change (IPCC) report [1], climate change is predicted to continue impacting the hydrology of river basins in Europe, particularly hydroclimatic extreme events. There is no generalized consensus about the areas that will be affected by the increase in the 100-year return period of discharge. Some studies project an increase in Continental Europe while others expect a decrease in parts of Northern and Southern Europe by 2100 [2,3]. Studies more focused on individual catchments indicate an increase in extreme discharge in regions such as Finland [4], France [5,6], and the Rhine basin [7,8]. On the other hand, meteorological droughts are to occur with greater intensity and for longer periods in Southern Europe. It is also highlighted that even in regions where an increase in summer precipitation is predicted, soil moisture may become more limited and hydrological droughts more severe because of the increasing evapotranspiration [9]. More recently, the updated IPCC report [10] projected an increase in air temperature over the Mediterranean from 1.5 • C to 2 • C compared to the present [11]. This temperature rise is intimately related to the increase in the drought risk in that area and the intensification of the hydrological cycle, facilitating more vapor in the atmosphere, activities. As such, by extending the work in Canuto et al. [24] to other transnational Iberian catchments, this study aims now to simulate the natural flow regime of the rivers in different watersheds of the Iberian Peninsula from 1979 to 2014. The estimation of the natural flow regime allows a better understanding and assessment of the impacts caused by changes made in the watersheds. On the other hand, the non-consideration of some processes, such as reservoirs management and water treatment plants discharge, helps reduce the sources of uncertainty related to the processes themselves and the input variables. Thus, in this application, reservoirs were not considered. Streamflow predictions were calibrated/validated through the comparison with measurements taken in hydrometric stations not influenced by this type of structure. Additionally, data available in the management plans of each hydrographic region was also used for model calibration/validation purposes. Hence, this paper aims: (i) to simulate streamflow in the Iberian transnational catchments using the MOHID-Land model; (ii) to quantify the prediction errors associated with the coarse-scale approach adopted in previous studies [21,22,24]; (iii) to discuss the limitations of the modeling approach and the challenges that still need to be faced for improving the hydrological model when implemented at such a coarse scale.

Description of the Studied Area
The Iberian Peninsula, in southwest Europe, covers a total area of 583,832 km 2 , mainly occupied by Portugal (89,060 km 2 ) and Spain (492,175 km 2 ) [25,26]. The Peninsula is divided into 24 river basin districts, with these being defined as "the area of land and sea, made up of one or more neighboring river basins together with their associated groundwaters and coastal waters" and constituting the main units for management of river basins according to the EU directive Establishing a framework for Community action in the field of water policy [27].
In this application, the study area covered a total area of 438,178 km 2 (75% of the Iberian Peninsula), comprehending 19 river basin districts located in the west and north part of the Iberian Peninsula (Figure 1a). The main river basins were Douro (98,107 km 2 ), Tagus (86,283 km 2 ), Guadiana (78,234 km 2 ), and Guadalquivir (57,228 km 2 ), with the first three being transnational basins divided between Portugal and Spain (Table 1). According to the Koppen climate characterization, the selected region is mainly characterized by four climate units (Figure 1b). In the north, the climate is temperate, with warm and dry summers in the northwest (Csb), and warm but not dry summers in the northeast (Cfb). The south and southwest regions have a temperate climate characterized by dry and hot summers (Csa), while the eastern part is characterized by an arid, steppe, and cold climate (BSk). The elevation of the study area ranges between 5 m and 3332 m (Figure 1a), with the highest value located in Sierra Nevada, Spain. In Portugal, the highest point is at an altitude of about 1990 m in Serra da Estrela. The more representative soil reference groups are Cambisols (43%), Regosols (22%), Leptosols (11%), Luvisols (10%), and Fluvisols (5%) [28]. The land use is mainly characterized by non-irrigated arable land (18.7%), broad-leaved forest (11.1%), agro-forestry areas (7.4%), sclerophyllous vegetation (6.9%), and transitional woodland-shrub (6.9%) [29]. Table 1. Characteristics of the studied river basin districts (Area; number of reservoirs and respective total capacity, C total ; total storage capacity, C storage ; natural surface flow; regularization index, RI). The studied domain further comprises a total of 1113 reservoirs, 71% of these with a total capacity lower than 10 hm 3 , 28% between 10 and 1000 hm 3 , and 1% higher than 1000 hm 3 [40,41]. As presented in Table 1 3 ) instead. On the other hand, the impact of reservoirs does not only depend on the storage capacity installed in a river basin district, but it is intimately related with the ratio between the surface flow and the storage capacity, which is here represented by the regularization index (RI) computed as follows:

River
The surface flow values presented in Table 1 correspond to the natural flow regime in each basin district. These values were obtained from the Portuguese and Spanish river The studied domain further comprises a total of 1113 reservoirs, 71% of these with a total capacity lower than 10 hm 3 , 28% between 10 and 1000 hm 3 , and 1% higher than 1000 hm 3 [40,41]. As presented in Table 1 3 ) instead. On the other hand, the impact of reservoirs does not only depend on the storage capacity installed in a river basin district, but it is intimately related with the ratio between the surface flow and the storage capacity, which is here represented by the regularization index (RI) computed as follows: The surface flow values presented in Table 1 correspond to the natural flow regime in each basin district. These values were obtained from the Portuguese and Spanish river basin management plans. River basin management plans are the reports developed under the implementation of the Water Framework Directive of the European Union [27] and aim to characterize the river basin district, diagnose the state of its water bodies in terms of water quantity and quality, and propose solutions to improve that state when necessary. In these sources, for both countries, the natural flows were estimated using the Témez model [42] for the Portuguese river basin districts and the SIMulación Precipitación-Aportación (SIMPA) model [43][44][45] for the Spanish territory. This latter approach was developed by CEDEX and is also based on the Témez model. Considering the regularization index, the river basin districts where the storage capacity has more impact on streamflow are the Portuguese Guadiana (266%), Spanish Guadiana (248%), Guadalete and Barbate (151%), Spanish Tagus (132%), and Guadalquivir (117%), with all having a storage capacity higher than the volume of water naturally produced in an average year. The alterations in the natural flow regime of these rivers caused by the presence of these dams were already reported in different studies [18,46,47]. Aside from the modifications in the flow regimes caused by dams, the transfer of water between Tagus and Segura watersheds (not considered in the modeling approach) also affects the natural regime flow in both watersheds [48], modifying the hydrological cycle in Tagus headwaters from where water is transferred [49].

Model Description
MOHID-Land [23,24] is a fully distributed, physically based model, which uses mass and momentum conservation equations to simulate water movement between four main compartments (atmosphere, porous media, soil surface, and river network) using a finite volume approach. To avoid instability problems and save computational time, the model time step is variable, acquiring higher values during dry seasons, and lower values in wet periods when water fluxes increase.
The atmosphere processes are not explicitly simulated, but this compartment allows the input of atmospheric data, which can be space and time-variant, as surface boundary conditions. The simulated domain is represented by a regular grid in the surface plane, with a user-defined grid size and cell resolution. The soil discretization (three-dimensional) considers this same grid in the horizontal plane and a division by layers with variable thickness in the vertical direction following a cartesian coordinate system. The river network is defined by a one-dimensional (1D) domain, by connecting the surface cell centers (nodes) with lower altitude in the digital terrain model.
The free surface flow is computed according to the Saint-Venant equation in its conservative form and considering the advection, pressure, and friction forces: where Q is the water flow (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), x i represents the xyz directions (−) and the subscripts u and v denote flow directions. This equation is solved for one direction (1D domain) in the river network and for two directions (2D domain), coinciding with the directions of the horizontal grid, in cells without drainage network (overland flow). In the boundary of the drainage network and overland flow, the water exchanges are calculated according to a kinematic approach, neglecting bottom friction.
In the porous media, the variable-saturated water flow is computed using the Richards' equation in three directions (3D domain): 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 corresponds to the sink term, representing the root water uptake (L 3 L −3 T −1 ). The soil hydraulic properties are defined following the van Genuchten Mualem functional relationships [50,51]. The water exchanges between the porous media and the river network are driven by the pressure differences in the interface of these mediums. The crop evapotranspiration rates (ET c ) are obtained from the product of the reference evapotranspiration (ET o ) rates computed according to the FAO Penman-Monteith method and a crop coefficient (K c ) [52]. ET c rates are then divided into potential soil evaporation (E p , LT −1 ) and crop transpiration (T p , LT −1 ) based on Ritchie et al. [53]. The T p rates define the maximum values of the sink term for root water uptake in the Richards equation. These may be reduced due to the presence of rootzone stressors following a macroscopic approach proposed by Feddes et al. [54]. On the other hand, the actual soil evaporation (E a , LT −1 ) is estimated by imposing a pressure threshold value to the potential evaporation values [55].
The connection between the surface runoff and the porous media is based on infiltration and exfiltration processes [56]. The infiltration rate can be simulated according to three approaches, namely, Darcy's law, the Green and Ampt method [57], or a modified version of the SCS curve number (SCS-CN) method [58]. In the SCS-CN case, the first step comprehends the estimation of the surface runoff, with the remaining amount of water infiltrating the soil; however, when the soil has saturated, a correction in the infiltration rate (and consequently in the runoff flow) is made to respect the physical limitations of the system.
More detailed information on the MOHID-Land model can be found in Canuto et al. [24], Oliveira et al. [59], and Ramos et al. [60].

Model Setup
The MOHID-Land model was implemented in the studied area using a constant horizontally spaced grid, with a resolution of 0.045 The river network was delineated based on the interpolated DTM, considering the steepest slope in all 8 directions of each cell, by connecting the cell centers (nodes) with the lowest height. The river cross-sections were defined according to the Strahler order of each reach [62,63]. A trapezoidal geometry was adopted for all cross-sections, with the height and top and bottom widths being set according to Canuto et al. [24].
The soil data were interpolated from the rasterized Harmonized World Soil Database [28], with an original resolution of 30 arc-second (~1 km). The main soil units identified in the model domain were Calcic Cambisol (28%), Calcaric Regosol (17%), Humic Cambisol (12%), Lithosol (7%), Calcaric Fluvisol (5%), and Albic Luvisol (5%) (Figure 2d). In-depth, the soil was defined by three horizons comprehending six grid layers. The surface horizon corresponded to a layer of 0.3 m thickness, the middle horizon was composed of a layer with 0.5 m thickness, and the bottom horizon comprehended four layers with depth increasing thicknesses of 0.42, 0.50, 10.0, and 10.0 m. For each soil horizon, the Mualem-van Genuchten model parameters [50] were obtained from the soil texture classes described in the Harmonized Soil Database using the HYPRES class pedotransfer functions [64]. The soil's maximum and minimum depth were defined as 30 and 0.3 m, respectively, with the soil depth in each cell being adjusted according to its slope. The initial soil condition was set as 80% saturated (from the bottom to the surface) and in the unsaturated zone water content was set to field capacity.  The land use in the studied domain was obtained from the CORINE Land Cover (CLC) 2012, with a resolution of 100 m [29]. The information of the CLC map was interpolated to the horizontal grid by assigning to each CLC class: (i) a Manning coefficient according to van der Sande et al. [65]; (ii) a corresponding vegetation class from the MOHID-Land database ( Figure 2e); an annual K c value following Canuto et al. [24] (Figure 2f). The hydrologic soil groups were defined by combining land use and soil texture data to derive the respective CN values following the Soil Conservation Service [58] tables ( Figure 2c).
The meteorological data were obtained from the SAFRAN model, which was developed, calibrated, and validated by Quintana-Seguí et al. [66]. SAFRAN is a meteorological analysis system based on an optimal algorithm that combines observations and a first guess, such as the outputs of a global meteorological model. The resulting dataset is in gridded format, with a resolution of 5 km, and includes hourly data of precipitation, surface air temperature, wind speed, relative humidity, and solar radiation from 1979 to 2014. The information of these meteorological properties was interpolated to the MOHID-Land grid following a triangular method.

Model Evaluation
Model performance was evaluated by comparing simulated and measured monthly streamflow data at 73 stations not influenced by reservoirs operation (Figure 2a), i.e., only hydrometric stations under natural flow regime were considered. The simulation period was from 1985 to 2013 (35 years). The measured dataset was obtained from MITERD [67] and SNIRH [68] for the stations in Spain and Portugal, respectively. Table A1, in Appendix A, presents the number of records as well as the minimum, maximum, mean, and standard deviation values of the measured monthly streamflow in each station. The calibration period was defined from January 1985 to December 1999, while the validation period went from January 2000 to December 2013. A warm-up period was also considered from September 1979 to December 1984. The calibration procedure considered the modification of selected model parameters, one at a time, within reasonable ranges, to minimize deviations between simulated and observed streamflow in the selected hydrometric stations. Following Oliveira et al. [59], who identified the most sensitive parameters affecting simulations of streamflow in MOHID-Land, the modified parameters were: the vertical saturated hydraulic conductivity; the multiplying factor relating the vertical and horizontal saturated hydraulic conductivities (f h ); the surface and channel Manning coefficients; the crop coefficients; the dimensions (height and top and bottom widths) of the river cross-sections. The calibrated parameters were then used for validation of streamflow predictions, with the deviation between model simulations and observed data being assessed again at the same hydrometric stations using the validation dataset.
Model performance was evaluated using four statistical parameters, namely, the coefficient of determination (R 2 ), the percent bias (PBIAS), the root mean square errorobservation standard deviation ratio (RSR), and the Nash-Sutcliffe model efficiency (NSE), which are computed, respectively, as follows: where Q i obs and Q i sim are the observed and simulated flow on day I, respectively, Q mean obs and Q mean sim are the observed and simulated mean flow for the analyzed period, respectively, and p is the total number of days in this period. When NSE > 0.5, RSR ≤ 0.7, PBIAS ± 25%, and R 2 > 0.5 the results of the modelled streamflow are considered satisfactory [69].
Finally, model results were also compared with the official values available in Portuguese river basin management plans and in a report produced by Centro de Estudios y Experimentación de Obras Públicas (CEDEX) in the scope of Spanish river basin management plans (Table 1). This type of validation allows the assessment of the general values of streamflow obtained by watershed instead of considering only the headwater areas defined by the hydrometric stations presented before; however, it is important to denote that the values presented in the river basin management plans are also results of models and, consequently, they are also subjected to uncertainty, but in the end, they constitute official information. The values available in the river basin management plans and used in this study for comparison with model results were the average annual precipitation and the average annual natural surface runoff in each river basin district. The corresponding difference, in percentage, between the MOHID-Land modeled and official results were computed as follows: where V model and V official are the modeled and official average annual values of precipitation/surface runoff, respectively, and ∆ is the difference between these values.

Model Parametrization
The best fit between simulated and observed streamflow values was obtained for a simulation in which the Manning coefficient ranged in the different soil units from 0.02 to 0.298 s m −1/3 (Figure 2b), and the annual K c value for the different crops varied from 0.15 to 1.02 (Figure 2f). The curve number was set to values between 68 and 80 ( Figure 2d), with these values not having a significant impact on the streamflow simulation, as demonstrated by Oliveira et al. [59]. The multiplying factor f h was set to 10. The soil hydraulic parameters were adjusted to the values presented in Table A2 of Appendix B, with the saturated water content (θ s ), the residual water content (θ r ), the empirical shape parameters α and η, and the vertical saturated hydraulic conductivity (K s ) ranging between 0.35 and 0.766 m 3 m −3 , 0.01 and 0.025 m 3 m −3 , 0.198 and 3.83 m −1 , 1.086 and 1.377, and 9.26 × 10 −7 and 6.94 × 10 −6 m s −1 , respectively. The connectivity/tortuosity parameter (L, m) adopted the value 0.5 as proposed by Mualem [48]. A summary of the ranges for the calibrated parameters is presented in Table 2. Finally, the calibrated dimensions for the river cross-sections were defined as presented in Table 3. The calibrated values of the Manning coefficient, K c , CN, f h , and the dimensions of the river cross-sections were thus in agreement with those used in Canuto et al. [24] and Oliveira et al. [59] for simulating streamflow in different catchments of the Iberian Peninsula. The soil hydraulic parameters were also in the range of those proposed by Ramos et al. [70] for the different texture classes of soils in Portugal.   Figure 3 presents the histograms of the distribution of the statistical parameters obtained by comparing simulated and measured streamflow data at 73 hydrometric stations selected for calibration/validation. In general, the statistical parameters showed a large dispersion. In the calibration period, the R 2 , RSR, PBIAS, and NSE values ranged from 0 to 0.91, 0.36 to 9.83, −1160% to 81%, and −96 to 0.87, respectively. The R 2 histogram showed that the most populated classes were those between 0.4 and 0.6 (at 21 stations) and 0.6 and 0.8 (21 stations), thus expressing the model's capability in explaining most of the variability of the observed data at 57% of the selected hydrometric stations. For the RSR parameter, values higher than 1 were obtained at 29 stations, being the most frequent and indicating a large error of the estimate. Yet, at 31 stations the RSR values were within the 0.4 to 1.0 range, thus showing a more acceptable performance. The class with a deviation lower than −25% was the most populated (33 stations) for the PBIAS statistical parameter, revealing a small tendency of overestimation of the observed data in these stations. Lastly, the NSE values lower than 0.15 were the most frequent at 32 stations; however, the NSE values were within the acceptable range at 22 stations, being higher than 0.35 in 11 of those stations and higher than 0.55 in the remaining 11 stations. In the validation period, the R 2 values ranged between 0.01 and 0.88, the RSR between 0.4 and 19.22, the PBIAS was in the range −3500% and 90%, and the NSE varied from −368 to 0.84. In this period, the most populated classes for R 2 , RSR, PBIAS, and NSE were, respectively, ]0. 6 values were within the acceptable range at 22 stations, being higher than 0.35 in 11 of those stations and higher than 0.55 in the remaining 11 stations. In the validation period, the R 2 values ranged between 0.01 and 0.88, the RSR between 0.4 and 19.22, the PBIAS was in the range −3500% and 90%, and the NSE varied from −368 to 0.84. In this period, the most populated classes for R 2 , RSR, PBIAS, and NSE were, respectively, ]0.6,0    [24] for the Guadiana basin, showing the difficulties in correctly portraying the spatial dynamics of evapotranspiration, soil moisture, and land-use change impacts on simulations of streamflow in coarse-scale applications of a distributed-based model such as MOHID-Land in highly anthropogenic watersheds. As such, the use of these models for conducting scenario analysis should be carried out with great care due to the great uncertainty associated with model predictions, and calibration of model parameters should definitely not be limited to one or two locations as sometimes is observed in the literature since these are hardly representative of the hydrological processes across the catchments.  A more general analysis allowed the comparison of the annual average precipitation from the SAFRAN model, and the respective annual natural surface flow generated by the MOHID-Land model with the corresponding information available in the Portuguese and Spanish river basin management plans (Table 4). With the natural surface flow values available in these sources already presented (Table 1) and explained in the description of the studied area, it is also important to denote that the average annual precipitation presented in river basin management plans was estimated based on observations from meteorological stations. The official precipitation and flow annual average values were estimated considering a period of 28 years for Portugal (1985-2013) and 38 years for Spain (1980-2018).  [38] The comparison between precipitation values used as input to MOHID-Land and the river basin management plans showed that the main differences occurred in the Portuguese river basin districts, namely, Vouga, Mondego, and Lis (∆ = 39%); Cávado, Ave, and Leça (∆ = 27%); Minho and Lima (∆ = 17%); Douro (PT), Tagus and West Rivers (PT), and Algarve Rivers (∆ = 15%) (Figure 5a). Precipitation data used in the MOHID-Land model were thus higher than that considered in the river management plans; however, the pattern found for the precipitation differences was not reproduced for flow, with the more significative flow differences occurring in Spanish river basin districts, especially in Guadiana (ES) (∆ = 169%), Guadalate and Barbate (∆ = −122%), and Guadalquivir (∆ = −121%) (Figure 5b). This time, flow data in the MOHID-Land model were smaller than in the river management plans at these watersheds.

Model Performance
Although SAFRAN was classified as a robust dataset, differences in precipitation values could result from limitations of the SAFRAN algorithm, already highlighted by MITERD [67], namely, the overestimation of the number of precipitation days, the missing of high precipitation events, and the errors often found at the border of the homogenous climatic zones; however, the high precipitation differences reported in Portugal result from the fact that SAFRAN has been extended to this country without considering any meteorological station in this area to better apply the SAFRAN methodology, which favors high station density, thus resulting in the reported difference errors. Although SAFRAN was classified as a robust dataset, differences in precipitation values could result from limitations of the SAFRAN algorithm, already highlighted by MI-TERD [67], namely, the overestimation of the number of precipitation days, the missing of high precipitation events, and the errors often found at the border of the homogenous climatic zones; however, the high precipitation differences reported in Portugal result from the fact that SAFRAN has been extended to this country without considering any meteorological station in this area to better apply the SAFRAN methodology, which favors high station density, thus resulting in the reported difference errors.
The differences found for natural surface flow do not have the same magnitude as the precipitation differences, with the former being much higher than the latter in central and south Spanish river basin districts where the differences in precipitation are not significant. Since the only water input in this implementation was rainfall, these results show that the hydrological model might be evaporating or retaining too much of the precipitated water in the soil or sub-surface layers instead of accounting it to surface flow. In Portugal, Minho and Lima, Douro (PT), Tagus and West Rivers (PT), Sado and Mira, and Guadiana (PT) also presented a similar behavior, but less pronounced. Finally, model results in some river basin districts, such as Eastern Cantabrian and Andalusia Mediterranean Basins, correctly reproduced the average annual natural surface flow or, at least, maintained a correspondence between surface flow and precipitation when the latter was overestimated by SAFRAN. This occurred in the river basin districts where the drainage network had more and/or shorter branches (Figure 1a).

About Model Predictions
Canuto et al. [24] extensively discussed the challenges of modeling streamflow in highly complex, anthropogenic catchments using a coarse-scale approach and the MOHID-Land model. These referred to the calibration procedures to meet the set of parameters that best describe the simulated landscape processes, but also a set of causes related to field measurements, and to model input and model structure errors. While most of those causes are still valid for this study, a few more can be added.
The choice of a 5 km grid resolution hindered the correct delineation of the drainage network in the headwaters of the watershed where the analyzed hydrometric stations were located. That resolution grid was chosen to save computational time since this factor has a significant impact on simulation time [59]. Errors in the delineation of the river network happened even after the main river lines were forced to flow through the correct cells of the grid domain by using the burn-in process in the MOHID-Land model. This The differences found for natural surface flow do not have the same magnitude as the precipitation differences, with the former being much higher than the latter in central and south Spanish river basin districts where the differences in precipitation are not significant. Since the only water input in this implementation was rainfall, these results show that the hydrological model might be evaporating or retaining too much of the precipitated water in the soil or sub-surface layers instead of accounting it to surface flow. In Portugal, Minho and Lima, Douro (PT), Tagus and West Rivers (PT), Sado and Mira, and Guadiana (PT) also presented a similar behavior, but less pronounced. Finally, model results in some river basin districts, such as Eastern Cantabrian and Andalusia Mediterranean Basins, correctly reproduced the average annual natural surface flow or, at least, maintained a correspondence between surface flow and precipitation when the latter was overestimated by SAFRAN. This occurred in the river basin districts where the drainage network had more and/or shorter branches (Figure 1a).

About Model Predictions
Canuto et al. [24] extensively discussed the challenges of modeling streamflow in highly complex, anthropogenic catchments using a coarse-scale approach and the MOHID-Land model. These referred to the calibration procedures to meet the set of parameters that best describe the simulated landscape processes, but also a set of causes related to field measurements, and to model input and model structure errors. While most of those causes are still valid for this study, a few more can be added.
The choice of a 5 km grid resolution hindered the correct delineation of the drainage network in the headwaters of the watershed where the analyzed hydrometric stations were located. That resolution grid was chosen to save computational time since this factor has a significant impact on simulation time [59]. Errors in the delineation of the river network happened even after the main river lines were forced to flow through the correct cells of the grid domain by using the burn-in process in the MOHID-Land model. This tool allows the subtraction of a constant value to the elevation in cells identified as river path (usually, this information is available in shapefile format). Since the drainage network is delineated to make the river lines follow the lower elevation cells, this process should force those lines crossing the right cells, correcting the river path; however, after observing the details of the drainage network in the watershed's headwaters, it is possible to identify several cases where the adoption of this coarse resolution caused a connection between cells that corresponded to different river lines. As exemplified in Figure 6, this situation led to significant differences between the real and the simulated drained area of a hydrometric station, with an obvious impact on the amount of water flowing through that sub-basin.
should force those lines crossing the right cells, correcting the river path; however, after observing the details of the drainage network in the watershed's headwaters, it is possible to identify several cases where the adoption of this coarse resolution caused a connection between cells that corresponded to different river lines. As exemplified in Figure 6, this situation led to significant differences between the real and the simulated drained area of a hydrometric station, with an obvious impact on the amount of water flowing through that sub-basin. Oliveira et al. [59] tested the MOHID-Land model with two different grid resolutions (0.5 and 1 km) in the Ulla basin, Northern Spain, and found that the coarser application would cause a decrease of more than 70% in streamflow values as a result of a less detailed representation of the watershed. These results suggest that, besides the errors previously mentioned in the delineation of the drainage network, the grid resolution in this type of model can lead to inferior streamflow generation by itself, especially when the adopted resolution cannot simulate the processes in the watershed with sufficient detail. The same had also been experienced in Canuto et al. [24]. These conclusions are also corroborated by Refsgaard [71], who applied the MIKE SHE model in the Karup basin (440 km 2 ), Denmark, with four different resolutions, namely, 0.5, 1, 2, and 4 km, and concluded that there was a maximum grid size for a good model performance of streamflow simulations and result deteriorated for resolutions lower than this maximum. Sreedevi and Eldho [72] applied the SHETRAN model in the humid tropical Netravathi basin (3657 km 2 ), Southern India, reporting that for three used resolutions (4, 2, and 1 km), there was no gain in the streamflow performance.
Furthermore, Oliveira et al. [59] demonstrated that cross-sections geometry could significantly impact MOHID-Land simulations of streamflow. In that study, the differences in the highest and lowest values (extremes) of streamflow showed an increase of 39% and 30%, respectively, after increasing cross-sections height by 100% when compared to the baseline simulation. Canuto et al. [24] also modified their cross sections geometry in the calibration process, reporting a ratio between the calibrated and default cross-sections area ranging from 1.5 to 10.5, thus demonstrating the sensitivity of streamflow simulation in MOHID-Land to channel geometry. Since the domain covered in the present case study is extensive, the variation of model inputs can be very large, including the Oliveira et al. [59] tested the MOHID-Land model with two different grid resolutions (0.5 and 1 km) in the Ulla basin, Northern Spain, and found that the coarser application would cause a decrease of more than 70% in streamflow values as a result of a less detailed representation of the watershed. These results suggest that, besides the errors previously mentioned in the delineation of the drainage network, the grid resolution in this type of model can lead to inferior streamflow generation by itself, especially when the adopted resolution cannot simulate the processes in the watershed with sufficient detail. The same had also been experienced in Canuto et al. [24]. These conclusions are also corroborated by Refsgaard [71], who applied the MIKE SHE model in the Karup basin (440 km 2 ), Denmark, with four different resolutions, namely, 0.5, 1, 2, and 4 km, and concluded that there was a maximum grid size for a good model performance of streamflow simulations and result deteriorated for resolutions lower than this maximum. Sreedevi and Eldho [72] applied the SHETRAN model in the humid tropical Netravathi basin (3657 km 2 ), Southern India, reporting that for three used resolutions (4, 2, and 1 km), there was no gain in the streamflow performance.
Furthermore, Oliveira et al. [59] demonstrated that cross-sections geometry could significantly impact MOHID-Land simulations of streamflow. In that study, the differences in the highest and lowest values (extremes) of streamflow showed an increase of 39% and 30%, respectively, after increasing cross-sections height by 100% when compared to the baseline simulation. Canuto et al. [24] also modified their cross sections geometry in the calibration process, reporting a ratio between the calibrated and default cross-sections area ranging from 1.5 to 10.5, thus demonstrating the sensitivity of streamflow simulation in MOHID-Land to channel geometry. Since the domain covered in the present case study is extensive, the variation of model inputs can be very large, including the channel crosssections geometry. The database published by Andreadis et al. [73], who derived the river widths and depths using HydroSHEDS river topology dataset and simple geomorphic relationships between those dimensions and the drained area and the discharge, allows understanding the variation of channel geometry in the domain. Figure 7 shows the ratio between the estimated width and depth (a) and the drained area (b) using that database. The referred figure clearly shows that the ratio between width/depth does not vary linearly with the drained area and that for similar drained areas, this ratio is higher in the north than in the south; however, these differences were not able to be implemented in the model as the available preprocessing tool only allows the assignment of cross-section geometries to the entire modeled domain. Hence, it was not possible to represent with detail the existing differences in the channel river cross-sections geometry throughout the domain, preventing a detailed calibration/validation process and, consequently, deteriorating results. For this reason, applications of the MOHID-Land model at the catchment scale instead of at the regional scale are highly recommended.
shows the ratio between the estimated width and depth (a) and the drained area (b) using that database. The referred figure clearly shows that the ratio between width/depth does not vary linearly with the drained area and that for similar drained areas, this ratio is higher in the north than in the south; however, these differences were not able to be implemented in the model as the available preprocessing tool only allows the assignment of cross-section geometries to the entire modeled domain. Hence, it was not possible to represent with detail the existing differences in the channel river cross-sections geometry throughout the domain, preventing a detailed calibration/validation process and, consequently, deteriorating results. For this reason, applications of the MOHID-Land model at the catchment scale instead of at the regional scale are highly recommended.

Figure 7.
Ratio between cross-sections width and depth (a) and the drained area (b) in the studied domain, according to [72].
Finally, it is important to denote that the solution presented here is the best from a set of tests performed during the calibration process; however, since each simulation takes about 15 days to perform, it is possible that with more efficient implementation, the calibration process would have gone further, and better results would be achieved. This also demonstrates that in fully distributed physically based models such as MOHID-Land the calibration process should not be dismissed, especially when the modeled domain is characterized by a significant lack of data [74]. The importance of the calibration process in MOHID-Land model applications is demonstrated in Oliveira et al. [59], Canuto et al. [24], Epelde et al. [75], and Brito et al. [21], where the first three studied domains were extensively calibrated, and satisfactory results for streamflow simulation were obtained, while the latter never referred the existence of a calibration process which along with the choice of hydrometric stations influenced by reservoirs to evaluate the model performance appears to be the reason for the model not being able to represent the observed values.

Conclusions
The MOHID-Land model was used to simulate the natural flow regime in Iberian transnational catchments using a process-based, fully distributed approach implemented in a coarse (5 km) resolution grid domain. While results were satisfactory, they also evidenced the difficulties in simulating the streamflow in highly modified basins following such a complex modeling approach. As such, the goodness-of-fit indicators computed at Finally, it is important to denote that the solution presented here is the best from a set of tests performed during the calibration process; however, since each simulation takes about 15 days to perform, it is possible that with more efficient implementation, the calibration process would have gone further, and better results would be achieved. This also demonstrates that in fully distributed physically based models such as MOHID-Land the calibration process should not be dismissed, especially when the modeled domain is characterized by a significant lack of data [74]. The importance of the calibration process in MOHID-Land model applications is demonstrated in Oliveira et al. [59], Canuto et al. [24], Epelde et al. [75], and Brito et al. [21], where the first three studied domains were extensively calibrated, and satisfactory results for streamflow simulation were obtained, while the latter never referred the existence of a calibration process which along with the choice of hydrometric stations influenced by reservoirs to evaluate the model performance appears to be the reason for the model not being able to represent the observed values.

Conclusions
The MOHID-Land model was used to simulate the natural flow regime in Iberian transnational catchments using a process-based, fully distributed approach implemented in a coarse (5 km) resolution grid domain. While results were satisfactory, they also evidenced the difficulties in simulating the streamflow in highly modified basins following such a complex modeling approach. As such, the goodness-of-fit indicators computed at 73 hydrometric stations located in basins headwaters showed large variability, with the Nash-Sutcliffe modeling efficiency (NSE) reporting acceptable results for only 22 and 28 of those locations during calibration and validation, respectively. Nevertheless, the MOHID-Land application presented here shows the capacity of this model to simulate water processes at the regional scale instead of considering the division and simulation by watershed. This work also highlights the difficulties in calibrating the model in such a vast domain, as well as the importance of adopting a resolution with sufficient detail for representation of the physical processes without losing computational efficiency. Computational efficiency is essential for a successful calibration process.
The implementation of the MOHID-Land model at a regional scale can contribute to better water governance of the Iberian water resources, namely those shared by Portugal and Spain. While the modeling approach still needs to be improved, the implementation of