Investigating Parameter Transferability across Models and Events for a Semiarid Mediterranean Catchment

: Physically based distributed hydrologic models (DHMs) simulate watershed processes by applying physical equations with a variety of simplifying assumptions and discretization approaches. These equations depend on parameters that, in most cases, can be measured and, theoretically, transferred across di ﬀ erent types of DHMs. The aim of this study is to test the potential of parameter transferability in a real catchment for two contrasting periods among three DHMs of varying complexity. The case study chosen is a small Mediterranean catchment where the TIN-based Real-time Integrated Basin Simulator (tRIBS) model was previously calibrated and tested. The same datasets and parameters are used here to apply two other DHMs—the TOPographic Kinematic Approximation and Integration model (TOPKAPI) and CATchment HYdrology (CATHY) models. Model performance was measured against observed discharge at the basin outlet for a one-year period (1930) corresponding to average wetness conditions for the region, and for a much drier two-year period (1931–1932). The three DHMs performed comparably for the 1930 period but showed more signiﬁcant di ﬀ erences (the CATHY model in particular for the dry period. In order to improve the performance of CATHY for this latter period, an hypothesis of soil crusting was introduced, assigning a lower saturated hydraulic conductivity to the top soil layer. It is concluded that, while the physical basis for the three models allowed transfer of parameters in a broad sense, transferability can break down when simulation conditions are greatly altered.


Introduction
Physically based distributed hydrologic models (DHMs) are useful tools for simulating watershed hydrologic response in applications that include water resources estimation and management, flood forecasting, and support of numerical weather prediction models [1][2][3][4]. To account for the spatial variability of physiographic watershed characteristics, DHMs discretize the spatial domain into smaller elements and, in each of them, solve physics equations of the different processes. These equations depend on parameters that could potentially be measured in the field and whose values fall within often known admissible ranges. As a result, when different DHMs are applied to the same basin or the same DHM is applied to different events or even diverse sites, similar parameter values should, in principle, be adopted for the simulation of the same hydrologic processes. In practice, however, parameter transferability across DHMs, or for the same DHM at different locations or times, is still a challenging task for several reasons. First, the availability of land surface and hydrometeorological data used to force and validate a model is often limited, thus constraining the selection of model parameters when trying to match model outputs with observations [5]. Second, DHMs differ greatly in their representation of spatial domain (e.g., grids versus irregular meshes), in their level of sophistication in parameterizing physical processes, and in the numerical schemes used to solve their governing equations [6]. Third, the same DHM or different DHMs can be applied at different spatial and temporal resolutions [7]. As a result, the same parameter may assume different values if it is adopted in models with different structures and resolutions [8]. Finally, there is a gap between theory and hydrologic models stemming from the pragmatic manner in which models are often developed, limiting the ability to generalize about hydrologic behaviors and resulting in poor parameter transferability [9].
Several studies have examined the issue of parameter transferability in space (i.e., across basins) and time (i.e., across different hydrologic conditions or events) for the same DHM. Heuvelmans et al. [10] applied the semi-distributed SWAT agrohydrological model to three catchments in Belgium and found that the transfer of parameters across catchments with different environmental conditions could be problematic, especially for the simulation of low flows. Hogue et al. [11] applied the Noah land surface model to two semiarid sites in Arizona and showed that, in the absence of calibration data, a proxy basin set of calibrated parameters can be applied with a moderate decline in performance. Shrestha et al. [12] used the Blockwise use of TOPMODEL with Muskingum-Cunge routing (BTOPMC) rainfall-runoff model in Nepalese mountain catchments and found that calibrated parameters can be spatially and temporally transferred within catchments located in similar physiographic regions. Ren et al. [13] evaluated the potential of transferring optimized parameter values of the Community Land Model (CLM) land surface model across 431 basins in the United States and found that these are more transferable within sites that belong to the same class as regards climate, soil, and land cover/land use than between different classes, so long as the soil and climate conditions did not vary substantially within the class. Melsen et al. [7] applied the Variable Infiltration Capacity (VIC) land surface model at different spatial and temporal resolutions for a basin in Switzerland and found that parameters transferred more readily across spatial than temporal resolutions, suggesting that this could be caused by a substantial under-representation of the spatial variability of basin properties and hydrometeorological forcings in their simulations.
Model intercomparison initiatives (e.g., [14][15][16]) should be an ideal platform for addressing parameter transferability across models, but often in these studies each hydrologic model is independently calibrated. There are nonetheless several intercomparison studies that have focused on parameter transferability across DHMs. In Maxwell et al. [17] and Kollet et al. [18], different sets of seven coupled, physically based surface-subsurface hydrologic models were intercompared for a series of benchmark test cases, adopting common parameter values for the simulation of the same physical processes. Their results showed generally good agreement for the simple test cases, while the more complex cases brought out differences that were attributed to the varying level of detail in the mathematical and numerical representation of physical processes in the different models. Koch et al. [19] compared simulations of three DHMs against soil moisture patterns for a catchment in Germany using a surrogate model (HYDRUS 1D) to establish the common values for the soil hydraulic parameters in order not to bias the intercomparison. The parameter transfer from HYDRUS 1D to the more complex multi-dimensional models proved to be of limited success because the HYDRUS 1D results diverged from those of the three DHMs, and moreover, the transferred parameter values appeared to ignore the drying out of the soil. Rosero et al. [20] conducted a sensitivity analysis to explore the complex parameter interactions in three versions of Noah-the standard version, a version enhanced with a simple groundwater module, and a version augmented by a dynamic phenology module. They showed that transferability between sites and models is not straightforward, since the three model versions have different optimal distributions for the same physical parameters. Yen et al. [21] applied three different versions of SWAT to a watershed in Maryland and found that parameters were transferable temporally or spatially but not structurally (i.e., across the model versions).
The aim of this study is to assess, for a single study basin, the potential of parameter transferability across both models and events. Three DHMs that are conceptually and structurally very different in representing physical processes and terrain features are applied for two contrasting periods-a one-year simulation period representative of average hydrometeorological conditions for the study region and a two-year period under very dry conditions. The study site is the Rio Mannu basin (RMB), a small semiarid Mediterranean catchment located in southern Sardinia (Italy). The models are the TIN-based Real-time Integrated Basin Simulator (tRIBS; [22]), which was implemented for the RMB as part of a prior study [23], the TOPographic Kinematic Approximation and Integration model (TOPKAPI; [24]), and the CATchment HYdrology model (CATHY; [25]). Some of the main features of each model that highlight their differences are as follows: (i) tRIBS simulates a coupled water and energy balance based on irregular meshes, uses a simplified 1D scheme for water infiltration, and estimates actual evapotranspiration (ET) as a fraction of the potential ET on the basis of a threshold level of moisture deficit; (ii) CATHY is a grid-based model that simulates coupled surface and subsurface flow, uses the full 3D Richards equation for subsurface flow, and a boundary condition switching procedure to partition potential ET into actual ET when surface nodes reach a threshold level of moisture deficit; and (iii) TOPKAPI combines the kinematic approach with the topography of the basin, uses a 1D scheme for infiltration, and estimates actual ET as a function of the potential ET on the basis of the actual soil moisture content. Given that tRIBS was previously calibrated in the RMB, the same soil and land cover parameter settings are used for the other two DHMs in order to explore the validity of the hypothesis of direct parameter transferability for two contrasting simulation periods, under the assumption of parameter scale independence. The results are discussed in the context of the process representations for each model, and model performance is evaluated using the Nash-Sutcliffe index [26] and the mean absolute peak time error (MAPTE) [27].

Study Area and Dataset
The study site is the Rio Mannu di San Sperate basin located in southern Sardinia, Italy ( Figure 1). This catchment has been extensively studied through numerical simulations and field campaigns during the European research project Climate Induced Changes on the Hydrology of Mediterranean Basins (CLIMB [28]). In particular, a series of field campaigns was conducted during the project in 2011, during which a total of 50 soil samples of 80 cm depth were collected throughout the watershed and analyzed to characterize soil texture and to reduce the uncertainty in the soil texture classification [23]. Another factor that contributed to the selection of this study site is that this basin has been affected by prolonged drought periods during the last 30 years that have led to water restrictions for the agricultural sector, with significant financial losses for the region of Sardinia.
The RMB drains an area of 473 km 2 , with elevation varying from 66 m to 963 m a.s.l. and a mean slope of about 17%. The main channel is about 39 km long, with a concentration time of 12 h. The basin is characterized by a typical Mediterranean climate with extremely dry summers and rainfall from October to April. Figure 2 shows the mean monthly values of discharge, precipitation, and temperature in the catchment, computed from daily data during the period from 1925 to 1935, while additional subdaily meteorological data in more recent periods were used for calibration of the weather generator (see later in this section). The discharge regime is characterized by low flows (less than 1 m 3 /s) for most of the year and sporadic flood events in the fall and winter, due mainly to frontal precipitation systems [29,30]. According to the available data, the mean annual precipitation is 680 mm and the mean monthly temperature ranges from 9 • C in January to 25 • C in July-August.      Geospatial datasets characterizing the catchment land surface properties include a digital elevation model (DEM) at 10 m resolution and a CORINE land cover (LC) map [32] for the year 2008 reclassified into eight groups. As shown in Table 1 and Figure 3b, the dominant land use classes are agriculture (48%), concentrated in the valley area that is part of the Campidano plain, and sparse vegetation (26%). The most important crops grown are wheat, maize, beans, and artichokes [33]. The pedologic map of Aru et al. [34] was used to characterize the soil properties, and its classes were aggregated on the basis of field campaigns conducted as part of the aforementioned European project (Table 1 and Figure 3a; [35]). A soil depth map was obtained by combining the DEM with soil texture information, following Wigmosta et al. [36]. Table 1. Land cover and range of soil texture classes used as input for the three hydrologic models, with the corresponding percentage of basin area (adapted from Mascaro et al. [23]). Geospatial datasets characterizing the catchment land surface properties include a digital elevation model (DEM) at 10 m resolution and a CORINE land cover (LC) map [32] for the year 2008 reclassified into eight groups. As shown in Table 1 and Figure 3b, the dominant land use classes are agriculture (48%), concentrated in the valley area that is part of the Campidano plain, and sparse vegetation (26%). The most important crops grown are wheat, maize, beans, and artichokes [33]. The pedologic map of Aru et al. [34] was used to characterize the soil properties, and its classes were aggregated on the basis of field campaigns conducted as part of the aforementioned European project (Table 1 and Figure 3a; [35]). A soil depth map was obtained by combining the DEM with soil texture information, following Wigmosta et al. [36]. Table 1. Land cover and range of soil texture classes used as input for the three hydrologic models, with the corresponding percentage of basin area (adapted from Mascaro et al., [23]).

Range of Soil Texture Classes % Basin Area Land Cover Class % Basin Area
Sandy  The hydrometeorological records for this catchment are limited and sparse, and were collected at different resolutions and over discontinuous time periods. Daily discharge data at the outlet section were collected from 1925 to 1935. During this period, 12 rain gauges provided daily precipitation records, while one thermometric station located in Cagliari (the circle in Figure 1b) registered daily minimum (Tmin) and maximum (Tmax) temperatures. To provide forcings at resolutions needed for the application of tRIBS, Mascaro et al. [23] developed two downscaling strategies to generate an ensemble of 50 equally probable hourly precipitation grids at a spatial The hydrometeorological records for this catchment are limited and sparse, and were collected at different resolutions and over discontinuous time periods. Daily discharge data at the outlet section were collected from 1925 to 1935. During this period, 12 rain gauges provided daily precipitation records, while one thermometric station located in Cagliari (the circle in Figure 1b) registered daily minimum (T min ) and maximum (T max ) temperatures. To provide forcings at resolutions needed for the application of tRIBS, Mascaro et al. [23] developed two downscaling strategies to generate an ensemble of 50 equally probable hourly precipitation grids at a spatial resolution of~13 km and basin-averaged hourly potential evapotranspiration (ET 0 ). Coarse precipitation measured by the 12 rain gauges was statistically downscaled in space and time using the Space Time RAINfall (STRAIN) multifractal model [37,38], adapted for orographically induced heterogeneity as described in Badas et al. [39]. The STRAIN model was applied to disaggregate precipitation from the coarse resolution of 104 km, 24 h to the fine resolution of 13 km, 1 h, which are the same ranges of temporal and spatial scales for which Badas et al. [39] proved the existence of scale invariance laws. Each member of the 50 ensemble simulations represents a statistically equiprobable scenario characterizing the sampling uncertainty from the same coarse-scale conditions. In particular, for each member of this ensemble the total daily precipitation volume was conserved, while the spatiotemporal distribution randomly changes [23]. The downscaling model was calibrated with high-resolution (1-min) rain gauges available over the period 1986-1996 in the coarse domain of 104 km x 104 km that contains the RMB. Basin-averaged hourly estimates of ET 0 were instead obtained through a downscaling method based on dimensionless functions that reproduce, for each month, the average daily cycle of ET 0 using T min and T max as inputs.
This routine was calibrated and tested by applying the Penman-Monteith equation [40,41] with hourly meteorological data collected in Cagliari from 1995 to 2010. More details on the downscaling tools can be found in Mascaro et al. [23].

Hydrologic Models
The tRIBS, TOPKAPI, and CATHY models are physically based DHMs that differ greatly in their representation of physical processes, numerical complexity, and terrain features. Table 2 summarizes the characteristics of each hydrologic model and highlights some of their main differences. Topographic representation differs between the three models; terrain features are represented in tRIBS via triangulated irregular networks (TINs), which define the model domain through associated Voronoi polygons, while TOPKAPI and CATHY use a regular grid. The infiltration scheme in tRIBS is based on a two-dimensional modified Green-Ampt model [42], while a kinematic wave routing model is used to simulate transport of water in the channel network. In each basic computational element of the tRIBS model (i.e., the Voronoi polygon), the governing equations are solved using a finite difference control-volume approach [22]. The model considers a soil column of heterogeneous, sloped, and anisotropic soil in which the saturated hydraulic conductivity decreases with normal depth [42]. For the TOPKAPI model, four nonlinear reservoir differential equations, derived from the integration in space of the kinematic wave model [43], are used to describe subsurface, overland, and channel flow. In TOPKAPI, subsurface flow can occur in two layers-superficial and deep. The superficial layer is characterized by thin thickness and high hydraulic conductivity because of the macroporosity, and it plays a key role in direct flow contributions to the drainage network and in the activation of the saturated area, which causes surface flow. The deep layer is characterized by a greater thickness and reduced hydraulic conductivity due to the more compact soil and plays a key role in determining infiltration and base flow. CATHY simulates instead in a detailed manner the interactions between subsurface and surface water [25]. The surface module is based on the resolution, using the finite difference method, of a one-dimensional diffusion wave approximation of the Saint Venant equation for overland and channel routing [44]. The subsurface module solves, using the finite element method, the three-dimensional Richards equation that describes flow in variably saturated porous media [45]. Each cell of the surface grid is subdivided into two triangles, which are projected along the vertical into a number of layers whose value and thickness can be chosen by the user to create the 3D tetrahedral finite element mesh [25]. In particular, the model needs a finer layer resolution close to the surface to capture the water flux partitioning. In tRIBS, potential ET can be computed by solving the energy balance through the Penman-Monteith approach [40,41], or alternatively, the model can be forced with external input, while the actual ET is estimated as a fraction of the potential ET, based on the soil moisture available in the soil. Potential ET in TOPKAPI is estimated using the Thornthwaite and Mather [46] formula, or it can be provided as an external input, while actual ET is computed with the radiation method [47]. CATHY requires potential ET values as external input, while a boundary condition switching procedure is used to partition potential (atmospheric) fluxes into actual fluxes across the land surface (infiltration, exfiltration as evaporation, and exfiltration as return flow) and changes in surface storage (ponding heads) [25].
The three hydrologic models account for spatial variability of soil properties and can be forced by spatially variable precipitation and meteorological data at hourly resolution. The tRIBS model has been applied across a large range of scales (from small to regional basins) in the areas of hydrometeorology, climate change, and ecohydrology (e.g., [48][49][50]). Recently, Mascaro et al. [23] calibrated the tRIBS model in the RMB, and Piras et al. [35,51] applied it to evaluate the hydrologic impacts of climate change. TOPKAPI has been successfully implemented as a research and operational hydrologic model in several catchments worldwide (e.g., [52][53][54][55]). The CATHY model has been used in many exploratory studies, benchmarking exercises, and real catchment applications (e.g., [18,31,[56][57][58][59]).

Setup of the DHMs
We first summarize the procedures adopted by Mascaro et al. [23] to set up, calibrate, and validate the tRIBS model in the RMB. We then describe how the parameters of TOPKAPI and CATHY that are not shared with tRIBS were assigned for this present study.

tRIBS
The spatial domain for the tRIBS model was created by generating a TIN with an equivalent resolution of~50 m from the 10 m DEM, following the approach described in Vivoni et al. [60]. The model was calibrated against discharge observations at the outlet for a one-year period (1930) and validated for a two-year period (1931)(1932). These three consecutive years were chosen by Mascaro et al. [23] among the available years 1925-1935 because the published rating curves did not vary significantly over this three-year period. This strategy is based on Klemes [61], who suggested the use of a wet period for calibration and a drier period for validation, especially when the hydrologic model is used to evaluate the effects of climate change [62]. A spin-up interval of two years was used for tRIBS, following Vivoni et al. [63], to set the initial conditions. The total observed rainfall during the calibration year and the two validation years, obtained from the 12 rain gauges within the catchment boundaries via the Thiessen polygon method, was, respectively, 902 mm during 1930 and 577 mm and 450 mm during 1931 and 1932, while the total potential evaporation was 898 mm during 1930 and 933 mm and 895 mm during 1931 and 1932. The calibration was limited to the two most sensitive model parameters: saturated hydraulic conductivity at the surface in the uppermost part of the soil (K s ) and conductivity decay parameter (f), which dictates the variation of K s with soil depth [42]. For each soil class (see Figure 3a), the calibrated values of K s and f were selected within physically plausible ranges taken from Rawls et al. [64]. All other soil and vegetation parameter values were assigned from the literature [64][65][66][67] and are reported in Mascaro et al. [23]. The resulting parameter values for the main soil classes are presented in Table 3, while the routing parameters are reported in Table S1. These are considered reference values for the TOPKAPI and CATHY models. The use of a CORINE land cover map from 2008 as a reference for the 1930-1932 simulation period was also assessed by Mascaro et al. [23], who found minimal long term changes in vegetation cover and urbanization over the study region.

TOPKAPI and CATHY
The spatial domains of TOPKAPI and CATHY consist of regular grids. These were obtained by aggregating the 10 m DEM at resolutions of 250 m. The coarser resolution for the regular-grid TOPKAPI and CATHY models compared to the 50 m (on average) resolution for the irregular-grid tRIBS model was a necessary compromise between accuracy and computational requirements. Under the hypothesis of direct parameter transferability that underpins this study, an assumption of parameter scale independence was also made on account of the grid discretization differences between the models. The parameter values for both models were assigned based on the same calibration and validation periods used for tRIBS, using the same soil texture and land cover maps. Both models were forced with the same hourly precipitation and potential evapotranspiration datasets used for tRIBS. In TOPKAPI, the soil was discretized into a superficial layer with a depth of 0.5 m and a deep layer that ranges from 0.5 m to 1.5 m. To parameterize TOPKAPI, for the superficial layer the same values for the saturated hydraulic conductivity at the surface K s , saturated moisture content θ s , residual moisture content θ r , and air entry pressure head ψ a used for tRIBS (see Table 3) were assumed, while for the deep layer the K s values were set one order of magnitude lower, following suggestions from the model developers. The routing parameters for the TOPKAPI model are reported in Table S2, where the maximum channel width and the Manning coefficient for the channel Strahler order 1 were set on the basis of values assumed for tRIBS. A spin-up interval of two years was used, as for the tRIBS model, to set the initial conditions. For the CATHY model, the soil was discretized into 12 layers, with thinner layers at the surface where the partitioning of precipitation into surface water and infiltration takes place, for a total depth for each cell assigned according to the soil depth map. The thinner surface layers allow for better accounting of the interactions between surface and subsurface water. The bottom and lateral boundaries of the domain were assumed to be impermeable, as were those of tRIBS and TOPKAPI. The same values for K s , θ s , θ r , ψ a , f, and the van Genuchten and Nielsen [68] pore size distribution index m used for tRIBS were assumed. As with tRIBS and TOPKAPI, therefore, in the CATHY model the saturated conductivity also decreased with depth. The threshold level of moisture deficit (in terms of pressure head), P min , for the CATHY model was set to −3 m, which corresponds to the threshold in soil moisture, θ*, used for the tRIBS model to identify the stress threshold for evaporation and transpiration. The channel network was identified from the DEM of the catchment using an upstream drainage area threshold of 5 km 2 on the basis of visual similarity between the extracted network and the streamlines depicted on topographic maps. The structural parameters for the channel and overland flow networks were identified based on values assumed for tRIBS (Manning coefficient and water surface width for channel flow) and reported in literature studies as a basis [69][70][71]. The assigned values obtained are reported in Table S3. For CATHY, the initial conditions (pressure heads at each node of the domain) were derived by draining the catchment from full saturation and with no atmospheric forcing until the flow at the outlet matched the observed discharge at the beginning of 1930.

Metrics for Model Evaluation
Hydrologic model performance was compared during the two considered periods, the one-year simulation period (1930) representative of average hydrometeorological conditions for the basin and the two-year period (1931)(1932) under very dry conditions. The three hydrologic model results were qualitatively discussed in terms of magnitude of peak, low flows, and receding limb shape, and quantified in terms of mean absolute peak time error (MAPTE) [27] and total runoff volumes. In addition, the Nash-Sutcliffe index (NS) [26] was evaluated using the observed and simulated discharge volumes at different aggregation times (daily, weekly, and monthly; performance is best for highest NS values, i.e., closest to 1).

Results and Discussion
Different sets of simulations using as atmospheric input the ensemble of 50 downscaled precipitation series that were assumed as forcing for the tRIBS model from the previous study of Mascaro et al. [23] that were carried out for the TOPKAPI and CATHY models for the 1930 and 1931-1932 periods. Figure 4 shows the comparison of simulated and observed discharge in 1930 for the three models presented as 90% confidence intervals across the 50 members. The two insets offer a more detailed comparison for two subperiods with significant flood events. The difference between the downscaled ensemble average (MAP D ) and observed (MAP O ) mean areal precipitation at the daily scale is also plotted for each inset. It can be seen that, despite the uncertainty in hydrometeorological inputs, the three hydrologic models reproduce the overall shape and timing of the major flood events. However, there are some differences, described also in Mascaro et al. [23] for the tRIBS application. The differences are due in part to the use of a stochastic downscaling tool that, in redistributing daily rainfall volumes from a large domain to smaller areas and times, does not always preserve the exact spatial localization of the storms, and in part to the different temporal resolution between simulated (sub-hourly or hourly) and observed (one measurement per day) discharges. As a consequence, we observe that none of the models is able to reproduce the peaks labeled as M (missed) due to a previous period of underestimated precipitation (negative MAP D -MAP O ). CATHY and TOPKAPI simulate quite well the recession curves, while tRIBS has a tendency to simulate steeper recession limbs [62], probably due to higher terrain resolution adopted in this model, which better captures steeper terrain slopes, thus leading to higher water velocities in the routing schemes of hillslopes and channels. Moreover, the timing of flood peaks can also be affected by the differences in observed and simulated MAP, as illustrated by the event labeled D (delayed). The tRIBS model misses the peak of hydrograph seven times, and CATHY and TOPKAPI five and three times, respectively. We can see in Figure 4 that CATHY is unable to reproduce the final discharge peak in 1930 because, with the drying out of the soil after the warm months, this model with its detailed subsurface representation allows greater infiltration and thus less runoff generation for this subperiod. Table 4 reports the MAPTE values for the three hydrologic models. We see that the CATHY model produces the minimum average of all absolute peak time errors during the 1930 period, while the TOPKAPI model yields the highest value. considered, NS increased and reached a mean value of 0.55 at monthly resolution for the tRIBS model, 0.61 for TOPKAPI, and 0.62 for CATHY. To better illustrate the model performance, Figure 5 shows the comparison between the observed flow duration curve (FDC) and the 90% confidence intervals from the 50 ensemble simulations for the three hydrologic models. The shape of the observed FDC is well reproduced for the major flood events for all models. However, tRIBS underestimates the discharge values corresponding to the 5-10% exceedance level due to the previously mentioned tendency to generate steeper recession limbs.  In terms of total runoff volume, the ensemble mean for the three hydrologic models, together with the standard deviation across the 50 members, is reported in Table 5. The minimum, mean, and maximum discharge volume statistics for the 50 ensemble members during 1930 are reported in Table 6. The lowest values of NS were obtained at daily resolution because, at this scale, the direct correspondence between observations and simulations was more affected by the different sampling time step and by the mismatch in the disaggregated forcing [23,62]. When larger time scales were considered, NS increased and reached a mean value of 0.55 at monthly resolution for the tRIBS model, 0.61 for TOPKAPI, and 0.62 for CATHY. To better illustrate the model performance, Figure 5 shows the comparison between the observed flow duration curve (FDC) and the 90% confidence intervals from the 50 ensemble simulations for the three hydrologic models. The shape of the observed FDC is well reproduced for the major flood events for all models. However, tRIBS underestimates the discharge values corresponding to the 5-10% exceedance level due to the previously mentioned tendency to generate steeper recession limbs.

Wet Period Simulations
Water 2019, 11, 2261 11 of 21 derived from 50 ensemble simulations for each of the models. In the insets, a zoom on two periods with significant flood events is shown in order to better highlight the differences between the daily downscaled ensemble average (MAPD) and observed (MAPO) mean areal precipitation time series. The peaks missed or delayed are labeled as M and D, respectively.

Dry Period Simulations
The results for the drier 1931-1932 period are shown in Figure 6. The tRIBS and TOPKAPI models perform well for the discharge time series over 1931 and most of 1932. In the subperiod October-December 1932, the two models simulate a number of peaks that were not observed due to precipitation overestimation in the previous periods, and they sometimes underestimate the discharge, for similar reasons. CATHY is unable to generate any discharge except at the beginning of 1931, likely for the same reason regarding subsurface representation mentioned earlier. To circumvent this and attempt to generate discharge in the CATHY simulations, we conducted extensive tests by adjusting the soil parameter P min . As shown by Camporese et al. [72], CATHY is quite sensitive to this value that sets the boundary condition switching between potential (atmosphere-controlled) and actual (soil-limited) evapotranspiration. In particular, P min was increased from the value transferred from tRIBS (P min = −3 m) to −1 m and −0.5 m, but none of these efforts was successful in reproducing the observed discharge. Figure 7 reports the results using a member of the 50 ensemble simulations during the wet 1930 and the drier 1931-1932, with the different values of P min . It was therefore concluded that the subsurface representations between tRIBS (on which the calibrated values were based) and CATHY are too different for a direct transfer of parameter values between the two models to be valid under all circumstances. In particular, whereas tRIBS accounts for subsurface flow via a single layer wherein an approximation to the 2D Richards equation is solved, CATHY solves the full Richards equation, in 3D and over several layers (12 for this study). It is likely that this finer, multi-dimensional representation allows for greater infiltration during rain events, owing to the increased numerical resolution and the greater allowance for lateral moisture redistribution. We believe this difference in physical representation between the models to be a more important factor in this case than numerical discretization differences (e.g., coarser resolution DEM compared to the TIN resolution). This dry period case, very different from the simulation period used for model calibration, therefore reveals some of the key structural differences between the models and, consequently, limitations with regards to parameter transferability between models. In further investigating the difficulty to generate discharge with CATHY under the very dry conditions that prevailed during the 1931-1932 period, we hypothesized a process of soil crusting that, numerically, can be simulated by lowering the hydraulic conductivity of the topmost soil layer during the dry years. This represents a variation (i.e., non transferability) of a parameter value for CATHY not just with respect to the other models but also for different events. While the hypothesis of soil crusting provides a physical basis for a temporal change in parameter value, a feature readily incorporated in the CATHY model, further study of the differences between the three models would be needed to explore in greater depth this additional aspect related to parameter transferability.  The physical basis and modeling of soil crusting phenomena have been investigated by Assouline [73], who provides a review of various approaches for modeling infiltration into sealed soils. In earlier work, Hillel and Gardner [74,75] hypothesized that a sealed soil could be modeled as a soil profile capped with a saturated thin layer of low permeability. This simple approach has been applied in different studies [76][77][78][79] and is suited to our purposes here. Figure 8 shows the comparison of simulated and observed discharge for the 1931-1932 period for the CATHY model with the assumption of soil crusting, while Table 5 summarizes the total simulated runoff volume for the three hydrologic models. From Table 4, we can notice that the tRIBS model produces the minimum average of all absolute peak time errors during the 1931-1932 years, while the TOPKAPI model yields the highest value, as it did for the 1930 period.      Table 7 reports the Nash-Sutcliffe index for the three hydrologic models. With the assumption of soil crusting, the NS values obtained with the CATHY model are comparable to those of the other models. Figure 9 shows the comparison between the observed flow duration curve and the 90% confidence intervals from the 50 ensemble simulations for the three hydrologic models. This figure also confirms the presence during the 1931-1932 period of an observed infiltration excess runoff, characterized by a fast baseflow recession curve, which CATHY overestimates for the exceedance percentage range 2-7% due to the higher groundwater contribution from the bedrock aquifer, tRIBS underestimates due to its tendency to simulate steeper recession limbs, and TOPKAPI reproduces quite well.
Water 2019, 11, 2261 16 of 21 Table 7. As Table 6, but for the 1931-1932 simulation period (validation period for tRIBS). For the CATHY model an hypothesis of soil crusting was assumed.

Conclusions
This study examined the potential of direct parameter transferability by assessing the comparative runoff response of three physically based distributed hydrologic models for a small Mediterranean catchment during two contrasting periods. One of the models, tRIBS, had been calibrated and validated against discharge records as part of a prior research study [23]. In the present study, two additional DHMs-TOPKAPI and CATHY-were applied using the same data and parameter sets as the tRIBS model. While the three hydrologic models responded similarly during the wetter year used for the tRIBS calibration, significant differences were found for the much drier two-year period used for the tRIBS validation, in particular for the CATHY model, which produced very low streamflow. To obtain satisfactory results for the CATHY model, an hypothesis of soil crusting was assumed through which the topmost soil layer was assigned a lower saturated hydraulic conductivity [80]. This study provides insights on parameter transferability across distributed hydrologic models characterized by different structures and levels of sophistication in their process descriptions. Although broadly speaking the three physically based models allowed transfer of parameters, it was shown that transferability can break down when simulation conditions are greatly altered.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/11/11/2261/s1, Table S1. Routing parameters for the tRIBS model; Table S2. Routing parameters for the TOPKAPI model; Table S3. Hydraulic geometry parameters for the surface routing module of the CATHY model.