Improvement of surface run‐off in the hydrological model ParFlow by a scale‐consistent river parameterization

We propose an improvement of the overland‐flow parameterization in a distributed hydrological model, which uses a constant horizontal grid resolution and employs the kinematic wave approximation for both hillslope and river channel flow. The standard parameterization lacks any channel flow characteristics for rivers, which results in reduced river flow velocities for streams narrower than the horizontal grid resolution. Moreover, the surface areas, through which these wider model rivers may exchange water with the subsurface, are larger than the real river channels potentially leading to unrealistic vertical flows. We propose an approximation of the subscale channel flow by scaling Manning's roughness in the kinematic wave formulation via a relationship between river width and grid cell size, following a simplified version of the Barré de Saint‐Venant equations (Manning–Strickler equations). The too large exchange areas between model rivers and the subsurface are compensated by a grid resolution‐dependent scaling of the infiltration/exfiltration rate across river beds. We test both scaling approaches in the integrated hydrological model ParFlow. An empirical relation is used for estimating the true river width from the mean annual discharge. Our simulations show that the scaling of the roughness coefficient and the hydraulic conductivity effectively corrects overland flow velocities calculated on the coarse grid leading to a better representation of flood waves in the river channels.

and the subsurface described by precipitation, evapotranspiration, infiltration into the subsurface, exfiltration to the surface, and overland flow. Physically-based distributed hydrological models simulate the evolution in time and space of the state variables and fluxes, and can be used to calculate and predict discharge and/or water levels at each grid-column and time-step and to compensate for the lack of sufficient observations. Observations are, however, required for model calibration, initialization, and validation.
In this paper, we focus on overland flow as part of the hydrometeorological system. Several methods exist for simulating overland flow such as the Saint-Venant equations or approximations thereof like kinematic or diffusive waves. Even simpler routing methods exist based only on source and sink terms (Clark et al., 2017(Clark et al., , 2015de Rooij, 2017). Some models calculate overland flow at the grid scale, whereas others use routing schemes for separate channels not related to the grid resolution. We target the prior model types for which grid scales are usually much larger than the size of the real rivers. Realistic overland flow simulations despite low spatial resolutions would limit the usually high computational demand of such models (Clark et al., 2015). Much higher resolutions-and thus substantially increased IT resources-are required when the scale dependence of lateral fluxes are explicitly taken into account (Wood et al., 2011).
The distributed model ISBA/MODCOU (Masson et al., 2013, Ledoux, Girard, Marsily, & Deschenes, 1989, Häfliger et al., 2015 is one example of the model type discussed; the model is used to simulate discharge and river water levels over France. Other examples are the global-scale Catchment-based Macroscale Floodplain (CaMa-Flood; Yamazaki, de Almeida, & Bates, 2013) and the ISBA-TRIP  models, which are run at relatively low resolution but use subscale digital elevation model (DEM) information to characterize floodplain areas within individual grid cells. The limited resolution of such regional or large scale models requires the parameterization of subgrid scale processes especially for overland flow and river routing. Neal, Schumann, and Bates (2012) studied the effect of subgrid scale channel routing on flood dynamics by using a hydraulic model for the Niger River in Mali, which allows for river channels of any width; they showed that more detailed (smaller) channels considerably improved model performance.
An alternative to explicit sub-grid channel routing is the scaling of relevant parameters in grid-scale river routing models. Schalge et al. (2017) show that without any scaling, the grid-scale routing scheme at a resolution of 400 m delays the discharge peak at the mouth of the Neckar catchment (south-western Germany) by 3 days, whereas peaks are underestimated compared with observations. Sub-scale parameterizations have been developed based on the subgrid scale topographic index, for example, by Niedda (2004), which has been refined, for example, by Fang, Bogena, Kollet, and Vereecken (2016), and shown to improve river discharge. Other methods employ probabilistic approaches (Piccolroaz et al., 2016) to account for uncertainty related to overland flow.
In this study, we propose an improvement of the overland-flow parameterization for distributed hydrological models with constant horizontal grid resolution, which employ the kinematic wave approximation for both hillslope and river channel flow. Without subscale parameterization, such models simulate reduced river flow velocities for streams narrower than the horizontal grid resolution. Moreover, surface areas, through which these wider model rivers may exchange water with the subsurface, are larger than the real river channels potentially leading to unrealistic vertical flows. We propose an approximation of the subscale channel flow by scaling Manning's roughness in the kinematic wave formulation via a relationship between river width and grid cell size following a simplified version of the Barré de Saint-Venant equations (Manning-Strickler equations, see Section 2.2 for details). The too-large exchange areas between model rivers and the subsurface are compensated by a grid resolution-dependent scaling of the infiltration/exfiltration rate across river beds for the top layer. Both the roughness coefficient and the hydraulic conductivity can be objectively scaled by assuming a rectangular river channel cross section with variable width embedded in the fixed-size river grid cells. The approach is simple to implement and requires no additional computational demand.
We implemented the approach in the distributed hydrological model ParFlow (Kollet & Maxwell, 2006), which simulates surface run-off and subsurface flow processes, and was also used in Schalge et al. (2017). ParFlow uses the same equally-spaced horizontal grid for 3D subsurface flow and 2D overland flow. For computational reasons, the horizontal resolution is on the order of several hundred metres if run over medium (~50,000 km 2 ) or large-scale (>150,000 km 2 ) catchment sizes. To our knowledge, such an approach has not been tested before to improve river behaviour, though a similar approach has been recently used to infer a scaling law for hydraulic conductivity for entire catchments by Foster and Maxwell (2019).
The ParFlow model is introduced briefly in Section 2, with a focus on shallow overland flow and variably saturated groundwater flow and on its integration in the fully coupled land surfaceatmosphere modelling framework TerrSysMP. Our scaling approaches for Manning's coefficient and hydraulic conductivity are described in Section 3. After presenting the experimental design in Section 4, we demonstrate in Section 5 the usefulness of the approach in synthetic numerical experiments and in an application to the real-world Neckar catchment in southern Germany. Section 6 summarizes the results and concludes with implications and suggestions for future studies.

| THE PARFLOW MODEL
ParFlow is a 3D variably saturated groundwater flow model with a free-surface overland flow boundary condition (Jones & Woodward, 2001;Kollet & Maxwell, 2006). In a continuum approach, ParFlow is used to simulate the coupled surface-subsurface hydrodynamics, including the redistribution of soil moisture and the flow of groundwater and surface water in an integrated fashion.

| Variably saturated groundwater flow
The mass balance of variably saturated flow following Richards (1931) relates the temporal variation of the degree of saturation S w [−] in soil to the water (Darcy) flux divergence q [L T −1 ]. Assuming no sinks or sources in the domain this yields, is the relative permeability depending on pressure head ψ p , and z [L] is the gravitational potential with the positive z-axis pointing upward.
The hydraulic conductivity K sat depends on soil texture, which is parameterized in ParFlow by the van Genuchten relationship (Van Genuchten, 1980): (2b)

| Shallow overland flow
The solution of the Barré de Saint-Venant equations describes the routing of water over a rough surface. In general, the resulting flood wave has advective and diffusive properties. The kinematic wave approximation, which is used in ParFlow, neglects the diffusive part.
Under this assumption, Manning's equation can be used to establish a flow depth-discharge relationship for a rectangular river channel where Q is the river discharge [L 3 T −1 ], v is the average flow velocity In ParFlow, overland flow is incorporated via an overland flow boundary condition assuming pressure and flux continuity at the surface-subsurface interface (Kollet & Maxwell, 2006). Thus, the system of equations of variably saturated groundwater and overland flow is coupled via the boundary condition at the ground surface leading to with ||A,B|| the greater of A and B. Pressure head continuity (ψ s = ψ p = ψ) states that the pressures of the surface and subsurface domains vary continuously across the land surface.

| RIVER PARAMETERIZATION
ParFlow does not differentiate between hillslope run-off and river flow, for example, by accounting for specific channel geometries, and the same horizontal grid resolution applies both for the subsurface and surface water domains. Because especially the latter depends on unresolved subgrid scale surface topography, negative effects on water flow may result. This could be compensated either by incorporating subgrid scale river channel geometries, including the relevant exchange fluxes with the subsurface, for example, via conductance concepts or by scaling grid-scale parameters according to the used resolution. In order to utilize the overland flow boundary condition in ParFlow, we explore the latter option and derive grid resolutiondependent scalings of the Manning's coefficient n and the saturated hydraulic conductivity K sat .

| Scaling the Manning's coefficient
The proposed scaling of the Manning's coefficient n adjusts the flow velocity v in a grid cell with a river channel to the flow velocity in a "true" river channel of smaller width W. n can be objectively scaled by assuming a rectangular river channel cross section with a certain width contained in the larger river grid cell ( Figure 1). Let the true situation (given, e.g., by a high resolution simulation) be a river channel of width W 1 with a Manning's coefficient n org and a flow velocity v 1 . In a coarser-scale model, the river has a (wider) width W 2 equal to the width of the grid cell. The original Manning's coefficient n org will then result in a lower flow velocity v 2 (compared with v 1 ) because of the smaller water depth in the wider channel. Whereas keeping the model river width at W 2 (the size of a coarse grid cell), we scale (reduce) the original Manning's coefficient n org to n scale , such that the resulting flow velocity v 3 is equal to the true velocity v 1 . In order to obtain n scale , we impose for all three cases the same discharge Q in the upstream part of the respective virtual channel; thus from Equation (3) follows where ψ 1 is the true ponding depth of the river. With the same Manning's coefficient n org and discharge Q. For a larger (model) river width W 2 follows in the same manner a smaller ponding depth ψ 2 of the model river. Because W 2 > W 1 and ψ 2 < ψ 1 , it follows v 2 < v 1 .
We require from the scaled Manning's coefficient n scale , that the original (true) flow velocity is conserved (v 3 = v 1 ). Using n scale in (5a) and (5b) yields which can be solved for the scaled Manning's coefficient n scale FIGURE 1 Illustration of the scaling concept for the Manning's coefficient (see text). W 1 is the true channel width and W 2 the (wider) width of the grid cell and thus the model river width. The flow velocities v i are indicated by the length of the grey arrows and the vertical lines represent the position of peak discharge at each time step tn In Equation (7) is the scaling coefficient for Manning's coefficient, which results in the original (true) river flow velocity for a model channel width W 2 independent of channel slope S f and discharge Q.

| Scaling the hydraulic conductivity
Because the model river width is typically larger than the true river width, a larger surface area will exchange water with the subsurface compared with the real river. This can be corrected by using a scaled (reduced) hydraulic conductivity K sat scale (scaled K sat , case 4 in Figure 2). Conserving the original (true) infiltration/exfiltration fluxes Q in/ex (Q in/ex case 4 = Q in/ex case 1 ) requires the following: with A [L 2 ] as the horizontal area of the river in the grid cell Assuming a homogeneously distributed hydraulic conductivity K sat over the entire grid cell, we get and finally with ϰ ¼ W 1 W 2 as the scaling coefficient for the hydraulic conductivity.

| EXPERIMENTAL DESIGN
In order to verify the proposed scaling of n and K sat , an idealized and a real-world test case were implemented in ParFlow. The idealized test case considers a simple domain where the truth (high resolution) is known. Thus, we can directly validate both the approach and its implementation in ParFlow. The real-world test case is set up for the Rhine-Neckar area in order to evaluate the approach by comparison with real observations. Here, we are faced with realistic heterogeneity, which leads to additional nonlinear feedbacks of the scaled parameterizations with the model. Although a direct comparison with observations is difficult due to additional uncertainties in the real land surface processes and meteorological forcings, we can use the results discussed in Schalge et al. (2017) to evaluate whether our scaling approach generates the expected changes in the catchment response.

| Idealized test case
The idealized test case (      the full domain with n = 5.52 × 10 −5 hr m −1/3 . The hydraulic conductivity K sat ranges from 0.001 to 0.1 m hr −1 in the surface layer.

Initialization and experimental set-up
The initial condition for pressure head ψ is obtained from a 5-year

| Scaling for the Neckar catchment area
In order to scale the hydraulic parameters via Equations (7) and (11) for each river grid cell, the true river widths (W 1 ) are needed, which we set to an effective river width W 1 , estimated from known upstream discharge values via an empirical relationship by Leopold and Maddock (1953) as a function of the mean annual discharge Q, with k = 7.12 and m = 0.53. The parameters k and m were fitted, tested, and validated by Häfliger et al. (2015) for the Garonne catchment (55,846 km 2 ), whose catchment area is close to the simulated domain (57,850 km 2 ), covering all the Rhine-Neckar area. Figure 5 shows the resulting effective river width (spatial and cumulative distribution) for the Rhine-Neckar area, which varies between 10 and 150 m, with a mean width of about 34.5 m and a standard deviation The scaled Manning's coefficient n scale (Figure 6) ranges from 0.5 × 10 −5 hr m −1/3 in mountain rivers to 2.7 × 10 −5 hr m −1/3 in the lowland. The scaling parameter ϰ for the hydraulic conductivity K sat org ranges from 0.025 in the mountains to 0.33 in the lowlands. For example, the hydraulic conductivity K sat org for a 10-m and a 150-m wide river would be scaled by 0.025 and 0.3, respectively, in order to result in in/exfiltration rate comparable with those of the true rivers. The resulting scaled coefficients K sat scale range from 0.5 mm hr −1 in mountain rivers to 6 mm hr −1 in low lands.

| Idealized testcase
In all figures showing results of the idealized test case setup for different spatial resolutions, the Nash-Sutcliffe Criterion (NSC) (Nash & Sutcliffe, 1970) was computed by taking the results of the high resolution (40 m) setup as observations. Figure 7 compares the hydrograph at the channel outlet simulated at the 40 m resolution for the three initial water table conditions for the low rain case. The initial condition has only a small impact on the discharge characteristics, whereas the effect of a higher infiltration can be observed for the low water table case. For the high rain case, the results for all setups become indistinguishable (not shown). The effect of the different scaling applied at 400 m resolution is shown in Figure 10. For a high water table and a low rain rate, the scaling of the hydraulic conductivity alone has little impact; only

| Semi-idealized Neckar test case
The results of the semi-idealized experiment with a single half-day precipitation event when scaling the Manning's coefficient alone  Figure 14). The dry period in the beginning of the simulation removes potential footprints of past rain events. Scaling the Manning's coefficient enhances the discharge peak after the rain event (+25%) and shifts it earlier by 2 days. The additional scaling of the hydraulic conductivity further increases the discharge peak by 5% to 15%. Overall, these results compare well with those obtained for the ideal case for a high rain rate.

| DISCUSSION AND CONCLUSION
We suggest an improvement of river hydrographs simulated with the kinematic wave approach in distributed low hydrological models by

FIGURE 11
Hydrographs for the case with low rain rate (0.005 mm/ hr) and WT = 0 m initial conditions for the 40 m case without scaling (ex1) and the 200 m cases without scaling (ex31), with Manning's scaling (ex37), with K sat scaling (ex43) and both Manning's and K sat scaling (ex49)

FIGURE 12
Hydrographs for the case with high rain rate (0.01 mm/ hr) and WT = 0 m initial conditions for the 40 m case without scaling (ex4) and the 400 m cases without scaling (ex58), with Manning's scaling (ex64), with K sat scaling (ex70) and both Manning's and K sat scaling (ex76)

FIGURE 13
Hydrographs for the case with low rain rate (0.005 mm/ hr) and WT = −5 m initial conditions for the 40 m case without scaling (ex3) and the 400 m cases without scaling (ex57) Although we see an improvement in the response to the rain event, we have to acknowledge that the Neckar catchment is heavily managed, effects of which are not included in the simulations. Thus, we suggest to apply our approach in catchments with minimal human influence (see e.g., Caballero et al., 2007or Pedinotti et al., 2012. The choice of the threshold value for the computation of effective river widths W 1 in real catchments influences the number of river grid cells and the length of the river reaches in which scaling operations are performed. A higher threshold will reduce the scaled river cells and the length of the scaled river reaches, and vice versa. Thus, the choice of the threshold will also influence the propagation of the simulated flood waves. In ParFlow, smaller rivers may disappear completely during dry periods. Scaling of even these smaller rivers can lead to an improvement by keeping the river network as a whole; but the base flow of the larger rivers will only be marginally impacted. Although there is no objective lower limit, we recommend to refrain from river scaling, if either the river is not present over more than 50% of the simulated time in the unscaled experiment or if the river is that small that nonlinearities introduced by the riverbed morphology become important (width at and below 2 m). All this depends on the model resolution; as we could see in the idealized study, scaling becomes less successful for real channel width less than 1/10th of the grid resolution. Thus, the approach is likely not suitable for regions with pronounced seasonal variations in rainfall (monsoons or similar) which cause rivers to behave differently on longer timescales than typical floods in the mid-latitudes. Similarly, rivers that only appear due to heavy rain events cannot be scaled properly as there is not enough data to infer what the typical river width would be in that case.
Through the scaling of the saturated conductivity also soil moisture, base flow, and the infiltration/exfiltration between the river and the aquifer will be impacted. Following Darcy's equation (Equation (1b)), the vertical infiltration or exfiltration between river and aquifer depends on the hydraulic head gradient between the top soil layer of the river grid cell and the layer directly below. By scaling Manning's coefficient, the ponding depth is lower, which reduces the infiltration/exfiltration flux between river and aquifer and might even reverse its direction. These processes should be analysed in more detail for a real river network in search for an improved K sat scaling, because our current approach neglects nonlinear interactions between river and subsurface. Our method might also lead to only little benefits, when the riverbed morphology-especially the width-discharge relation-largely differs from the empirical values we used, or when features like non-resolved waterfalls exist.
One last point of concern related to the way the rivers are simulated. Because of the low spatial resolution used in ParFlow, typical shapes of river valleys are not well represented, which affects also the stream-aquifer interaction since the model groundwater dynamics depends on the grid resolution (Refsgaard, 1997).
The developed methodology can be easily implemented in all models, which do not explicitly resolve the true river width for river routing when the width is known with reasonable accuracy. Practically only a preprocessing step is required, which does not increase computational demand during runtime. To our knowledge, no similar methods have been tried so far as most approaches rely on dedicated channel parameterizations that are much more complex to implement.