Interpreting the Manning Roughness Coefficient in Overland Flow Simulations with Coupled Hydrological-Hydraulic Distributed Models

There is still little experience on the effect of the Manning roughness coefficient in coupled hydrological-hydraulic distributed models based on the solution of the Shallow Water Equations (SWE), where the Manning coefficient affects not only channel flow on the basin hydrographic network but also rainfall-runoff processes on the hillslopes. In this kind of model, roughness takes the role of the concentration time in classic conceptual or aggregated modelling methods, as is the case of the unit hydrograph method. Three different approaches were used to adjust the Manning roughness coefficient in order to fit the results with other methodologies or field observations— by comparing the resulting time of concentration with classic formulas, by comparing the runoff hydrographs obtained with aggregated models, and by comparing the runoff water volumes with observations. A wide dispersion of the roughness coefficients was observed to be generally much higher than the common values used in open channel flow hydraulics.


Introduction
The estimation of overland flow is one of the most relevant processes when assessing the hydrological response of a basin. Despite hydrological processes being extremely complex and, nowadays, still not fully understood, their evaluation is essential to characterize the basin flow production under rainfall events. The characterisation of overland flow is a crucial aspect for flood hazard management, which is the most catastrophic natural hazard and results in the greatest material and life damages all over the world [1][2][3]. The flow characteristics are key inputs to perform particular risk analyses as well as for integrating into flood risk early warning systems [4][5][6][7].
Numerical modelling is a widely extended technique that allows studying this kind of complex phenomena for which, in general, no analytical solutions exist. Hydrological modelling has been traditionally performed by means of the so-called aggregated models [7,8], which rely on the simplification of the study area as a single whole unit (generally the basin) or as the contribution (or "aggregation") of different parts (subbasins) constituting it all linked through a routing method. These models are usually based on conceptual or empirical formulations to describe several physical processes related to precipitation, infiltration, vegetation interception, overland storage, surface runoff, etc., and normally use a representation of the basin, or subbasins, as a homogenous unit [9].
In this aggregated approach, one of the fundamental parameters needed for the transformation of the rainfall into runoff is the time of concentration (t c ). At present, of first performing a hydrologic simulation with an aggregated model, and subsequently using a hydraulic model to characterize the flow. This coupled model technique usually implies solving the equations of mass and momentum conservation in two dimensions, i.e., the two-dimensional shallow water equations (2D-SWE). In these, the hydrological processes (rainfall-runoff transformation and losses) are usually considered as source terms of the mass conservation equation, and their contribution to momentum is considered only in a few cases.
Both aggregated and coupled hydrological-hydraulic distributed models depend on a large number of parameters (e.g., soil and land characteristics, vegetation, topography, atmospheric conditions, etc.) [34][35][36][37], but also depend on the practitioner's expertise. Nevertheless, coupled distributed models constitute a conceptual change with respect to aggregated models as they are based on a physical approach. This means that the routing of overland flow is dependent on the terrain roughness and is considered in the model with hydraulic equations, commonly the Manning equation. Thus, the parameters: time of concentration (t c ) and lag time (t lag ) are no longer part of the modelling process, with the Manning roughness coefficient being the only parameter apart from geometry that conditions the runoff propagation.
For hydraulic applications, the Manning roughness coefficient has been widely studied; numerous publications and guidelines for its estimation in natural or artificial channels exist [38][39][40]. However, in coupled hydrological-hydraulic distributed modelling, the value of roughness coefficients must incorporate not only the material roughness itself, but also all the sub-grid geometric features that affect friction, or energy dissipation, and thus the water depth and flow velocity [38]. It is worth mentioning that the flow patterns of overland flow in hydrological models are totally different than for standard uniform steady open channel flow, for which the Manning equation was originally introduced. In hydrological applications, the water depth values can be of a few millimetres or less. Thus, the roughness coefficient values for these extremely shallow overland flows characteristic of hydrological applications can, by no means, be those in the aforementioned guidelines, or at least, those must be questioned and verified. For overland flows, such as those on hillsides, the coefficient values tend to be higher than common ones used in river or stream flows [31,41,42]. Moreover, in coupled hydrological-hydraulic distributed models, the terrain roughness is a fundamental parameter that governs water propagation. It, together with the geometry discretisation, takes the role of determining the basin response as the time of concentration does in aggregated hydrological models.
This work aims to provide insight on the role of the Manning roughness coefficient in overland flow simulations by means of two-dimensional coupled hydrological-hydraulic distributed numerical models, and, thus, to assist in the construction of more accurate distributed rainfall-runoff models. For that purpose, three different strategies were used to characterize the hydrological response of a basin. The first one is a set of four small basins where the concentration time was estimated using empirical formulas and the Manning coefficient in the distributed model was adjusted for the S-shaped hydrograph to fit with the concentration time. In the second case, the resulting hydrograph obtained by means of traditional aggregated hydrological modelling in a basin was used to calibrate the terrain roughness in a 2D-SWE-based distributed model. Finally, the last strategy consisted of the calibration of the roughness coefficients of a monitored basin to adjust the computed results in four observed rainfall events.

Materials and Methods
As already introduced in the previous section, two different modelling approaches were used in the present work: (1) traditional aggregated hydrological modelling and (2) coupled hydrologic-hydraulic distributed modelling. The two techniques are briefly presented in the following sections.

Aggregated Hydrological Modelling
Case Study 2, which consists of a 1016 km 2 basin located in Mexico, was firstly analysed with an aggregated hydrological modelling approach through the HEC-HMS software [43]. A standard methodology was applied, based on the local administration guidelines, to discretize the basin. The unit hydrograph (UH) methodology was used for the rainfall-runoff transformation processes, which was obtained from the dimensionless unit hydrograph of the Soil Conservation Service [16,44]. The t lag parameter and the value of the time of concentration (t c ) were estimated by means of the Kirpich formula. Precipitation losses were accounted for using the SCS-CN method [44,45], also referred to as the NRCS-CN method [46] after the U.S. Soil Conservation Service was renamed as Natural Resources Conservation Service. All parameters of the numerical model are detailed in the Supplementary Material.
Iber solves the 2D Shallow Water Equations (2D-SWE) using a conservative scheme based on the Finite Volume Method (FVM) on a structured or/and unstructured mesh of triangles and/or quadrilaterals. It uses the Roe scheme [60], which consists of the Godunov method together with Roe Approximate Riemann Solver [33].
The shear stress terms due to friction are generally incorporated in the SWE via the concept of the friction slope. This term expresses the contribution of the momentum change due to the energy dissipation produced by flow-boundary interactions as well as by a sub-grid of obstructions and by the losses due to the flow turbulence if, as it is generally done in coupled hydrological-hydraulic distributed modelling, no turbulence model is used. The friction slope (S o ) is calculated using the Manning formula, which, in 2D, results in the following equations for the X and Y directions: where n is the Manning roughness coefficient; U x and U y are the two components of the depth averaged velocity vector in X and Y direction, respectively; |U| is the modulus of the flow velocity; and h is the depth. The DHD numerical scheme [25], developed ad hoc for hydrological modelling purposes, was used in all simulations with Iber. Briefly, this scheme merges the hydrostatic pressure gradient with the bed slope in a single term that depends on the free surface gradient. With this approach, when the free surface is horizontal, an exact balance between the bed slope and the hydrostatic pressure gradient is obtained naturally. As a result, more stable, efficient, and faster simulations (up to 1.5-times) are achieved as compared with other traditional FVM schemes [25,58,59].
Two additional options were also used in the simulations. On the one hand was an option that removes local depressions generated during the discretisation process along the flow path [61]. This tool is based on the technique proposed by Jenson and Domingue [62] in which each depression is refilled considering the lowest elevation of its neighbour cells. The application of this methodology is convenient even if the Digital Terrain Model (DTM) was previously treated with a Geographical Information System (GIS) software due two factors: (1) common GIS tools obtain the flow paths based on the connections between cells through their edges defined by the vertices, while the FVM uses the cells elevation to compute the flow; and (2) the topographical discretisation in the numerical model can use elements larger than the DTM raster cell size, and in the obtention, new depressions might be generated. Thus, this option ensures a proper definition of the flow path to work with the FVM whilst avoiding spurious depressions. On the other hand, to handle the transition from wet to dry conditions, and vice versa, Iber implements a wet-dry method by defining a water depth threshold (ε wd ) below which a cell is considered to be dry. Low values of ε wd , such as 0.0001 m, as used herein in all cases, guarantee mass conservation and improve rainfall-runoff transformation as well as overland flow propagation [63].
Other particular parameters, as well as the numerical discretisation, are detailed in the description of each case study and in the Supplementary Material.

Case Studies
The role of the Manning coefficient roughness in the numerical modelling of the overland flow, performed with coupled hydrological-hydraulic distributed numerical models, was addressed, not referring to specific case studies around the world but to the numerical strategy performed. To that end, well documented studies that analyse the hydrological response of a basin were used from different points of views.
In Case Study 1, theoretical aspects related to the definition of the time of concentration and the role of the bottom roughness on the hydrological response were analysed. In Case Study 2, a distributed numerical model was calibrated by varying the roughness coefficient to fit the results to the hydrograph resulting from an aggregated model forced with synthetic hyetographs based on historical data. For that purpose, the parameters involved in the aggregated model were determined following the local administration recommendations. Finally, in Case Study 3, four different well-documented rainfall events in a monitored and regulated basin were used to force a distributed model aiming to analyse the role of the roughness coefficient in the overland propagation. The roughness coefficient was calibrated to adjust the simulated water elevation to the observed one during each event. A total of six basins, with different topographic and hydrologic characteristics, were analysed ( Figure 1).

Case Study 1: Adjustment of the Roughness Coefficient Based on the Time of Concentration
In Case Study 1, the hydrological behaviour of four basins presented by Grimaldi et al. [10] is analysed. The four basins, located in USA, are: Brazos basin (Cow Bayu), San Antonio (Escondido Creek), Trinity basin (North Creek), and Brazos basin (North Elm Creek). Table 1 summarizes the main geometric characteristics. In Brazos (Cow Bayu) and Sant Antonio, there are large urban areas and both basins are regulated by a reservoir located at their lower parts. Trinity basin is a semi-rural area with five reservoirs that play a significant role in their hydrological behaviour. Finally, the largest basin, Brazos (North Elm Creek), is located in a rural area with two reservoirs (see Figure S1). In this case study, a distributed coupled model was built for each basin using Iber. Details of the parameters and structure of the numerical models are described in the Supplementary Material. No aggregated hydrological modelling was used.
The numerical models were operated with constant rainfall intensities of infinite duration, analysing the output hydrographs (S-shaped). The required time for the hydrograph to stabilize to constant discharge, i.e., the time for the basin to achieve a stationary state, should coincide with the time of concentration if the WMO definition [14] is adopted. For the adjustment of the roughness coefficients based on t c , the six common expressions used by Grimaldi et al. [10], plus the Témez [17] formula commonly used in Spain, were selected. At this point, the question of ambiguity on the concept and formulas to estimate t c arises again as pointed out by Beven [12] and commented upon in the introduction. If the variable used to compare methodologies is the time the S-curves need to stabilize, but the formulas used to estimate the t c are empirically derived using the SCS-UH method or other similar means, there is a conceptual discrepancy as explained by Beven [12].
The different formulas used to estimate t c together with the results are presented in Table 2. It also shows the minimum, maximum, and mean values of the time of equilibrium that are going to be compared with the results of the distributed method.

Case Study 2: Adjustment of the Roughness Coefficient Based on the Peak Time and Discharge from Aggregated Hydrological Models
Case Study 2 corresponds to the Mexican basin called Marquelia, located in the Guerrero region ( Figure 1). The basin precipitation is monitored by eight rainfall gauge stations within or near it, and a hydrometric station at the outlet (see Figure S3a). The area of the Marquelia basin is 1016 km 2 and its mean slope is 0.75%. Land uses are characterized by large extensions of secondary vegetation (48.51%) and forested (24.23%) areas. Other areas are covered by rivers (10.70%), pasturelands (10.60%), and rainforest agriculture extents (5.42%). According to the World Reference Base for Soil Resources [67], the basin edaphology is mainly composed of large areas of regosol (71%) and other types of soil such as phaeozem (11.3%) and cambisol (9.6%). From the land uses classification, a Curve Number (CN) of 76.7 was estimated as an area weighted average value in the whole basin.
Historical rainfall data was used to generate synthetic hyetographs associated with a 50-years return period (T50). Alternate block synthetic hyetographs [16] for a T50 event were calculated for each rain gauge station using the Intensity-Duration-Frequency (IDF) curves [68] and the time of concentration of 12.5 h resulting from the application of the Kirpich formula. Historical data of annual maximum precipitation were corrected, as proposed by Campos-Aranda [69], to account for the effect of using a fixed observation time-interval [70] and the non-uniformity of rainfall in the whole area of the basin.
On one hand, the output hydrograph was first calculated with aggregated modelling using HEC-HMS. The basin was discretized in seven subbasins. In each one, t c was estimated with Kirpich formula and t lag as 0.6t c . The subbasins were connected by three river reaches where the flow was propagated using the Muskingum routing method [16].
On the other hand, Iber was used for coupled hydrological-hydraulic distributed modelling. In this case, the terrain roughness was calibrated to adjust the hydrological behaviour (peak time and peak discharge) to that of the aggregated model. Six different Manning coefficient (n) values related to the land uses distribution of the basin (see Figure S3b) were used. All information regarding the topography, land uses, and precipitation data are presented in the Supplementary Material.

Case Study 3: Adjustment of the Roughness Coefficient Based on Observed Storm Events
Case Study 3 focuses on La Muga basin (Spain). Four well-documented rainfall events [58,59,71] were used to assess the hydrological response of the basin by varying the roughness coefficient values in the distributed model.
The study area consists of the upper part of La Muga basin, which has an extension of 181.2 km 2 and a reservoir at the outlet of 61 hm 3 of storage capacity. Weather conditions of the basin are characterized by a wide variability of the rainfall regime, which is influenced by marine conditions with small thermal as well as rain variations [72]. Extreme weather conditions-such as heavy rain events with high precipitation intensities and water accumulations concentrating in a few days or hours, and water scarcity associated to long dry-weather periods-are typical of the area [73,74].
A land use analysis [75] revealed that the basin is mainly covered by large forest extensions (92.5%), while a terrain characteristics analysis indicated a low permeability and underground storage capacity [72]. Thus, the hydrological response of the basin can be characterized with a unique land use and soil texture.
In the same way as for Case Study 2, the NRCS-CN method was applied to transform rainfall into runoff. Previous studies [58,59,71] allowed determining the CN value as a function of each studied rainfall event and the soil moisture antecedents evaluated with remote sensing [76,77]. The resulting CN values varied from 50 to 81. The four events were labelled with the starting date (year/month/day) and the duration in days as follows: 20130304_3d, 20131115_4d, 20141128_2d, and 20150320_6d. Table 3 summarizes the cumulated rainfall for each event, the maximum registered rainfall intensity, and the CN value. A more extended description of these events is detailed in the Supplementary Material. Based on the uniformity of land uses and soil characteristics in the basin, a unique Manning coefficient value was used. The model calibration was performed by comparing the calculated and observed water elevation at the Boadella reservoir dam, located at the basin's lower end (see Figure S5).

Case Study 1
Four distributed hydrological models were built, one for each of the four basins. The domain was discretized with a calculation mesh of a regular grid generated directly using the DTM data, being quadrilateral elements with sides of approximately 9.5 m. A constant rainfall intensity of 25 mm/h was assigned to the whole domain. In this case, no infiltration losses were considered; thus, the rainfall intensity was considered to be effective precipitation. The value of the Manning coefficient (n) was varied from 0.01 to 0.2 s/m 1/3 , in intervals of 0.01 s/m 1/3 . The purpose was to obtain the best approximation for the minimum, the maximum, and the mean times of concentration (t c ) calculated with the empirical formulas ( Table 2). The maximum discharge generated with the aforementioned rainfall intensity is, in this case, a constant discharge. For the estimation of t c , the discharge was considered to be constant when its variation between two consecutive results time steps was lower than 0.1% of the discharge. Results time steps of 60 s were used. Figure 2 shows the simulated hydrograph in each basin for the selected Manning coefficients, together with its value. All of them are consistent with the S-curve shape characteristic of the rational method [16,78], which stabilizes at the same value obtained with the rational method, i.e., the product of the effective rain intensity and the basin area (see Table S1). For the Brazos (Cow Bayu) ( Figure 2a) and San Antonio (Figure 2b) basins, a very low value of the Manning coefficient (0.01 s/m 1/3 ) was needed for the S-curve to fit with the minimum t c , which, in this case, corresponds to the Kirpich formula. In Trinity (Figure 2c) and Brazos (North Elm Creek) (Figure 2d) basins, it was not possible to reproduce the minimum t c without considering any friction since the hydrograph needed more time than the minimum t c to stabilize. For this reason, they are not plotted in Figure 2. This induces the thought that the Kirpich formula, the one that provides such minimum values, is quite unrealistic for these basins, as water propagation is slower even with no friction. This is in agreement with the assertion of Michailidi et al. [11] that "the Kirpich formula is a special case of a very general expression that is valid under very limited conditions" and that this formula was obtained for basins where channel flow was predominant.
The particular shape of each S-curve depends on the geometry and hydrological behaviour of each basin. The shape of the basin and the river network play an important role not just in the hydrograph definition, but also the drainage area. For the analysed basins, the smaller ones (Figure 2a,b) follow a more standard S-shaped hydrograph; in contrast, in the bigger ones (Figure 2c,d), the contribution of each part of the basin is most clearly appreciated and the hydrographs' shape differ more from a smooth S-shape. In particular, the hydrographs of Trinity basin (Figure 2c) show some more flat areas (around 2 h for n = 0.05 s/m 1/3 and between 3 and 5 h for n = 0.14 s/m 1/3 ) after an initial rapid discharge increase, which can be attributed to the effect of small basins having reached their maximum discharges.
Numerical results show a wide range of the Manning coefficient, from values lower than 0.01 up to 0.2 s/m 1/3 in order to achieve similar responses of the basin to those obtained by common time of concentration formulas and an aggregated model. As expected, a low value of the Manning coefficient provides faster response, while high values in a slower runoff advance. It is worth mentioning that a unique roughness coefficient value was used all over the basin; thus, a more detailed roughness discretisation would result in higher values at hillsides and areas outside of the hydrographic network, and lower ones in rivers and floodplains.

Case Study 2
The hydrological analysis of the Marquelia basin was performed first with an aggregated approached using HEC-HMS and the UH method, and then with a distributed model using Iber. The first model, built following the existing recommendations in Mexico, was taken as the reference for adjusting the roughness coefficients in the second one. In this, the domain was discretised with seven subbasins and three reaches. The time of concentration for each subbasin, t c,i , was determined by the Kirpich formula. The Muskingum method was used as the routing method [16], the parameters-K and χ-ere estimated following the recommendations of Fuentes et al. [79]. A more extended description is detailed in the Supplementary Material.
On the other hand, for the distributed approach, an irregular calculation mesh made of triangular elements, between 100 and 500 m of side, was built. Smaller elements were used in the river network while greater ones were used on the hillsides, aimed at optimizing the computational time and results accuracy balance. The topographical data used was the 15 m cell size DTM provided by the Mexican National Institute of Statistics and Geography, INEGI [80]. The Manning coefficients that associated the spatial distribution of the land uses map ( Figure S3b) were calibrated. Additional information of the model discretisation is detailed in Supplementary Material.
The values of the Manning roughness coefficient that provide the best adjustment in terms of peak discharge are depicted in Table 4. They range from 0.069 s/m 1/3 , corresponding to urban settlements, to 0.207 s/m 1/3 corresponding to forested area. All these values widely exceed the traditionally recommended values for hydraulic computations in artificial or natural channels and floodplains areas [38][39][40][81][82][83]. Table 4. Manning roughness coefficient according to the land uses discretisation shown in Figure S3b that provides the best adjustment in the calibration process.  Figure 3a plots the hydrographs resulting from the aggregated approach (dashed line) and the distributed approach (continuous line). The peak discharge and the hydrograph volume differ, in absolute value, by less than 0.2% and 3%, respectively, being in both cases lower in the ones obtained in the distributed approach. The aggregated model provides a faster hydrological response with almost two partial peaks discharges at 8.5 h and 11.5 h, while the maximum peak discharge is produced 1.6 h after the one obtained by the distributed model. In this sense, the distributed model also produces a partial peak at around 8.5 h but with a lower discharge. Despite the whole basin being split into seven subbasins and three river-reaches, the hydrological response of this semi-aggregated model does not represent the basin response as accurately as the distributed one. In this case, the shape of the hydrograph is clearly influenced by the rainfall-runoff method (UH). This is especially denoted in the generation of a partial peak during the first hours and a fast decay of recession curve after the maximum peak discharge, the discharge being zero after 27 h.

Land
On the other hand, the distributed numerical modelling based on the 2D-SWE solution, which allows considering the full topographical and land uses characteristics of the basin and subbasins, provides more accurate results for the hydrograph shape. A sensitivity analysis was carried out by considering variations of ± 20% in the Manning roughness coefficients (Table 4). Figure 3b also presents the resulting hydrographs of this analysis (dotted and dashed lines), increasing the peak discharge while decreasing the time of the peak when the roughness coefficient decreases, and vice versa. It is worth noting that the peak discharge varies from −9.5% (+20% of n values) to 11.9% (−20% of n values), while the time of the peak and the hydrograph volume vary by exactly ± 8.8% and ± 0.5%, respectively. Thus, variations on the roughness coefficient are directly related to the time when peak discharge is produced, but the value of the peak discharge is not proportionally increased/reduced due to this variation.

Case Study 3
The model of La Muga basin was based on a spatial discretisation of triangular mesh elements with variable size that were finer in rivers (20 m) and coarser in hillslopes (200 m), constructed from a high resolution DTM of 2 × 2 m cell size [84]. As normal river discharge is small in comparison with the bankfull, the DTM, from a LiDAR flight, is a good approximation of the riverbed geometry.
The model was operated with variations on the Manning roughness coefficient (n), from 0.02 to 0.16 s/m 1/3 with steps of 0.02 s/m 1/3 , resulting in different hydrological responses. The computed water elevation in the reservoir at the basin outlet was compared with the field observations along the events and at the end of each one. Figure 4 shows the simulated water elevation for the four rainfall events analysed.
In Figure 4a, the results of 20130304_3d event are shown, which concentrates on the rainfall during day 2 and the first-half of day 3. In it, the simulated water elevation at the end of the event is similar to the observed one for all scenarios, with relative errors being less than 1.5%, and the value of n being 0.1 s/m 1/3 -one that achieves the best adjustment (0.005% of error). For this event, significant differences on the arrival time of the flood can be observed, with a time gap of around 12 h between the hydrographs corresponding to the minimum and maximum n values. The evolution of the water elevation during the second-half of day 2 clearly shows faster hydrological responses for low values of n. Figure 4b shows the results of the simulations for the 20131115_3d event. The simulations fit well with the observations in terms of final water resources (volume of water in the reservoir), with relative errors being less than 0.2%. In this case, where rainfall intensities around 50 mm/h were recorded in the middle of day 2 and at the beginning of day 3, the evolution of the water elevation was not very sensible to variations of the Manning coefficient, especially for n values greater than 0.06 s/m 1/3 . The 20141128_2d event (Figure 4c) shows the worst adjustment of all events in the reproduction of the water elevation evolution, but, by contrast, it also achieves a good performance in terms of water resources at the end of the event (relative error less than 0.6%). In this case, the three peaks of rainfall intensity registered generated an almost constant increment of water resources stored in the reservoir; however, the numerical model did not reproduce this trend, probably because in the model the rainfall distribution was considered to be homogeneous all over the basin, unlike the real case [59]. Only for a Manning coefficient between 0.06 and 0.12 s/m 1/3 , the water elevation at the end of the event shows a good adjustment with relative error being less than 0.16%. Finally, in the 20150320_6d event shown in Figure 4d, which lasted 6 days, although more than 80% of the rainfall was concentrated during day 2, the numerical model reproduced well the water resources volume at the end of the event, with relative error being less than 0.09%. In this case, the evolution of the water elevation was similar to the observations, but for a slightly sharper increment of the water elevation between day 2 and 3. This fact is probably because the rainfall registered at the gauge station overestimated the real precipitation as pointed out by Sanz-Ramos [59].
The results show that the dependency of the time arrival of the water front to the Manning coefficient values is more evident for low n values (less than 0.08 s/m 1/3 ). In general, the hydrological response of the model is more sensible for low values of n than for high ones, especially from the point of view of water resources. As expected, low values of n produce larger overland flow discharges in shorter times.
No remarkable differences were observed in the model results for values of n greater than 0.08 s/m 1/3 . The precision of the results is aligned with the assumption of a unique Manning roughness coefficient value for the whole catchment, which in this case, was appropriate for hydrological purposes. On the other hand, for values of n in the lower range, oscillations on the water surface elevation appear.
Several indicators were used to assess the adjustment between the observations and the simulated results-root-mean square error (RMSE), mean absolute error (MAE), and the coefficient of determination (R 2 ). Table 5 summarizes the performance of the model for the four rainfall events in terms of these indicators. A good adjustment is seen for 20130304_3d and 20131115_3d events, with low values of MAE and RMSE, but also with high values of R 2 . The performance of 20141128_2d event is, in general, the poorest-the simulated water resources are clearly underestimated during almost all the event. In contrast with the other events, as the hydrological response is slower in comparison with the observations, the indicators show the best adjustments are for low values of n. Finally, the 20150320_6d event presents the worst R 2 values. A peer-to-peer results analysis (RMSE and MAE) reveals that the model has, in general, a good adjustment, especially for values of n greater than 0.06 s/m 1/3 . It is worth noting that the influence of the rainfall antecedents (arrival time of the water front) and the spatial distribution of the precipitation (non-expected increase of water elevation at the end of the episode) can condition the previous statistics. Nevertheless, from the previous considerations, it can be asserted that the performance of the model is, in general, good enough for the purpose of the simulations.

On the roughness coefficient values and their effect on the hydrological response
In free surface hydraulic numerical modelling, shear terms due to friction express the contribution of the momentum change due to the energy dissipation produced not only by flow-boundary interactions, but also by sub-grid turbulence if, as is generally done, no turbulence model is used. Since for hydrological purposes, the effect of turbulence is generally neglected, only the friction of the fluid with the terrain represents the energy dissipation. I.e., in 2D-SWE-based coupled hydrological-hydraulic distributed numerical modelling, the bottom roughness conditions not only the flow routing in the river network, but also the overland flow on the areas where the rainfall-runoff process takes place (hillsides).
Hence, the Manning roughness coefficient can be considered as a "property" of the terrain material where the overland flow is propagated. In 2D-SWE-based models, anisotropic properties of the shear stress terms (see Equation (1)) can be due not only to the velocity field but also to the Manning coefficient [85][86][87]. Thus, Manning coefficient could be considered as a vector to represent the anisotropic properties of the terrain (Equation (2)), where n x and n y are the two components of the Manning coefficient vector in X and Y direction, and |n| its modulus.
However, there is no strong evidence that confirms that the Manning's anisotropy should be considered in hydrological studies because the overland flow usually propagates in the same direction, especially over the hillsides where the rainfall-runoff process is generated. By contrast, the time-space variability of the Manning coefficient, including the dependency of the main hydraulic variables (depth and/or velocity), must be considered in general; the coupled hydrological-hydraulic distribute modelling is the unique numerical strategy that can deal with this.
On the other hand, despite the fact that in distributed numerical modelling the topographical discretisation of the domain is more accurate than for aggregated numerical modelling, the hydraulic behaviour of the water for very low depths (millimetric scale or lower) is not properly reproduced by the solution of the SWE. In this situation, the flow scale is in the same order of magnitude of the terrain particles scale, and the fluid tends to move more slowly than for larger river-like flow scales. Thus, in order to avoid too-fast hydrological responses, the bottom roughness must be increased in relation to the values used in common hydraulic situations. This explains why high values of the bottom roughness, the Manning coefficient value (n), are required in Case Study 2 and Case Study 3 to obtain a good performance of the numerical model from a hydrological point of view. In contrast, very low values of n are needed in Case Study 1 for an accurate adjustment of the time of concentration obtained with the Kirpich formula. A n value of 0.01 was needed to reach a similar t c in Brazos and San Antonio basins, whereas in the other basins, this value was even lower than the previous one. This is in line with various authors' discussions on the appropriateness of the Kirpich formula in situations other than those for which it was obtained [10][11][12].
The generally high values of the roughness coefficient obtained in the three case studies are far from the values extensively used for hydraulic purposes in flood routing [40], but are within the range of values generally found in the literature for hydrological purposes [26,30,31,42,[88][89][90][91], especially for hillslopes and flood plains, where the rainfallrunoff process predominates.
On the role of the domain discretisation in the basin hydrological response In general, coarse meshes provide faster simulations but poorer resolution and accuracy. In contrast, finer meshes provide more accurate results while the computational effort increases. Thus, practitioners must deal with achieving a proper balance to obtain good results with a reasonable computational effort.
Improvements in topographical and rainfall data acquisition and the increase of the computational capacity, in particular those related to Graphical Processing Unit (GPU) computing in particular [89,[92][93][94], are continuously leading to the advance of the numerical modelling, especially in hydrology. This new framework pushes modellers to generate increasingly detailed coupled hydrological-hydraulic distributed numerical models, even using, if necessary, the full information of the DTM (as shown for example in Case Study 1) as the computation mesh and distributed rainfall data. However, mesh resolution affects the model hydrological response (peak discharge) and rainfall intensity affects the value of the time of concentration of the basin.
To illustrate the last two facts, a double sensitivity analysis was carried out on Marquelia basin. The results show that the peak discharge increases when the number of elements of the numerical model also increases (Figure 5a). An asymptotic trend towards 7000 m 3 /s is observed for calculation meshes with more than 1.9 million elements. This peak discharge increment is probably due to the better representation of the topography, in general, as well as the better definition of the river-channel, in particular, which improve the flow propagation representation. This last domain discretisation (1.9 M elements) represents approximately 19 elements per ha, which is still far from the values of flood studies where more than 1000 elements per ha is quite common [95]. On the other hand, the time of response of a basin depends not only on the Manning coefficient but also on rainfall intensity [11]. Higher rainfall intensities imply shorter hydrological responses, i.e., the time of concentration becomes smaller. In this case, the t c was estimated as the time from the end of the excess rainfall to the point of inflection of the hydrograph where the recession curve begins [96]. Using the same model discretisation for Case Study 2 and varying the return period of the rainfall hyetograph from 2 to 10,000, it is possible to observe how the time of concentration tends to a constant value of around 15 h for the statistically less frequent events (Figure 5b). This fact is aligned with Grimaldi et al. [10], who suggest that t c becomes quasi-invariant with respect to the rainfall intensity for flood events with a high return period and who also signalled the time of concentration dependency with the cell-size of the DTM when using NRCS method. Thus, the use of coupled hydrological-hydraulic distributed models, where the time of concentration is no longer part of the modelling process, might reduce the uncertainties related to the rainfall-runoff processes.
On the numerical approach Currently, there are three different strategies to carry out an overland flow analysisconceptual, aggregated, and distributed. Conceptual and aggregated approaches are the simplest and are based on the main basin characteristics, such as the extent, main river length, main river slope, and soil characteristics to estimate the rainfall-runoff threshold, etc., and summarize all hydrological processes for the estimation of the peak discharge associated to a rainfall intensity (e.g., Rational Method).
Aligned with the evolution in the state of the art of applied hydrology, the aggregated models integrate the previous conceptual models into numerical models. They are tools based on the physics of the problem that, by integrating empirical formulas for overland processes, first allowed simulating rainfall events, groundwater interactions, flow routing, etc. These models are still widely extended because they are fast, robust, and conceptually simple. However, as for other numerical tools, to achieve good results they must be calibrated. In them, one of the main parameters to calibrate is the time of concentration (t c ), which can be estimated by empirical formulas of one or the other authors. Additionally, practitioner's expertise is crucial to determine the number of subbasins to define the rainfall-runoff and routing methods, etc.; thus, in all of these processes, the error derived from selecting one or the other parameter is propagated till the final result: the hydrograph.
Finally, the more recently-developed distributed models, and in particular the coupled hydrological-hydraulic 2D-SWE-based tools, provide more detailed hydrologic as well as hydraulic results. In this case, the practitioner's expertise may be somewhat less relevant because there are less parameters to calibrate; in fact, only the terrain roughness. The domain discretisation is also important to achieve suitable results. It plays a similar role to the subbasins' definition in aggregated approaches but, in this case, the modeller has only to decide the mesh element size (minimum, maximum, or mean) and not how many subbasins/reaches the domain must be divided into.
As shown along the document, the value of the terrain roughness in 2D-SWE-based distributed approach is usually out of the range for hydraulic purposes as overland flow tends to move more slowly than when the flow is propagated along a river. Additionally, a distributed approach allows considering the terrain properties and characteristics element by element.

Conclusions
The bottom roughness, evaluated here via the Manning roughness coefficient, is one of the main calibration parameters for hydrological modelling based on coupled hydrologicalhydraulic distributed numerical tools. The mesh resolution also plays an important role on the overland flow propagation, especially on the definition of the peak discharge associated to a rainfall event.
Due to the ambiguities in the definition of the time of concentration and the high dispersion of the values obtained with different formulae, the adjustment of the Manning coefficients with the time of concentration is uncertain and of little practical interest. The Manning coefficients that have to be used in coupled hydrological-hydraulic modelling in order to reach a time of concentrations similar to those obtained by empirical formulas show an extremely high dispersion. Moreover, they are out of the common range of values used for hydraulic purposes. In particular, it was shown that for the Kirpich formula, the required value of the roughness coefficient for the time of concentration adjustment is close or lower than 0.01 s/m 1/3 , which would correspond to the hydraulic value for very smooth material, much different from those existing in natural basins.
The global picture of the presented results shows that higher values of the Manning coefficient than those commonly used for hydrodynamic purposes are required when computing the overland flow of a basin using coupled hydrological-hydraulic distributed numerical tools. This is because the roughness coefficient must integrate the momentum change due to the energy dissipation produced not only by flow-boundary interactions, but also by sub-grid turbulence, in the river network and in the areas where the rainfall-runoff is mainly generated (hillsides).
In coupled hydrological-hydraulic distributed models, the roughness coefficient has a clear role on the arrival time of the water front, playing a similar role like the lag time in conceptual and aggregated models. However, there is no exact relation of roughness values with commonly used hydrological parameters, which means that trying to estimate the roughness coefficient by comparing distributed and aggregated model results is not a good idea. From the authors' points of view, there is still little experience, though increasing rapidly, in the verification of coupled hydrological-hydraulic distributed models; thus, the recommendation would be to use this approach in cases with enough data for model calibration and validation while experience in this kind of model continues to increase and eventually more recommendations and guidelines appear. Example of the levelling process of the terrain in the area that represents a dam in Trinity basin. Calculation mesh of Iber before (a) and after (b) the levelling process, Figure S3: Characterization of the Case Study 2: (a) Topographical map of the basin, influence area of the most closed rain gauges stations (green stars) according to a Thiessen polygon discretization (black lines) and location of the hydrometric station (yellow pentagon); (b) Map of land uses and location of the outlet (red dot), georeferenced in UTM Zone 14N; (c) Map of soil edaphology; and (d) Map of Curve Number according to SCS-CN method, Figure S4: Case study 2. Synthetic hyetographs corresponding to a 50-years return period for the eight most closed rainfall station to Marquelia basin, Figure S5: Case Study 3. Topographical map of the upper part of La Muga basin and location of the gauge station (yellow pentagon) and the outlet (red dot), which is referenced in UTM Zone 31N, Figure S6: Case Study 3. 5-minutal time-resolution hyetographs rainfall intensity (grey bar) and outlet discharge from the reservoir (black line), Table S1: Maximum discharge (in m 3 /s) for the four basins of the Case Study 1 estimated by the Rational Method for a constant rainfall intensity of 25 mm/h, Table S2: Characteristics of the discretization of the Marquelia basin and parameters of the aggregated numerical approach, Table S3: Spatial characteristics of the land uses of the Marquelia basin and Manning roughness coefficient best fit to the 50-years return period synthetic hydrograph.