Improving Hillslope Link Model Performance from Non-Linear Representation of Natural and Artiﬁcially Drained Subsurface Flows

: This study evaluates the potential for a newly proposed non-linear subsurface ﬂux equation to improve the performance of the hydrological Hillslope Link Model (HLM). The equation contains parameters that are functionally related to the hillslope steepness and the presence of tile drainage. As a result, the equation provides better representation of hydrograph recession curves, hydrograph timing, and total runoff volume. The authors explore the new parameterization’s potential by comparing a set of diagnostic and prognostic setups in HLM. In the diagnostic approach, they conﬁgure 12 different scenarios with spatially uniform parameters over the state of Iowa. In the prognostic case, they use information from topographical maps and known locations of tile drainage to distribute parameter values. To assess performance improvements, they compare simulation results to streamﬂow observations during a 17-year period (2002–2018) at 140 U.S. Geological Survey (USGS) gauging stations. The operational setup of the HLM model used at the Iowa Flood Center (IFC) serves as a benchmark to quantify the overall improvement of the model. In particular, the new equation provides better representation of recession curves and the total streamﬂow volumes. However, when comparing the diagnostic and prognostic setups, the authors found discrepancies in the spatial distribution of hillslope scale parameters. The results suggest that more work is required when using maps of physical attributes to parameterize hydrological models. The ﬁndings also demonstrate that the diagnostic approach is a useful strategy to evaluate models and assess changes in their formulations.


Introduction
Flood forecasts that are calculated using regional distributed hydrological models are becoming more common and relevant because they also provide information about internal watershed processes in large domains, along with predicted hydrographs for all streams in the river network. These forecasts are expected to be accurate at the region's ungauged watersheds [1] as a consequence of appropriate spatial representation of processes and parameters in the model. Current hydrological models correctly identify many aspects of the streamflow hydrographs, thereby providing acceptable forecasts. However, they still struggle to reproduce the hydrograph recession. According to [2], modelers need to pay more attention to storm runoff's slow flow, which is a crucial component of the recession. For regional models, recession becomes more challenging because its non-linearity increases with the spatial scale [3][4][5]. Landscape properties such as topography, soil, and the stream network seem to be involved in the recession variability [6][7][8]. Additionally, human landscape and land-use els, recession becomes more challenging because its non-linearity increases with the spatial scale [3][4][5]. Landscape properties such as topography, soil, and the stream network seem to be involved in the recession variability [6][7][8]. Additionally, human landscape and land-use interventions, such as tile drainage, can help to restore river health [9]. Nevertheless, these interventions also have the potential to alter streamflow and its recession [10,11]. A good representation of the recession may influence the estimation of the flood duration, and flood occurrence during successive events.

Issues with the Hillslope Link Model (HLM) in Iowa
The Iowa Flood Center (IFC) produces flood forecasts for the state of Iowa using the Hillslope Link Model (HLM) [12][13][14][15]. The operational HLM represents the hillslope subsurface flux, using a linear-reservoir equation. According to [15], the current HLM configuration accurately estimates peak flows with an overall acceptable performance in Iowa. However, the model has some limitations in capturing the hydrograph recessions and the total runoff volume at some locations, as reported in [14]. The discrepancies between simulated and observed recessions are more common in watersheds with tiling drainage. Sample streamflow simulation results, using the IFC HLM operational model for three Iowa watersheds are presented in Figure 1a-c (in red). The model performance is described in terms of Kling Gupta Efficiency [16], which is defined as KGE = 1 − √( − 1) 2 + ( − 1) 2 + ( − 1) 2 , where , , and denote correlation, the ratio of standard deviation ( / ), and the ratio of the mean ( / ) between simulated and observed streamflow, respectively. The model's limitations are most evident in the watersheds located in the north and west regions of Iowa, where the model has low performance in terms of the Kling Gupta Efficiency (KGE) index ( Figure 1d). We associate the model's poor performance in the region of north-central Iowa, known as the Des Moines Lobe, with the widespread use of artificial subsurface drainage (known as tile drains) in the region [11]. To address these issues, the authors of [17] developed a subsurface non-linear equation that can represent subsurface flow from hillslopes with different steepness and soil conductivities, as well as the presence of tile drainage. The blue lines in Figure 1a-c show the resulting hydrographs using the non-linear equation with parameters corresponding to no tile and a steepness of 2% [17]. Compared with the linear equation of the operational HLM, the non-linear equation tends to improve the total streamflow volume and the simulated recession shapes. However, we still observe discrepancies (Figure 1a,b) that are To address these issues, the authors of [17] developed a subsurface non-linear equation that can represent subsurface flow from hillslopes with different steepness and soil conductivities, as well as the presence of tile drainage. The blue lines in Figure 1a-c show the resulting hydrographs using the non-linear equation with parameters corresponding to no tile and a steepness of 2% [17]. Compared with the linear equation of the operational HLM, the non-linear equation tends to improve the total streamflow volume and the simulated recession shapes. However, we still observe discrepancies (Figure 1a,b) that are attributed to issues with the parameter values and spatial representation of processes. The model proposed in [15] serves as motivation for this study, to evaluate the performance of such model modifications.

The Diagnostic-Prognostic Approach
According to the authors of [18], the development of a hydrological model is subject to the hypothesis-testing process. This process evaluates, rejects, and replaces model components. We performed a diagnostic-prognostic analysis of the model at 140 USGS gauges in Iowa to test the utility of the non-linear equation to represent the hillslope subsurface flux. In this case, we adapted the diagnostic-prognostic approach developed in studies on evapotranspiration [19][20][21]. Our diagnostic setups have simplified, spatially uniform parameter values, while the prognostic scenarios use maps to determine parameter values. The diagnostic-prognostic approach offers complementary information about the model [22] and the required independence to perform model comparisons [23].
According to [15], an insightful way to improve models starts with model performance verification, followed by structure modification. We expanded on this approach by using the diagnostic-prognostic analysis to add tools to verify the model's processes and required parameters. Our objective is to identify the non-linear model parameters for Iowa, their uncertainties, and the model limitations.
The paper is organized as follows: we first describe the HLM model and the equations governing the hillslope processes, including the linear equation and non-linear equations to represent subsurface and drainage tile fluxes. Next, we describe the diagnostic and prognostic setups. Then, we compare the diagnostic and prognostic approach results using 140 USGS stations and we analyze the parameters' influence on the model performance. Finally, we provide conclusions based on our experiment results.

Model Description
The Hillslope Link Model (HLM) represents the hydrological processes at the hillslope scale (Figure 2a,b) and routes the streamflow through the channel network (Figure 2c). At the hillslopes, HLM has three storages, ponded surface (S p [m]), topsoil (S T [m]), and subsurface storage (S s [m]). The water from the ponded storage can either infiltrate the topsoil (q pT [m · min −1 ]) or flow as runoff to the channel link (q pL [m · min −1 ]). The water in the topsoil percolates (q Ts [m · min −1 ]) to the soil storage. Finally, the water in the soil storage seeps into the channel link as subsurface runoff (q sL [m · min −1 ]). Evaporation occurs from the three storages as a removal of volume from the model. Once in the river network, HLM transports the channel water (q [m 3 · s −1 ]) downstream. A detailed description of the hillslope and stream routing process can be found in [11,14].
The surface runoff, infiltration, and percolation rates are linked through the reference speed v r and the shape of the hillslope. Each hillslope has a parameter k 2 [min −1 ] (Equation (1)) that depends on the hillslope link length (L i [m]) and area (A h m 2 ), along with the reference velocity v r . The parameter k 2 is the inverse of the runoff residence time in the hillslope. The runoff q pL and the infiltration q pT are linked to k 2 through Equations (2) and (3), respectively. Moreover, the infiltration also depends on the topsoil depth (T l ) assumed to be equal to 0.1 [m]. Additionally, the percolation rate q Ts is computed as a proportion of k 2 , expressed by k i . Usually, k i is 2% of k 2 ; however, its value may change depending on the soil and topographical properties.
q Ts = k 2 · S T · k i (4)  The surface runoff, infiltration, and percolation rates are linked through the reference speed and the shape of the hillslope. Each hillslope has a parameter 2 [min −1 ] (Equation (1)) that depends on the hillslope link length ( [m]) and area ( ℎ [m 2 ]), along with the reference velocity . The parameter 2 is the inverse of the runoff residence time in the hillslope. The runoff and the infiltration are linked to 2 through Equations (2) and (3), respectively. Moreover, the infiltration also depends on the topsoil depth ( ) assumed to be equal to 0.1 [m]. Additionally, the percolation rate is computed as a proportion of 2 , expressed by . Usually, is 2% of 2 ; however, its value may change depending on the soil and topographical properties.
The current HLM setup represents the subsurface flux to the channels ( [m ⋅ min −1 ]) with a linear equation (red line on Figure 2d). The equation releases water to the channel at a rate , when is greater than the no-flow threshold ( ), as follows, Ref. [17] developed a set of parameterizations for ordinary differential equations that adds a non-linear component to Equation (5)   The current HLM setup represents the subsurface flux to the channels (q sL [m · min −1 ]) with a linear equation (red line on Figure 2d). The equation releases water to the channel at a rate m, when S s is greater than the no-flow threshold (S o ), as follows, Ref. [17] developed a set of parameterizations for ordinary differential equations that adds a non-linear component to Equation (5) when S s is above threshold storage. The following exponential equation (continuous line on Figure 2d) is added to Equation (5) if S s is greater than the activation threshold β [m], where α is a parameter that depends on the hillslope properties, such as its steepness and the soil conductivity. In [17], the authors also developed an exponential equation that applies when the hillslope has tiles. The following equation (dashed line on Figure 2d) is added when S s is greater than the tile relative depth D d [m], In the described scheme, the subsurface flux becomes a set of equations that HLM activates, depending on the value of S s relative to the thresholds S o , β, and D d . The segmented subsurface runoff is as follows, The relative tile depth (D d ) is independent of β, so either could be larger depending on the tile configuration and the hillslope properties. Moreover, if there are no tiles, Equation (8) is limited to its two first expressions. More details on the subsurface equation development can be found in [17].

Model Setup and Data
We used both diagnostic and prognostic approaches to test the performance using the non-linear equation. We used the river network for the state of Iowa derived from a DEM of 90 m and decomposed into about 420,000 individual hillslopes, following the approach presented in [11]. The precipitation force corresponds to hourly Stage IV QPEs [24]. We forced the evapotranspiration using the mean annual monthly values derived from MODIS [25] for the region. Equation (8) offers a formulation for the subsurface flux that we want to validate for the Iowa domain. In this process, we can fix parameters uniformly over the space or distribute them spatially. A uniform setup assumes that each hillslope in the region uses the same model parameters, while a distributed setup assumes parameter variability as a function of landscape properties. Neither approach is without error because the parameters are only approximate and could depend upon unknown factors that are variable in space. The fixed setup is unrealistic, and the distributed setup may be subject to spatial errors. However, both approaches are complementary. Fixed setups could help assess the ability of Equation (8) to improve the accuracy of the simulated streamflow fluctuations. In contrast, a distributed setup helps to validate the parameter description given by the map(s). Considering this, we used both approaches to validate the new q sL equation and to explore the limits of the so-called predefined setups. In the distribute parameters case, we use the steepness of the hillslopes ( Figure 3a) and the tiles localization according to the Iowa Department of Natural Resources (DNR) (Figure 3b).

+ + >
The relative tile depth ( ) is independent of , so either could be larger depending on the tile configuration and the hillslope properties. Moreover, if there are no tiles, Equation (8) is limited to its two first expressions. More details on the subsurface equation development can be found in [17].

Model Setup and Data
We used both diagnostic and prognostic approaches to test the performance using the non-linear equation. We used the river network for the state of Iowa derived from a DEM of 90 m and decomposed into about 420,000 individual hillslopes, following the approach presented in [11]. The precipitation force corresponds to hourly Stage IV QPEs [24]. We forced the evapotranspiration using the mean annual monthly values derived from MODIS [25] for the region.
Equation (8) offers a formulation for the subsurface flux that we want to validate for the Iowa domain. In this process, we can fix parameters uniformly over the space or distribute them spatially. A uniform setup assumes that each hillslope in the region uses the same model parameters, while a distributed setup assumes parameter variability as a function of landscape properties. Neither approach is without error because the parameters are only approximate and could depend upon unknown factors that are variable in space. The fixed setup is unrealistic, and the distributed setup may be subject to spatial errors. However, both approaches are complementary. Fixed setups could help assess the ability of Equation (8) to improve the accuracy of the simulated streamflow fluctuations. In contrast, a distributed setup helps to validate the parameter description given by the map(s). Considering this, we used both approaches to validate the new equation and to explore the limits of the so-called predefined setups. In the distribute parameters case, we use the steepness of the hillslopes (  The formulation of Equation (8) relies on the percolation rate because the non-linear formulation depends upon the amount of water in the subsurface storage. The described dependence increases the relevance of the percolation parameter (k i ). The distribution of k i can be derived from maps of the soil profile properties. However, using an additional map may increase the errors affecting the comparison of both setups. For this reason, we choose to fix three different percolation rates for the diagnostic and prognostic setups ( Figure 4c).
Moreover, we used the same values for S 0 , β, and D d in both setups. In [17] S 0 oscillates between 1.4 and 1.55 in function of the slope. Nevertheless, additional features of the hillslope such as the bedrock depth, and soil type determine the value of S 0 . To avoid adding additional uncertainty sources, we fixed S 0 = 1.48. On the other hand, for β and D d we used the values described by [17] of 1.67 and 1.635, respectively. In Table 1 we summarize the described diagnostic and prognostic setups. The model validation consists of comparing fixed (diagnostic) and distributed (prognostic) HLM setups ( Figure 4). The diagnostic setup ( Figure 4a) shows how different formulations could significantly improve the model across the region. On the other hand, the prognostic setups (Figure 4b) show the improvements and limitations derived from the application of "known" spatial variables. The formulation of Equation (8) relies on the percolation rate because the non-linear formulation depends upon the amount of water in the subsurface storage. The described dependence increases the relevance of the percolation parameter ( ). The distribution of can be derived from maps of the soil profile properties. However, using an additional map may increase the errors affecting the comparison of both setups. For this reason, we choose to fix three different percolation rates for the diagnostic and prognostic setups ( Figure 4c).
Moreover, we used the same values for 0 , , and in both setups. In [17] 0 oscillates between 1.4 and 1.55 in function of the slope. Nevertheless, additional features of the hillslope such as the bedrock depth, and soil type determine the value of 0 . To avoid adding additional uncertainty sources, we fixed 0 = 1.48. On the other hand, for β and we used the values described by [17] of 1.67 and 1.635, respectively. In Table 1 we summarize the described diagnostic and prognostic setups.   In the diagnostic setup (Figure 4a), we created four base parametrizations, using Equation (8) for the Iowa domain. The parametrizations range from flat hillslopes (light blue line on Figure 4a) to steep or tiled hillslopes (red line on Figure 4a). In the tiled case, we used a steepness of 2% under the assumption that tiles are usually installed in flat terrains. Then, we combined the four parameterizations and the three k i rates to obtain 12 diagnostic (D) scenarios (D1 to D12 in Figure 4c). In the scenarios, D1 to D4 use k i = 0.02; D5 to D8 use k i = 0.03; and D9 to D12 use k i = 0.04. A summary of the diagnostic scenarios is provided in Table 1.

Prognostic Setups
In the prognostic setup, we distributed parameter values in function of the hillslopes steepness and the Iowa DNR map describing tile presence (Figure 4b). According to [17] the coefficient α of Equation (6) changes in function of the hillslope steepness (γ h ) with a linear equation. Using the following equation, we assigned the parameter α to each hillslope, obtaining variable q s L curves that oscillate between the blue bands shown in Figure 4b: Additionally, we include Equation (8) for tiled terrain. To do it, we assigned tiles to the hillslopes based in the map in Figure 4b. In the tile drainage case, we changed c in function of the slope using the following equation: Then, we developed the prognostic (P) scenarios P1, P2, and P3 using the distributed parameters α and c with the percolation rates k i of 0.02, 0.03 and 0.04, respectively (distributed setups in Figure 4c). A summary of the prognostic scenarios is provided in Table 1.

Insights from a Diagnostic-Prognostic Approach
The diagnostic and prognostic setups produced significant differences between the model outputs. In Figure 5, we present the simulated hydrographs at three watersheds simulated by the diagnostic scenario D4 in blue and the prognostic scenario P1 in red. In this case, the diagnostic setup assumes that all the hillslopes have tiles or are steep. On the other hand, the prognostic setup assumes there to be tiles only at some hillslopes and that the parameter α of Equation (6) varies with the steepness. In these three cases, the diagnostic (or fixed) setup produces a longer recession curve than the one obtained by the prognostic setup. The diagnostic case has a better match to the Iowa River at Tama (Figure 5b), while the prognostic setup exhibits a better match to the White Breast Creek (Figure 5a) and at the Cedar River ( Figure 5c). Figure 5 gives a brief description of the expected differences between the setups. Additionally, it shows that Equation (8) can improve the streamflow representation, given the correct set of parameters that are obtained.
According to Figure 5, the non-linear model can produce a good representation of the hydrograph falling limb and early recession, depending on the parameters. Considering the described sensitivity, we compare the event-based KGE of the non-linear setups and the linear model ( Figure 6). The KGE equation summarizes the correlation (γ), the mean ratio (µ), and the deviation ratio (σ). Our results suggest that the KGE performance depends heavily on the percolation rate (k i ). With k i = 0.02 (first row of Figure 6), all the non-linear setups tend to improve the linear model, with a significant performance decrease in some events. Conversely, values of k i equal to 0.03 and 0.04 do not exhibit a significant KGE change (second and third rows of Figure 6). Cases such as D5 and D11 exhibited a performance such as that obtained by the linear model. Other cases, such as D9, tended toward a general decrease in performance. D6, D8, and P2 exhibited a slight performance increase. The described results highlight the relevance of the percolation rate and the subsurface parameters. The comparison with the linear model shows that Equation (8) can significantly improve the model performance, depending on the parameters. 5a) and at the Cedar River (Figure 5c). Figure 5 gives a brief description of the expected differences between the setups. Additionally, it shows that Equation (8) can improve the streamflow representation, given the correct set of parameters that are obtained. According to Figure 5, the non-linear model can produce a good representation of the hydrograph falling limb and early recession, depending on the parameters. Considering the described sensitivity, we compare the event-based KGE of the non-linear setups and the linear model ( Figure 6). The KGE equation summarizes the correlation ( ), the mean ratio ( ), and the deviation ratio ( ). Our results suggest that the KGE performance depends heavily on the percolation rate ( ). With = 0.02 (first row of Figure 6), all the non-linear setups tend to improve the linear model, with a significant performance decrease in some events. Conversely, values of equal to 0.03 and 0.04 do not exhibit a significant KGE change (second and third rows of Figure 6). Cases such as D5 and D11 exhibited a performance such as that obtained by the linear model. Other cases, such as D9, tended toward a general decrease in performance. D6, D8, and P2 exhibited a slight performance increase. The described results highlight the relevance of the percolation rate and the subsurface parameters. The comparison with the linear model shows that Equation (8) can significantly improve the model performance, depending on the parameters.  Differences among the scenarios are highlighted by comparing the performance gauge by gauge. First, we choose the diagnostic (D) and prognostic (P) setup with the best performance at each gauge. For this, we used the KGE to select the setup that outperformed the others for most of the events. In Figure 7, we present the KGE distribution and the percentage of time each scenario was chosen. We found similarities between the diagnostic and prognostic chosen setups when grouped by the percolation rate values ( ). D4 and P1 ( = 0.02) have a similar KGE distribution, as do D8 and P2 ( = 0.03) and the group that includes D9, D11, D12, and P3 ( = 0.04). The similarities among the described groups highlights the relevance of . Moreover, some differences also highlight the rele- Differences among the scenarios are highlighted by comparing the performance gauge by gauge. First, we choose the diagnostic (D) and prognostic (P) setup with the best performance at each gauge. For this, we used the KGE to select the setup that outperformed the others for most of the events. In Figure 7, we present the KGE distribution and the percentage of time each scenario was chosen. We found similarities between the diagnostic and prognostic chosen setups when grouped by the percolation rate values (k i ). D4 and P1 (k i = 0.02) have a similar KGE distribution, as do D8 and P2 (k i = 0.03) and the group that includes D9, D11, D12, and P3 (k i = 0.04). The similarities among the described groups highlights the relevance of k i . Moreover, some differences also highlight the relevance of Equation (8) parameters.
fixed percolation rate. Columns correspond to the four fixed equations. The color bar shows the percentage of events that fall at each bin of the 2D histogram.
Differences among the scenarios are highlighted by comparing the performance gauge by gauge. First, we choose the diagnostic (D) and prognostic (P) setup with the best performance at each gauge. For this, we used the KGE to select the setup that outperformed the others for most of the events. In Figure 7, we present the KGE distribution and the percentage of time each scenario was chosen. We found similarities between the diagnostic and prognostic chosen setups when grouped by the percolation rate values ( ). D4 and P1 ( = 0.02) have a similar KGE distribution, as do D8 and P2 ( = 0.03) and the group that includes D9, D11, D12, and P3 ( = 0.04). The similarities among the described groups highlights the relevance of . Moreover, some differences also highlight the relevance of Equation (8)   The results presented in Figure 7 follow a spatial distribution. Figure 8 shows each USGS gauge colored by the diagnostic (Figure 8a) and prognostic ( Figure 8b) setups with the best performance. In both cases, the percolation rate defines the spatial distribution. We can identify how the chosen setups ( Figure 8) follow the Iowa landforms to some extent in the diagnostic case (see Figure 3a). Scenario D12 is recurrent over the Des Moines Lobe and the Northwest Iowa Plain. D9 recurs over the Missouri River Alluvial and Loess Hills landforms. We found that D4 dominates over the Southern Iowa Drift area. In the The results presented in Figure 7 follow a spatial distribution. Figure 8 shows each USGS gauge colored by the diagnostic (Figure 8a) and prognostic ( Figure 8b) setups with the best performance. In both cases, the percolation rate defines the spatial distribution. We can identify how the chosen setups ( Figure 8) follow the Iowa landforms to some extent in the diagnostic case (see Figure 3a). Scenario D12 is recurrent over the Des Moines Lobe and the Northwest Iowa Plain. D9 recurs over the Missouri River Alluvial and Loess Hills landforms. We found that D4 dominates over the Southern Iowa Drift area. In the remaining regions, we see a mix of scenarios. The spatial distribution is similar among the chosen prognostic scenarios (Figure 8b) and seems to be highly influenced by the percolation rates, represented here by tones of blue (k i = 0.02), red (k i = 0.03), and green (k i = 0.04).  According to Figure 8, the chosen diagnostic and prognostic scenarios share percolation rates. However, differences exist in the spatial performance improvement distribution ( Figure 9). Figure 9a,b show the diagnostic and prognostic scenarios of KGE improvement with respect to the linear model. With only two cases of negative KGE differences (red dots on Figure 9a), the diagnostic scenarios outperform the linear model at almost all the USGS gauges. Alternatively, in the prognostic case (Figure 9b), the count of negative KGE differences increases to 13, while the number of gauges decreases in cases where the improvement is more significant than 0.1 (yellow). We attribute the decrease in prognostic case performance to the parameter's spatial distribution. According to Figure 8, the chosen diagnostic and prognostic scenarios share percolation rates. However, differences exist in the spatial performance improvement distribution ( Figure 9). Figure 9a,b show the diagnostic and prognostic scenarios of KGE improvement with respect to the linear model. With only two cases of negative KGE differences (red dots on Figure 9a), the diagnostic scenarios outperform the linear model at almost all the USGS gauges. Alternatively, in the prognostic case (Figure 9b), the count of negative KGE differences increases to 13, while the number of gauges decreases in cases where the improvement is more significant than 0.1 (yellow). We attribute the decrease in prognostic case performance to the parameter's spatial distribution.
According to Figure 8, the chosen diagnostic and prognostic scenarios share percolation rates. However, differences exist in the spatial performance improvement distribution ( Figure 9). Figure 9a,b show the diagnostic and prognostic scenarios of KGE improvement with respect to the linear model. With only two cases of negative KGE differences (red dots on Figure 9a), the diagnostic scenarios outperform the linear model at almost all the USGS gauges. Alternatively, in the prognostic case (Figure 9b), the count of negative KGE differences increases to 13, while the number of gauges decreases in cases where the improvement is more significant than 0.1 (yellow). We attribute the decrease in prognostic case performance to the parameter's spatial distribution.  The prognostic scenario performance decrease occurs mostly over the east and west regions of Iowa. The most significant decrease is observed for the Northwest Iowa Plains landform (Figure 9b). In this region, the chosen diagnostic setups were D12 and D9 (Figure 8a), suggesting a mix between tiled terrain and flat hillslopes. Over the Southern Iowan Drift landform area, the k i value is the same for the diagnostic and prognostic scenarios. However, the prognostic scenario performance declines at several stations in this region. On the other hand, the Iowa Surface region exhibits more k i discrepancies between both scenarios, as well as a higher number of performance differences. The described results suggest a level of heterogeneity in the parameters shown by the diagnostic and prognostic scenarios. This heterogeneity creates difficulties when choosing the most adequate regional parameterization for the model, regardless of whether it is fixed (diagnostic) or distributed (prognostic). To address this issue, we compare the KGE (upper diagonal in Figure 10) and the mean ratio µ (lower diagonal in Figure 10) of the chosen scenarios. We estimate µ as the ratio between the simulated (µ sim ) and observed (µ obs ) flows with values near indicating a perfect match. According to Figure 10, the KGE and mean ratio of scenarios D4 and P1 outperform almost all the scenarios. Additionally, both scenarios have the highest percentage of events with KGE values above 0.4 (blue bars in Figure 10 histograms). Compared with the linear model, the D4 and P1 mean ratio correction is significant. In both plots (Linear-D4 and Linear-P1), there are almost no events where the linear setup outperformed D4 and P1. The scenarios D4 and P1 have the same k i (0.02) value; however, their subsurface parameters are different.
The parameters of D4 are fixed for all the domains following line 4 of Figure 4a. This parameterization represents highly conductive soils or the presence of tiles. On the other hand, P1 parameters follow the hillslope steepness with Equation (9), and the presence of tiles described by the map in Figure 4b. The described differences in the parameters seem to develop slight dissimilarities in performance. According to panel D4-P1 in Figure 10, the KGE performance is similar in both, although D4 has a better performance in some events. Moreover, the panel P1-D4 shows that the mean ratio description of both setups is similar. Considering that D4 assumes tiles everywhere, our results suggest a high presence of tile-like signatures.
scenarios. Additionally, both scenarios have the highest percentage of events with KGE values above 0.4 (blue bars in Figure 10 histograms). Compared with the linear model, the D4 and P1 mean ratio correction is significant. In both plots (Linear-D4 and Linear-P1), there are almost no events where the linear setup outperformed D4 and P1. The scenarios D4 and P1 have the same (0.02) value; however, their subsurface parameters are different.

Extended Metrics
According to the diagnostic and prognostic KGE comparisons, the performance differences between the two scenarios are relatively small. However, the KGE is subject to three parameters that do not necessarily reflect all the relevant changes in the simulated streamflows. With this in mind, we also compared the NSE (Nash Sutcliffe efficiency), the hit rate, and the lags (Figure 11a-c, respectively). The NSE contrasts the simulated data prediction skill with the mean value of the observations. An NSE of below 0 indicates that the mean value performs better than the model, and an NSE of 1 indicates a perfect simulation. The hits rate is the percentage of time that is shared by observations and simulations during floods. A hit rate of zero corresponds to missing all the floods, and a hit rate of one corresponds to a perfect match. The lag is a measure of the displacements applied over the simulated data to maximize the correlation. We made hourly displacements from −48 h to 48 h. Negative lag values correspond to cases in which the simulated data exhibit responses earlier than the observed, and positive values correspond to the opposite behavior. A simulated series with good performance must have lags near zero.
Results from Figure 11 show that D4 and P1 have some similarities and some relevant differences. Regarding the NSE (Figure 11a), D4 has a slightly better score. Regarding the hits rate (Figure 11b), it is hard to tell which scenario has a better performance. Furthermore, P1 has a higher fraction of hit rates approaching one, but it also has higher frequencies at some lower intervals. The number of lags is also similar (Figure 11c). Nevertheless, P1 tends toward negative lag values more than D4 does, representing more frequent early peak estimations. simulations during floods. A hit rate of zero corresponds to missing all the floods, and a hit rate of one corresponds to a perfect match. The lag is a measure of the displacements applied over the simulated data to maximize the correlation. We made hourly displacements from −48 h to 48 h. Negative lag values correspond to cases in which the simulated data exhibit responses earlier than the observed, and positive values correspond to the opposite behavior. A simulated series with good performance must have lags near zero. Results from Figure 11 show that D4 and P1 have some similarities and some relevant differences. Regarding the NSE (Figure 11a), D4 has a slightly better score. Regarding the hits rate (Figure 11b), it is hard to tell which scenario has a better performance. Furthermore, P1 has a higher fraction of hit rates approaching one, but it also has higher frequencies at some lower intervals. The number of lags is also similar (Figure 11c). Nevertheless, P1 tends toward negative lag values more than D4 does, representing more frequent early peak estimations.
In addition to the indexes, we compare the simulated peaks of the chosen diagnostic and prognostic scenarios. Because the gauged watersheds areas range between 40 and 18,000 km 2 , we performed a scale-independent comparison. To obtain scale-independent peaks ( ), we divided the peaks [m 3 ⋅ s −1 ] by the mean annual peak ̅̅̅̅ [m 3 ⋅ s −1 ]. Then, for each event of each link, we computed the difference between the simulated ( ) and the observed ( ) standardized peaks. The peak flow estimation of D4 and P1 exhibited a similar performance, with D4 being superior. The D4 scenario reaches a fraction of In addition to the indexes, we compare the simulated peaks of the chosen diagnostic and prognostic scenarios. Because the gauged watersheds areas range between 40 and 18,000 km 2 , we performed a scale-independent comparison. To obtain scale-independent peaks (Z), we divided the peaks Q p [m 3 · s −1 ] by the mean annual peak Q p [m 3 · s −1 ]. Then, for each event of each link, we computed the difference between the simulated (Z s ) and the observed (Z o ) standardized peaks. The peak flow estimation of D4 and P1 exhibited a similar performance, with D4 being superior. The D4 scenario reaches a fraction of 32% for differences near zero (Figure 12), while in P1 this value drops to 28% ( Figure 12). Additionally, P1 has a higher fraction of errors equal or greater than one. 32% for differences near zero (Figure 12), while in P1 this value drops to 28% ( Figure 12). Additionally, P1 has a higher fraction of errors equal or greater than one. We expected the diagnostic superiority because, in the prognostic case, we impose restrictions based on maps. The resulting differences among simulations emphasize the parameters' relevance and the need for their correct representation. Contrasted with the diagnostic scenarios, the prognostic scenarios tend to reduce the performance. The differences between D4 and P1 suggest that the landscape descriptors could have errors that lead to decreases in the modeling performance. Additionally, our results suggest that there may be more tiles than the ones represented by the map in Figure 4b.

Analysis of Parameter Values
The diagnostic and prognostic scenarios offer different ways to determine the values of parameters in space. In the diagnostic cases, we identified the best fixed-parameter combination for each gauged watershed. In the prognostic cases, we pre-defined a set of distributed parameters based on the available information. In a previous step, we defined the best diagnostic and prognostic setup for each gauge (Figure 8a,b, respectively). According to our results, a spatial distribution of the parameters seems to be explained by We expected the diagnostic superiority because, in the prognostic case, we impose restrictions based on maps. The resulting differences among simulations emphasize the parameters' relevance and the need for their correct representation. Contrasted with the diagnostic scenarios, the prognostic scenarios tend to reduce the performance. The differences between D4 and P1 suggest that the landscape descriptors could have errors that lead to decreases in the modeling performance. Additionally, our results suggest that there may be more tiles than the ones represented by the map in Figure 4b.

Analysis of Parameter Values
The diagnostic and prognostic scenarios offer different ways to determine the values of parameters in space. In the diagnostic cases, we identified the best fixed-parameter combination for each gauged watershed. In the prognostic cases, we pre-defined a set of distributed parameters based on the available information. In a previous step, we defined the best diagnostic and prognostic setup for each gauge (Figure 8a,b, respectively). According to our results, a spatial distribution of the parameters seems to be explained by k i and the parameters of Equation (8). Additionally, the chosen diagnostic scenarios outperform the chosen prognostic ones (Figure 8a,b). In some gauges, the performance differences are small; however, in others, the difference is relatively large. This is an interesting result because the only difference between both cases is the parameterization of Equation (8). Considering the described performance differences, we explore in more detail how they are related to the parameterizations of the diagnostic and prognostic setups.
We compare Equation (8) setup for the diagnostic and prognostic scenarios to evaluate whether variations in the parameters explain the observed performance differences. We made the comparison at each gauged watershed. For the comparison, the prognostic setup has a set of curves q sL (P) for a given watershed (light blue lines in Figure 13), and there is one diagnostic curve q sL (D) for the same watershed (dark blue line in Figure 13). Using Equation (11), we compare q sL (D) with the 50th percentile of q sL (P) for storages between 1.6 and 1.7 [m] (green region in Figure 13). We choose the 1.6-1.7 range because it describes the shift the linear to the exponential expression. Besides, the hillslopes subsurface storage (S s ) is usually in this range during storm events.
According to Figure 14, the differences of the parameters (Δ ) do not adequately explain the performance differences between the diagnostic and the prognostic scenarios. In some cases, low absolute values of Δ are linked to low KGE differences. However, the described behavior does not apply for large absolute values of Δ . Figure 14 shows many watersheds in which the differences of the absolute parameters are larger than 10% (x-axis), while the KGE absolute differences are low. Additionally, some cases with low absolute Δ exhibit large KGE differences. On the other hand, according to the colors of Figure 14 (non-absolute Δ ), positive values of Δ are related to low KGE differences; and negative values of Δ correspond to high KGE differences. Figure 13. Example of the q sL parameterization comparison between the diagnostic and prognostic scenarios. The light blue lines are q sL curves of the prognostic scenario for a given watershed. The red line corresponds to the 50th percentile of the prognostic q sL curves. The dark blue line is the q sL curve of the diagnostic scenario. We used the green region to perform a comparison between the scenarios.
According to Figure 14, the differences of the parameters (∆q sL ) do not adequately explain the performance differences between the diagnostic and the prognostic scenarios. In some cases, low absolute values of ∆q sL are linked to low KGE differences. However, the described behavior does not apply for large absolute values of ∆q sL . Figure 14 shows many watersheds in which the differences of the absolute parameters are larger than 10% (x-axis), while the KGE absolute differences are low. Additionally, some cases with low absolute ∆q sL exhibit large KGE differences. On the other hand, according to the colors of Figure 14 (non-absolute ∆q sL ), positive values of ∆q sL are related to low KGE differences; and negative values of ∆q sL correspond to high KGE differences.
In some cases, low absolute values of Δ are linked to low KGE differences. However, the described behavior does not apply for large absolute values of Δ . Figure 14 shows many watersheds in which the differences of the absolute parameters are larger than 10% (x-axis), while the KGE absolute differences are low. Additionally, some cases with low absolute Δ exhibit large KGE differences. On the other hand, according to the colors of Figure 14 (non-absolute Δ ), positive values of Δ are related to low KGE differences; and negative values of Δ correspond to high KGE differences.

Figure 14.
Parameters' absolute differences vs. KGE differences at the USGS gauges. The colors correspond to non-absolute differences in the parameters. We also compared the ∆q sL and the KGE differences in space. According to Figure 15, the coincidences between the KGE and ∆q sL do not show a strong regional pattern. We observe some similarities only in the Des Moines Lobe and the Northwest Iowa Plains regions. In both cases, some significant KGE differences correspond with large absolute ∆q sL values. Additionally, there is a match between low KGE and ∆q sL differences in the Iowan Surface region, with some exceptions. We also compared the Δ and the KGE differences in space. According to Figure  15, the coincidences between the KGE and Δ do not show a strong regional pattern. We observe some similarities only in the Des Moines Lobe and the Northwest Iowa Plains regions. In both cases, some significant KGE differences correspond with large absolute Δ values. Additionally, there is a match between low KGE and Δ differences in the Iowan Surface region, with some exceptions. It is hard to establish a relationship between the diagnostic and prognostic parameters and their performance differences. We attribute this lack of correlation to the model's non-linear transformations at the hillslope scale and throughout the network. It is expected that the parameters would alter the model's output. However, our results show that a pre-defined distribution of the parameters could lead to modeling errors that are hard to identify.

Conclusions
The Iowa Flood Center (IFC) has been making operational flood forecasts for the state of Iowa since 2010. IFC forecasters use the hydrological Hillslope Link Model (HLM), along with rainfall data. The HLM has been accurate in forecasting streamflow fluctuations at several scales [14]. However, the model has limitations in its representation of the recession curve, and it underestimates the total streamflow volume. Moreover, the limitations seem to increase in a tiled landscape. The authors of [17] attributed these limitations to the linear equation HLM uses to represent the subsurface flux and the lack of an equation representing the tiled terrain. To address this issue, [15] developed an exponential equation that can be parameterized to represent the function of the hillslope steepness and the presence of tiles.
This paper evaluated the exponential equation proposed by [17], which represents subsurface hillslope-link interaction in HLM. The equation can represent hillslopes with It is hard to establish a relationship between the diagnostic and prognostic parameters and their performance differences. We attribute this lack of correlation to the model's non-linear transformations at the hillslope scale and throughout the network. It is expected that the parameters would alter the model's output. However, our results show that a pre-defined distribution of the parameters could lead to modeling errors that are hard to identify.

Conclusions
The Iowa Flood Center (IFC) has been making operational flood forecasts for the state of Iowa since 2010. IFC forecasters use the hydrological Hillslope Link Model (HLM), along with rainfall data. The HLM has been accurate in forecasting streamflow fluctuations at several scales [14]. However, the model has limitations in its representation of the recession curve, and it underestimates the total streamflow volume. Moreover, the limitations seem to increase in a tiled landscape. The authors of [17] attributed these limitations to the linear equation HLM uses to represent the subsurface flux and the lack of an equation representing the tiled terrain. To address this issue, [15] developed an exponential equation that can be parameterized to represent the function of the hillslope steepness and the presence of tiles. This paper evaluated the exponential equation proposed by [17], which represents subsurface hillslope-link interaction in HLM. The equation can represent hillslopes with and without tile drainage. We performed the equation evaluation at 140 USGS gauges in Iowa. The analysis used hourly records between 2002 and 2018. In the evaluation, we compared the exponential equation with a linear equation. The comparison used a diagnostic and a prognostic approach to establish the parameters. In the diagnostic setup, we implemented 12 fixed parameter scenarios, while in the prognostic setup, we distributed the parameters with consideration of the hillslope steepness and presence of tiles. In both cases, we considered three fixed percolation rates. Results from this study indicate the following:

•
Compared with the linear equation, the exponential equation corrects the volume bias on the simulated streamflow. We attribute the correction to the active layer threshold on the exponential equation and the significant outflow increase once the storage is above this threshold. In contrast, in the linear equation, the water remains in the soil for extended periods because of the described absence of these processes.

•
Depending on the parameters, the exponential equation could improve the performance of HLM. We found that the exponential equation outperforms the linear equation for several parameter combinations with changes in the shape of the hydrograph, the simulated peaks, and timing. We also found significant differences using different combinations of the equation parameters and the percolation rate.

•
The percolation rate plays a significant role in the representation of the subsurface flux from the described combinations. We found spatial coincidences in the percolation rates when choosing the best diagnostic and prognostic scenarios. Additionally, the percolation rate induces changes comparable to those produced by the exponential equation's parameters.

•
Determining the distributed parameters of HLM remains challenging. In this paper, we used the diagnostic and prognostic approaches to analyze the parameters of HLM. The diagnostic approach assumes unknown conditions and fixed parameters over the space. On the other hand, the prognostic method is the more classical approach, in which the parameters are derived from maps of the landscape. In our experiments, the diagnostic setups tended to outperform the prognostic setups. Additionally, had difficulty in identifying a link between the diagnostic and prognostic parameters and their respective performances.
In the current work, we showed that a better representation of the processes and the correct parameters can improve a hydrological model. The improvement is supported by comparisons performed at 140 USGS gauges. Moreover, the differences between the diagnostic and prognostic setups suggest that identifying the parameters is still challenging. Despite the limitation related to the number of gauges, the diagnostic approach reveals the parameters' potential spatial distribution.
Two main factors may explain the differences in parameters and performance between the diagnostic and prognostic setups: errors in the landscape description and unrepresented processes in HLM. Uncertainties exist in the tile localization maps; likewise, limitations exist in the representation of the average steepness at the hillslope scale. On the other hand, we unrepresented processes in some regions of Iowa, such as potholes over the northwest and agricultural terraces in the west. It is difficult to identify which of these factors is more relevant to the implementation of a hydrological model. However, according to our results, the use of maps as landscape descriptors may lead to the detection of errors that are usually hidden in a posterior calibration process. Moreover, we found it difficult to identify the errors caused by the prescribed distributed parameters. Both issues could be addressed using diagnostic setups that help to identify the uncertainties derived from the parameters and their possible regional distributions. Data Availability Statement: The data supporting the results of this work can be found in the following repository.