Sensitivity and Performance Analyses of the Distributed Hydrology–Soil–Vegetation Model Using Geomorphons for Landform Mapping

: Landform classiﬁcation is important for representing soil physical properties varying continuously across the landscape and for understanding many hydrological processes in watersheds. Considering it, this study aims to use a geomorphology map (Geomorphons) as an input to a physically based hydrological model (Distributed Hydrology Soil Vegetation Model (DHSVM)) in a mountainous headwater watershed. A sensitivity analysis of ﬁve soil parameters was evaluated for streamﬂow simulation in each Geomorphons feature. As inﬁltration and saturation excess overland ﬂow are important mechanisms for streamﬂow generation in complex terrain watersheds, the model’s input soil parameters were most sensitive in the “slope”, “hollow”, and “valley” features. Thus, the simulated streamﬂow was compared with observed data for calibration and validation. The model performance was satisfactory and equivalent to previous simulations in the same watershed using pedological survey and moisture zone maps. Therefore, the results from this study indicate that a geomorphologically based map is applicable and representative for spatially distributing hydrological parameters in the DHSVM.


Introduction
Hydrological models have been widely used in the literature to study watersheds and to predict future impacts caused by environmental changes [1][2][3][4]. These models are developed according to specific purposes. To accurately represent watershed processes, the choice of the model should rely both on data availability and on its ability to address research goals [5]. Distributed physically based models have a large set of parameters and mathematical equations to quantify water balance [6]. As an example, models such as the DHSVM (Distributed Hydrology Soil Vegetation Model) [7] require information on the spatial distribution of hydrological parameters at a fine scale [8]. These models should be able to create realistic simulations when properly parameterized. However, in many regions, their application is limited due to the required data, so in case of data deficiency conditions, conceptual methods are commonly used [9].
Considering that physically based models have parameters with spatial variability and physical interpretation [10], and sufficient good-quality data are not frequently available, calibration processes can lead to overparameterization or equifinality [11]. These refer to, respectively, when the model has more parameters than can be identified by available data [12] and to when more than one set of parameter values can lead to similar model outputs [13]. In hydrological modelling, two strategies can be used: "bottom-up" and "top-down" approaches. In "bottom-up", field investigations and experiments are conducted at smaller scales [14]. On the other hand, "top-down" is based on process simplifications, such as those proposed by Savenije [15] using "Height Above the Nearest Drainage" (HAND). In the case of distributed modelling, unavailable quality data representing physical properties' heterogeneity can lead to simulations with high predictive uncertainty. Despite the complexity of the representation of hydrological processes, simpler conceptual models often outperform complex ones [15].
Among the limitations for the application of physically based hydrological models, there is the scarcity of field data to develop soil maps through a pedological survey in small watersheds [16]. To overcome such constraints, terrain models' products derived from the digital elevation model (DEM) [17][18][19] have been used for hydrological simulations. In this context, Cuartas et al. [8] and Alvarenga et al. [16] tested an elevation zone map from the HAND model as input for hydrological modelling with the DHSVM on headwater watersheds in Brazil in the Amazon and Atlantic Forest, respectively. The HAND calculates the vertical distance of a DEM cell to the nearest drainage point, creating a normalized elevation map [18]. From this, wetland areas, hill slopes, and plateaus can be identified by a user-defined height threshold, and the result can be used in hydrological modelling as these features exhibit closer characteristics in slope, soil pedogenesis, runoff mechanisms, among others [15].
Different from the HAND, Geomorphons' products consider the terrain geomorphology. They classify the landform from a DEM based on a simple ternary pattern recognition [19]. The advantage of this method over a simple pixel-by-pixel elevation analysis is that it allows classification on larger scales without changing the product's spatial resolution [19]. Previous studies have successfully used Geomorphons for soil pedological identification (digital soil mapping) [20], for correlating landforms with soil texture [21], and for predicting the spatial distribution of soil saturated hydraulic conductivity [22,23]. As these characteristics are related to hydrological processes in watersheds, this tool can potentially represent the spatial distribution of soil properties in hydrological modelling.
Previous studies have used the DHSVM to simulate a small-scale headwater watershed in the Mantiqueira Range using both pedological survey [4] and HAND maps [16]. The Mantiqueira Range is located in southeastern Brazil. This is an important upland region from a hydrological point of view since it hosts the headwater streams of many important and strategic basins of Brazil [24]. Specifically, the Lavrinha Creek Watershed (LCW) is an experimental mountainous watershed and has been a study location for field experiments [23,[25][26][27][28] and Geomorphons applications [21][22][23].
According to Gao et al. [29], the discretization of catchments into small interacting cells in physically based models breaks patterns of landscapes. In this context, a landform map can overcome such constraints. Thus, considering a model previously calibrated in a small-scale largely studied watershed, it enables a "bottom-up" approach in modelling. Accordingly, the hypothesis that the Geomorphons map should improve a fully distributed physically based model performance can be tested by representing landform heterogeneity in soil parametrization. This is an attempt to maintain large-scale patterns in the landscape while modelling a small headwater watershed. For this, the DHSVM was chosen to test this hypothesis based on the fact that it considers the spatial variability of soil parameters and can simulate small headwater watersheds considering spatial heterogeneity.
In this context, the objective of this study is to investigate a fully distributed physically based model (DHSVM) performance when using a landform map from Geomorphons as an input for a headwater watershed in the Lavrinha Creek Watershed. Additionally, to analyze the parametrization effects of landform features for the model's outputs, a sensitivity analysis of soil parameters in relation to the landforms was executed. This study can be particularly important for improving the assessment of hydrological responses in small-scale watersheds using the DHSVM.

Study Area
The Lavrinha Creek Watershed (LCW) is located in the Mantiqueira Range, southeastern Brazil, at the coordinates 22 • 08 28 S and 44 • 26 30 W (Figure 1). Its drainage area is 7.67 km 2 , between altitudes of 1137 and 1732 m above sea level. The watershed is a headwater tributary of the Grande River Basin, which supplies water for a population of approximately 8 million and has a potential installed electricity capacity of 7672 MW [30]. Inserted in the Atlantic Forest biome, the preserved Atlantic Forest occupies 63% of the LCW's area, and pasture mainly occupies the rest. Pasture expansion occurred because the region is traditionally used for dairy and beef cattle production.
Water 2021, 13, x FOR PEER REVIEW 3 of 18 a sensitivity analysis of soil parameters in relation to the landforms was executed. This study can be particularly important for improving the assessment of hydrological responses in small-scale watersheds using the DHSVM.

Study Area
The Lavrinha Creek Watershed (LCW) is located in the Mantiqueira Range, southeastern Brazil, at the coordinates 22°08′28" S and 44°26′30" W ( Figure 1). Its drainage area is 7.67 km 2 , between altitudes of 1137 and 1732 m above sea level. The watershed is a headwater tributary of the Grande River Basin, which supplies water for a population of approximately 8 million and has a potential installed electricity capacity of 7672 MW [30]. Inserted in the Atlantic Forest biome, the preserved Atlantic Forest occupies 63% of the LCW's area, and pasture mainly occupies the rest. Pasture expansion occurred because the region is traditionally used for dairy and beef cattle production.

Data
Weather data were monitored from 2005 to 2010 by a Campbell automatic meteorological station (Figure 1). The average annual rainfall through this period was 2046 mm, averaging maximum and minimum temperatures of 23 and 10 °C, respectively. According to the Köppen climate classification, the LCW has a Cwb climate (temperate subtropical highland climate): mild summers, along with dry and cold winters [28]. The rainy season

Data
Weather data were monitored from 2005 to 2010 by a Campbell automatic meteorological station (Figure 1). The average annual rainfall through this period was 2046 mm, averaging maximum and minimum temperatures of 23 and 10 • C, respectively. According to the Köppen climate classification, the LCW has a Cwb climate (temperate subtropical highland climate): mild summers, along with dry and cold winters [28]. The rainy season occurs between October and April, while the wet season occurs between May and September. From 2006 to 2010, the water level was measured by a Global Water Instrumentation linigraph, model WL16, and a stage-discharge rating curve was fitted from measured streamflow data in a potential model [4]. The forested areas in the watershed have a mean leaf area index (LAI) of 4.05 m 2 m −2 , a basal area of 24.5 m 2 ha −1 , and a canopy height of 9.58 m [28]. The pedological survey of Menezes et al. [26] identified three soil classes at the LCW: Haplic Cambisol (92%), Haplic Gleysol (7%), and Fluvisol (1%). Junqueira Júnior et al. [25] found values of saturated hydraulic conductivity ranging from 0.01 to 32.4 m day −1 , total pore volume between 44% and 76%, field capacity between 19% and 50%, and wilting point between 9% and 25%.

Geomorphons
"Geomorphons" is an algorithm that classifies the landscape according to the patterns of landforms based on the digital elevation model (DEM). This model uses the ternary pattern recognition principle [31], identifying the morphology of the terrain according to zenith and nadir angles within a user-defined look-up distance. The last refers to the number of pixels within the search radius, in which neighbor pixels in eight directions are labeled as 1 if its elevation is above the central pixel plus a flatness threshold (t), −1 if it is below the central pixel less t, and 0 if it is between ±t. From these patterns, landform can be classified into 1 of 10 most commonly recognizable Geomorphons: flat, peak, ridge, shoulder, spur, slope, hollow, footslope, valley, and pit [19]. Further details on Geomorphons calculations can be found in Jasiewicz and Stepinski [19].
Look-up distances and DEM resolutions do considerably change the resulting map. Therefore, an adequate look-up distance can be selected by comparing with a reference soil classification survey [20,21], or by validating the results with experimental data [22,23].

Distributed Hydrology Soil Vegetation Model (DHSVM)
The DHSVM [7] is a physically distributed-based model commonly used in small mountainous watersheds. Among its applications, the DHSVM has been used in tropical areas to study land use and land cover changes [1,4,32] and climate changes [2,3].
The spatial variability of physical parameters is defined by maps with grid cells of the same resolution as the DEM, where soil and vegetation parameters are assigned to each. Based on meteorological data, cell-to-cell water and energy balance are calculated at each time step. The stream channel is represented as a series of connected pixels, where surface and subsurface waters are intercepted and routed through the channel [33].
The DHSVM simulates evapotranspiration, water movement in unsaturated and saturated soils, groundwater recharge, and streamflow [33]. Percolation (q v ) in unsaturated zones is calculated by Darcy's law using Brooks-Corey Equation (1) [33]: where Ks is the soil vertical saturated hydraulic conductivity, φ is the soil porosity, θ r is the residual soil moisture content, and m is the pore size distribution index. Soil transmissivity is calculated assuming that soil lateral saturated hydraulic conductivity decreases exponentially with depth [33]. Version 3.1.2 was used in this study.

DHSVM Input Data
The meteorological data required by the DHSVM are precipitation (m), air temperature ( • C), wind speed (m s −1 ), short-and longwave radiation (W m −2 ), and relative humidity (%). The wind speed at 30 m was estimated using the Bras [34] equation and the shortwave radiation according to Swinbank [35]. Input meteorological data for simulation ranged from 2005 to 2010 at hourly time steps.
The DEM was derived from interpolating contour line maps from the Brazilian Institute of Geography and Statistics (IBGE). Since simulations were performed at an hourly time step, and overland flow is estimated in the model as a ∆x/∆t ratio, being ∆x the spatial resolution and ∆t the time resolution, DEM resolution was defined as 30-m to avoid overestimation of peak flows [36]. A land use map ( Figure 2a) accounts for 63% of the area with the Atlantic Forest and 37% with pasture. The watershed's stream channel was delimited by a threshold of the drainage area of 60,000 m 2 . A soil depth map was created based on the experimental results of Oliveira et al. [37].
The feature map resulting from Geomorphons was used to represent the geomorphology of the LCW. According to Silva et al. [21], there are similarities between the Geomorphons maps and the pedological survey at the LCW [26]. Thus, to analyze drainage porosity, Pinto et al. [23] used Geomorphons maps with a fuzzy logic methodology to interpolate observed field data, determining that a 30-m DEM resolution and a 750 m (25 pixels) look-up distance provided the best representation of that soil physical property. Therefore, for this study, a Geomorphons map was also generated by using a 25-pixel look-up distance from a 30-m spatial resolution DEM.

DHSVM Calibration and Validation
Hourly meteorological data from January 2005 to September 2006 were used to warm up the model, followed by 2 years for calibration (October 2006 to September 2008) and 2 years for validation (October 2008 to September 2010). Alvarenga et al. [4] calibrated the DHSVM for the LCW considering vegetation parameterization as shown in Table 1. As there were no changes in representing vegetation for this study, these values were not altered.  The feature map resulting from Geomorphons was used to represent the geomorphology of the LCW. According to Silva et al. [21], there are similarities between the Geomorphons maps and the pedological survey at the LCW [26]. Thus, to analyze drainage porosity, Pinto et al. [23] used Geomorphons maps with a fuzzy logic methodology to interpolate observed field data, determining that a 30-m DEM resolution and a 750 m (25 pixels) look-up distance provided the best representation of that soil physical property. Therefore, for this study, a Geomorphons map was also generated by using a 25-pixel look-up distance from a 30-m spatial resolution DEM.

DHSVM Calibration and Validation
Hourly meteorological data from January 2005 to September 2006 were used to warm up the model, followed by 2 years for calibration (October 2006 to September 2008) and 2 years for validation (October 2008 to September 2010). Alvarenga et al. [4] calibrated the DHSVM for the LCW considering vegetation parameterization as shown in Table 1. As there were no changes in representing vegetation for this study, these values were not altered.
Interaction among parameters was not considered in this study; the values were set for each soil parameter individually, after sensitivity analysis. Parameter value ranges were obtained by field experiments in the primary literature. The following soil parameters were refined to improve the simulation of streamflow: lateral saturated hydraulic conductivity (LSHC), vertical saturated hydraulic conductivity (VSHC), exponential decay of LSHC (ED), porosity (PR), field capacity (FC), and wilting point (WP). These parameters were selected based on previous sensitivity analyses of the DHSVM [4,8,38] and the availability of field observations [23,25] (Figure 3). Interaction among parameters was not considered in this study; the values were set for each soil parameter individually, after sensitivity analysis. Parameter value ranges were obtained by field experiments in the primary literature. The following soil parameters were refined to improve the simulation of streamflow: lateral saturated hydraulic conductivity (LSHC), vertical saturated hydraulic conductivity (VSHC), exponential decay of LSHC (ED), porosity (PR), field capacity (FC), and wilting point (WP). These parameters were selected based on previous sensitivity analyses of the DHSVM [4,8,38] and the availability of field observations [23,25] (Figure 3).  The statistical performance indexes used as objective functions for calibration were the Nash-Sutcliffe coefficient (NSE) (Equation (2)), the NSE calculated in logarithmic values (lNSE) (Equation (3)), the percentage bias (PBIAS) (Equation (4)), and the root mean square error (RMSE) (Equation (5)).
where n is the total number of data, O i is the observed streamflow at time i, O is the average observed streamflow, and S i is the simulated streamflow at time i. According to Moriasi et al. [39], statistical values of NSE greater than 0.5 and PBIAS of ±15% can be considered a satisfactory model performance.

DHSVM Sensitivity Analysis
In this study, the sensitivity analysis was performed by the single parameter perturbation method, also used in previous studies with the DHSVM [4,8,16,38,40]. To evaluate the parameter sensitivity in the simulated streamflow, soil parameter perturbations were simulated individually for each Geomorphons feature. According to Du et al. [38], the most sensitive soil parameters of the DHSVM on streamflow are LSHC, ED, PR, FC, and WP. Thus, simulations started with the parameter values reported in Table 2, obtained from field measurements. Due to spatial similarities with the soil classes [21], initial soil parameters in Geomorphons features followed the referenced values (Table 2): (I) "valley" feature corresponded to Haplic Gleysol, (II) "pit" feature to Fluvisol, (III) and the remaining features to Haplic Cambisol.
Following the range of soil parameter values observed in situ by Junqueira Junior et al. [25], the LSHC, PR, FC, and WP parameters varied from 0.7 to 1.3 times the reference value. Due to lack of experimental data of ED, this parameter varied from 0.01 to 100 times (Table 2).
For each perturbation, the variation of the simulated streamflow associated with an exceeding probability of 5% (Q5%) and 95% (Q95%) was calculated from the flow duration curve. Thus, the effects of parameter variation on the statistical performance coefficients NSE and lNSE were determined in the simulated period: 2006 to 2010 [41]. Figure 4 shows the landform surface map generated by the Geomorphons algorithm for the LCW. Six different features were identified by landform pattern recognition in DEM data-"ridge" (C1), "spur" (C2), "slope" (C3), "hollow" (C4), "valley" (C5), and "pit" (C6)-which occupy, respectively, 10%, 20%, 32%, 22%, 14%, and 2% of the LCW's drainage area. "Ridge" and "spur" occur mainly in elevated areas of the watershed, and "slope" and "hollow" occur along slopes, while "valley" and "pit" features are characteristic of the valley area of the catchment. Following the range of soil parameter values observed in situ by Junqueira Junior et al. [25], the LSHC, PR, FC, and WP parameters varied from 0.7 to 1.3 times the reference value. Due to lack of experimental data of ED, this parameter varied from 0.01 to 100 times ( Table 2).

Surface Map Geomorphons
For each perturbation, the variation of the simulated streamflow associated with an exceeding probability of 5% (Q5%) and 95% (Q95%) was calculated from the flow duration curve. Thus, the effects of parameter variation on the statistical performance coefficients NSE and lNSE were determined in the simulated period: 2006 to 2010 [41]. Figure 4 shows the landform surface map generated by the Geomorphons algorithm for the LCW. Six different features were identified by landform pattern recognition in DEM data-"ridge" (C1), "spur" (C2), "slope" (C3), "hollow" (C4), "valley" (C5), and "pit" (C6)-which occupy, respectively, 10%, 20%, 32%, 22%, 14%, and 2% of the LCW's drainage area. "Ridge" and "spur" occur mainly in elevated areas of the watershed, and "slope" and "hollow" occur along slopes, while "valley" and "pit" features are characteristic of the valley area of the catchment. Because of soil erosion processes in different landscapes, geomorphology plays an important role in soil's physicochemical properties' spatial distribution [42]. In the LCW, Silva et al. [21] evaluated the spatial soil texture variability of the surface layer (0 to 20 cm) overlapped to the Geomorphons feature map, and found a significant increase in sand and a decrease in clay content in the "valley" and "pit" features. Thus, these areas spatially correspond to, respectively, the Haplic Gleysol and Fluvisol soil classes. Because of soil erosion processes in different landscapes, geomorphology plays an important role in soil's physicochemical properties' spatial distribution [42]. In the LCW, Silva et al. [21] evaluated the spatial soil texture variability of the surface layer (0 to 20 cm) overlapped to the Geomorphons feature map, and found a significant increase in sand and a decrease in clay content in the "valley" and "pit" features. Thus, these areas spatially correspond to, respectively, the Haplic Gleysol and Fluvisol soil classes. Figure 5 shows the bias variation of Q5%, Q95%, NSE, and lNSE of the simulated data when soil parameters (ED, LSHC, PR, FC, and WP) are modified in each Geomorphons feature. The horizontal axis indicates the perturbation of each parameter (Table 2), while the vertical axis indicates the percent variation of Q5%, Q95%, NSE, and lNSE obtained from parameter perturbation responses.

Sensitivity Analysis
High (Q5%) and low (Q95%) streamflow regimes were significantly sensitive for four soil parameters (ED, LSHC, PR, and FC), except WP, because the local climate is mostly humid and the wilting point is rarely reached. Thus, according to the results shown in Figure 5, Q5% and Q95% simulated by the DHSVM were more sensitive to variations of the ED parameter. Responses to these parameter perturbations in Q5% ranged from −12% to 20%, while Q95% ranged from −52% to 21%. The low streamflow regime, indicated by Q95%, responded to LSHC and FC perturbations varying from −5% to 5%, while PR modifications reached a 10% increase in Q95%. These results agree with those of Du et al. [38], which also concluded that Q5% is more sensitive to PR variations since lower soil porosities, in general, reduce infiltration capacity, inducing infiltration excess overland flow.
It is important to highlight that parameter interactions may alter the output in the model. For this analysis, the parameter value was modified after the simulations of each parameter and class, seeking the highest NSE and lNSE values. Figure 6 shows a boxplot of the variations of Q5% and Q95% to the perturbations of the ED, LSHC, PR, FC, and WP parameters for each Geomorphons feature. The x-axis is each Geomorphons feature, and the y-axis represents the variation of Q5% and Q95% values simulated during sensitivity analysis of soil parameters. Boxplot amplitudes closer to 0 indicate a low degree of sensitivity of soil parameters in this feature. It is noticeable that "slope" (C3) and "hollow" (C4) features, followed by "valley" (C5), significantly affect Q5% responses. The Q95% streamflow simulation was more sensitive to parameter perturbations in "slope" (C3), followed by the "hollow" (C4) and "spur" (C2) features.
The "slope" and "hollow" features occur in complex terrain areas of the watershed, and this facilitates the formation of infiltration excess overland flow [43]. Therefore, it explains their influence in both high and low simulated streamflow regimes. "Valley", on the other hand, is characteristic of lowland areas and favors saturation excess overland flow [44]. The "spur" feature occurs in highland areas, important for recharging the water table. For this, soil parameters mostly changed the high streamflow regime. However, this feature was very influential in terms of statistical NSE and lNSE coefficients.
According to Gao et al. [44], mean values and variations of model parameters could be used to find a middle space between complexity and simplifications in the level of segmentation of landscape representativeness. Because several Geomorphons features exhibited a close response to soil parameter perturbation during sensitivity analysis, those features can be merged to reduce uncertainties of parameter estimation when data are unavailable. Water 2021, 13, x FOR PEER REVIEW 10 of 18  table. For this, soil parameters mostly changed the high streamflow regime. However, this feature was very influential in terms of statistical NSE and lNSE coefficients. According to Gao et al. [44], mean values and variations of model parameters could be used to find a middle space between complexity and simplifications in the level of segmentation of landscape representativeness. Because several Geomorphons features exhibited a close response to soil parameter perturbation during sensitivity analysis, those features can be merged to reduce uncertainties of parameter estimation when data are unavailable.

DHSVM Calibration and Validation Using Geomorphons
The final soil parameters in DHSVM simulation that provided the best fitting to the observed daily averaged streamflow in the LCW are presented in Table 3. Thus, simulated and observed hydrographs are shown in Figure 7a,b for calibration and validation, as well as the scattered points around the 1:1 line (Figure 7b,d). DHSVM simulation was able to capture the seasonal streamflow behavior in the LCW. However, it showed a spiky behavior around streamflow peaks, mainly during the validation period. In general, the observed streamflow presented a more rapid recession after the wet season. 0.09 ** 0.09 ** 0.09 ** 0.09 ** 0.12 ** 0.14 ** 0.09 *** 0.09 *** 0.09 *** 0.09 *** 0.12 *** 0.14 *** Vertical conductivity-0.15 * 0.13 * 0.13 * 0.27 * 0.31 * 0.33 * VSHC (10 −4 m s −1 ) 0.14 ** 0.12 ** 0.12 ** 0.26 ** 0.305 ** 0.32 ** 0.13 *** 0.11 *** 0.11 *** 0.25 *** 0.30 *** 0.31 *** The final soil parameters in DHSVM simulation that provided the best fitting to the observed daily averaged streamflow in the LCW are presented in Table 3. Thus, simulated and observed hydrographs are shown in Figure 7a,b for calibration and validation, as well as the scattered points around the 1:1 line (Figure 7b,d). DHSVM simulation was able to capture the seasonal streamflow behavior in the LCW. However, it showed a spiky behavior around streamflow peaks, mainly during the validation period. In general, the observed streamflow presented a more rapid recession after the wet season.   Table 4 presents the results of the performance statistics for the simulated streamflow, indicating a satisfactory fitting for both the calibration and validation with an NSE of, respectively, 0.68 and 0.63 and an lNSE of 0.53 and 0.70. PBIAS suggests an underestimation of 4% and an overestimation of 13%, while RMSE indicates an index error of 0.10 m 3 s −1 for calibration and 0.11 m 3 s −1 for validation. In a previous study at the LCW using a pedological survey map, Alvarenga et al. [4] found an NSE of 0.52 in both calibration and validation periods and an lNSE of 0.06 and 0.58, respectively. Regarding simulation using moisture zone features from the HAND, Alvarenga et al. [16] improved the DHSVM's performance in the LCW, with an NSE in calibration and validation of 0.57 and 0.55 and an lNSE of 0.10 and 0.60, respectively. However, by using Geomorphons features as input in the DHSVM, model performance increased considerably, especially in the lNSE coefficient for calibration, as a result of the detailed spatial variability of a Geomorphons map and previous sensitivity analysis of soil parameters.  Figure 8 shows the observed and simulated flow duration curves (FDCs) from different soil map inputs. The simulated streamflow with Geomorphons features underestimates the observed streamflow below a frequency of 10%. In this case, the simulated streamflow at 5% exceedance probability underestimates the observed streamflow by 0.11 m 3 s −1 . At 90% and 95% exceedance probabilities, the differences between observed and simulated values were 0.02 and 0.03 m 3 s −1 , respectively. In comparison with simulations using pedological and HAND maps, Geomorphons simulation at 90% and 95% probabilities was greater by 0.02 m 3 s −1 , while at 5%, probability was less by 0.06 m 3 s −1 . However, by comparing the maximum simulated streamflow, simulation with Geomorphons overestimated the observed data by 0.04 m 3 s −1 , while with a pedological map overestimated by 0.36 m 3 s −1 , and HAND map by 0.29 m 3 s −1 .
Water 2021, 13, x FOR PEER REVIEW 14 of 18 storage capacity and, consequently, sustained groundwater flow [46]. Therefore, although there is an improvement in performance statistics and peak and base streamflow simulations, the increased spatial variability of soil parameters did not affect the simulated watershed storage capacity. In addition to streamflow analysis, evapotranspiration (ET) was monitored by water balance. However, during the calibration and sensitivity analysis, no significant changes in evaporation estimations were detected. DHSVM simulations in terms of the annual water balance (ET = P − Runoff) indicated that evapotranspiration was 47.4% of the total precipitation (P). Alvarenga et al. [4,16] found a value of 48.1% using both pedological and HAND maps. Mello et al. [28] estimated the ration ET/P to be 49% for the whole catchment, and approximately 50% for the forested portion of the drainage area.

Soil Moisture Spatial Distribution
The spatial variability of simulated soil moisture and water table depth for wet and dry events are presented in the panels of Figure 9. The left column maps correspond to an event in 3 February 2008, at 15 h local time, after 109 mm of antecedent precipitation during the previous 12 h, and 260 mm throughout the previous 7 days (wet event). As an example of the conditions in the dry season, the right column of Figure 9 shows the maps for 11 August 2007, at 15 h, with no recorded precipitation in the previous 14 days (dry event). In addition, the black contour line of Figure 9 outlines the Geomorphons features "hollow", "valley", and "pit" also depicted in  Beckers and Alila [45] highlighted the DHSVM's limitations in capturing preferential flow path in soil layers that contribute to subsurface streamflow. The ED parameter is used to account for the saturated hydraulic conductivity difference in relation to soil depth. However, while increasing this parameter, the magnitude of the simulated peak streamflow increases, and recession is accelerated, causing underestimation of baseflows. Other studies have also reported a similar behavior in DHSVM simulations [4,8]. Despite this, simulations with the DHSVM improved using Geomorphons features when compared with previous studies [4,16] in terms of simulated peak and base streamflow ( Figure 8). These results indicate that an increasing spatial variability of soil parameters can be beneficial for simulating extreme events in the watershed.
For exceeding probabilities between 20% and 70%, DHSVM simulations with Geomorphons, pedological, and HAND maps exhibited a similar behavior: all of them produced FDCs with a slope greater than the observed. For this range of probabilities, the FDC is determined by percolation rates: a mild slope of this segment indicates a larger soil storage capacity and, consequently, sustained groundwater flow [46]. Therefore, although there is an improvement in performance statistics and peak and base streamflow simulations, the increased spatial variability of soil parameters did not affect the simulated watershed storage capacity.
In addition to streamflow analysis, evapotranspiration (ET) was monitored by water balance. However, during the calibration and sensitivity analysis, no significant changes in evaporation estimations were detected. DHSVM simulations in terms of the annual water balance (ET = P − Runoff) indicated that evapotranspiration was 47.4% of the total precipitation (P). Alvarenga et al. [4,16] found a value of 48.1% using both pedological and HAND maps. Mello et al. [28] estimated the ration ET/P to be 49% for the whole catchment, and approximately 50% for the forested portion of the drainage area.

Soil Moisture Spatial Distribution
The spatial variability of simulated soil moisture and water table depth for wet and dry events are presented in the panels of Figure 9. The left column maps correspond to an event in 3 February 2008, at 15 h local time, after 109 mm of antecedent precipitation during the previous 12 h, and 260 mm throughout the previous 7 days (wet event). As an example of the conditions in the dry season, the right column of Figure 9 shows the maps for 11 August 2007, at 15 h, with no recorded precipitation in the previous 14 days (dry event). In addition, the black contour line of Figure 9 outlines the Geomorphons features "hollow", "valley", and "pit" also depicted in Figure 4.
Water 2021, 13, x FOR PEER REVIEW 15 of 18 the saturation of top soil layers [4,32]. Soil moisture was greater near the "pit" feature ( Figure 4) for both the wet and dry events, which explains the low sensitivity of soil parameters in this feature. The shallow soil depth and smaller LSHC parameter values for the "ridge" and "spur" features explain the higher soil moisture and water table depth during a wet event.
Soil moisture and landforms are important characteristics for studying not only the water regime in watersheds. Soil moisture modelling in physically based models can be used for predicting empirical thresholds associated with triggering shallow landslides [48,49], as the precedent soil moisture content, associated with rainfall intensity, helps improve the predictive accuracy of shallow landslides [49]. Therefore, the relationships between geomorphology and soil moisture content in physically based modelling should be further explored.

Conclusions
The streamflow and soil moisture simulations of a headwater mountainous watershed in the DHSVM using Geomorphons map were satisfactory. The patterns in geomorphology were accounted for in the model, representing heterogeneity in soil properties. Results showed a higher soil moisture and a shallower water table depth close to the stream channel due to the saturation of top soil layers [4,32]. Soil moisture was greater near the "pit" feature ( Figure 4) for both the wet and dry events, which explains the low sensitivity of soil parameters in this feature. The shallow soil depth and smaller LSHC parameter values for the "ridge" and "spur" features explain the higher soil moisture and water table depth during a wet event.
Soil moisture and landforms are important characteristics for studying not only the water regime in watersheds. Soil moisture modelling in physically based models can be used for predicting empirical thresholds associated with triggering shallow landslides [48,49], as the precedent soil moisture content, associated with rainfall intensity, helps improve the predictive accuracy of shallow landslides [49]. Therefore, the relationships between geomorphology and soil moisture content in physically based modelling should be further explored.

Conclusions
The streamflow and soil moisture simulations of a headwater mountainous watershed in the DHSVM using Geomorphons map were satisfactory. The patterns in geomorphology were accounted for in the model, representing heterogeneity in soil properties.
Sensitivity analysis indicated a strong influence of the ED, LSHC, and PR parameters on the NSE coefficient. Streamflows at 5% (Q5%) and 95% (Q95%) probabilities were sensitive to the ED, PR, LSHC, and FC parameters. The Geomorphons feature with the greatest influence on the sensitivity parameters was "slope", followed by "hollow" and "valley", which are related to the streamflow generation processes of infiltration and saturation excess overland flow, which are important mechanisms of streamflow in complex terrain watersheds.
The DHSVM performance statistics using a Geomorphons map in terms of streamflow improved from previous studies using pedological and HAND maps. Thus, the simulated soil moisture content is in accordance with experimental data. Based on these results in the LCW, a geomorphological map is an accurate tool for spatially distributing soil parameters in distributed models and for simulating the spatial variability of runoff generation processes in the catchment.
Therefore, the use of Geomorphons maps in a complex terrain headwater watershed offers an alternative for efficiently applying the DHSVM as a water resource management tool. However, the large number of Geomorphons features can lead to over-parametrization if there are no sufficiently available data for calibration. As an alternative, future studies may analyze model performance by merging similar Geomorphons features to reduce parameterization uncertainties in a "top-down" approach.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.