Flow and Transport Modeling in Heterogeneous Sediments Using an Integral Approach

An integral approach which can simultaneously model turbulent flow and transport at the sediment–water interface has been recently developed and validated for homogeneous sediment which was achieved by comparing numerical results to flume experiments on flow and transport over a rippled streambed and through the sediment for neutral, gaining, and losing conditions. In the present study, we validated the approach for heterogeneous conditions by comparing numerical simulations of flow and transport in heterogeneous sediment to analytical solutions as well as flume experiments on flow and transport through rippled streambed consisting of heterogeneous sediment. For this complex setup, simulation and experimental results agree well showing that flow and tracer transport prefer paths through areas with bigger grain diameters and higher porosities. The effect of flow redirections under losing and gaining conditions on hyporheic flow and residence times is discussed.


Introduction
In the hyporheic zone, which is the transition zone between aquifer and river, groundwater, and surface water are closely coupled (Cardenas 2009;Boano et al. 2014;Pryshlak et al. 2015;Lewandowski et al. 2019). This zone includes riverbeds, riverbanks, saturated sediments under dry bars and riparian areas (Edwards 1998). A multidirectional, multiscale flow exchange between groundwater and surface water occurs and is controlled by factors such as stream morphology, ambient groundwater, and sediment heterogeneity (Buffington and Tonina 2009;Cardenas 2009;Pryshlak et al. 2015). Hyporheic exchange transports contaminants (McKnight et al. 2001;Medina et al. 2002;Feris et al. 2003), organic carbon, and nutrients (Valett et al. 1996;Hill and Cardaci 2004) from the surface water to the streambed. The streambed is populated by microbial biomass (Stonedahl et al. 2012). Biomass, hyporheic exchange fluxes, and residence times in sediments are crucial factors in determining geochemistry of stream ecosystems (Arnon et al. 2013;Trauth et al. 2015;Fox et al. 2016).
Transport, settling, and remobilization of sediment particles under different flow conditions result in heterogeneous sediment composition (Powell 1998). Compared with a homogeneous sediment, heterogeneity modifies flow paths, exchange fluxes, solute transport, and residence times (Cardenas et al. 2004;Salehin et al. 2004;Singha et al. 2008;Sawyer et al. 2011 andCardenas 2009;Stonedahl et al. 2018). The influence of morphology on hyporheic exchanges has been addressed through many physical and numerical setups, by using different bedforms and structures (e.g., Elliott and Brooks 1997;Stonedahl et al. 2010;Broecker et al. 2019). Due to its major impact on hyporheic exchange processes in streams, heterogeneity should be given high priority in presence of morphology (Freeze and Cherry 1979;Pryshlak et al. 2015). In studying hyporheic flow through morphological structures such as bed forms, heterogeneity adjusts flow by redirecting it through areas with higher hydraulic conductivities (Vaux 1968). Other major factors which affect hyporheic exchange flows are groundwater down-and upwelling. The effects of losing and gaining groundwater conditions have been physically (Fox et al. 2014;Gomez-Velez et al. 2014;Fox et al. 2016) and numerically modeled . Figure 1 shows a cross-section through the hyporheic zone with potential subsurface flow paths. Representation of hyporheic exchange flows is best achieved by considering stream morphology, sediment heterogeneity, and ambient groundwater simultaneously. Field investigation addressing these aspects simultaneously requires setups which are extremely challenging (Fox et al. 2016). Instead, a combination of flume experiments which consider the above-mentioned factors with a modeling tool which can be calibrated using the experiment and can be used to track flow and transport of solutes offers valuable insight in hyporheic exchange flows.
To visualize flow and solute exchange in hyporheic zones, Fox et al. (2014) investigated the propagation of a tracer into a rippled streambed with a unidirectional surface water flow over bedforms. The flume consisted of a homogeneous streambed (including bedform morphology) with overlying flowing water. At the flume bottom gaining or losing conditions were simulated. Broecker et al. (2021) have modeled the flume experiment of Fox et al. (2014) with a newly developed integral approach which has advantages compared with coupled approaches (e.g., Saenger et al. 2005;Chen et al. 2018). Instead of using the interface between stream and groundwater as the transitional zone where results of surface water modeling are applied as pressure boundary condition to a groundwater model (as is common with one-way coupling), integral approach allows the continuous exchange between two zones by using a modified form of Navier-Stokes equations which can include the porous zone directly in the flow model. Sobhi Gollo et al. ( , 2022 have shown that compared with prevalent coupled approaches, the integral approach allows capturing flow and pressure distributions at the interface between groundwater and surface water more realistically. Transport across a homogenous rippled streambed has been investigated by flume experiments of Fox et al. (2014) and simulations of Broecker et al. (2021). Fox et al. (2016) used a setup very similar to the one of Fox et al. (2014) and included heterogeneity. They generated a random field with three sediment types and tracked flow and transport through the heterogeneous sediment.
The aim of the present study is to validate the integral approach of Broecker et al. (2021) for tracking flow and transport in a heterogeneous hyporheic zone using the experiments of Fox et al. (2016). To do so, we first compare the model results to a one-dimensional analytical solution to check if the transport of a tracer in different porous zones can be captured correctly. We then reproduce the exact sediment system presented in Fox et al. (2016). Next, we investigate flow through the heterogeneous sediment and introduce tracer accordingly for flow and transport validation.

Computational Domain
The numerical model geometry is based on the circulating flume experiments of Fox et al. (2016). Via obtaining reasonable anisotropy values and meeting ergodicity requirements, conditions close to natural streambed are imitated in the geometry of Fox et al. (2016) experiment which makes it a representative bed texture to study the effect of streambed heterogeneity on the hyporheic exchange fluxes. To reduce the computational effort, the numerical model is smaller (2.375 m) in length compared with the original experimental setup (6.4 m). The model width (in y-direction, see Figure 2) is 0.29 m for both numerical and experimental setups. The inflow from the inlet enters a ramp before flowing over the rippled  streambed (as done in Fox et al. 2016). The ramp slope starts at the base of the inlet and ends at the height of 0.2 m and is 0.93 m long in x -direction ( Figure 2). The height of the water from ripple crests is set to 7 cm. Identical to Fox et al. (2016), the sediment is packed with 20 layers of heterogeneous sediment with 1 cm spatial resolution in z -direction. The ripples (dune-shaped forms) are each 20 cm long and 2 cm high. The ripple crest is positioned 15 cm from the left and 5 cm from the right foot of the ripple.
A mesh was generated using GMSH (Geuzaine and Remacle 2009) using different element types. To enable mapping the heterogeneity in the sediment, the sediment part (ripples excluded) is meshed with cuboid elements. Surface water and ripples are meshed with unstructured elements. The model consists of 72,000 elements in total. To account for the steep gradients at the surface water-groundwater interface, mesh elements are finer on the interface with minimum cell size of 1 × 10 −3 m 2 compared with maximum of 8 × 10 −3 m 2 in some distance above the interface. In x -z plane, each 1 cm × 1 cm × sediment cell from Fox et al. (2016) is represented in the numerical grid by four 0.5 cm × 0.5 cm × cells.

Boundary and Initial Conditions
The system is initially saturated with water and remains so during the modeling process for neutral, gaining, and losing conditions. The heterogeneous structure is generated using Fox et al. (2016) who used three different sediments (fine, intermediate, and coarse sand) with porosity values of 0.33, 0.44, and 0.44 and mean grain diameters of 0.38, 1.3, and 2.3 mm, respectively. Figure 3 shows the generated patterns according to the sediment pattern plan acquired through personal contact with the authors of Fox et al. (2016).
An overview of the boundary conditions is illustrated in Figure 4. Fox et al. (2016) defined a constant velocity of 0.15 m/s over the bedforms for a mean water level of 0.07 m and a flume width of 0.29 m which, for this setup, accounts for a fixed inlet volumetric flow rate of 0.003045 m 3 /s and inlet velocity of 0.0362 m/s. Slip conditions are defined for the top boundary. Ramp surface as well as left and right side of the sediment are defined with no slip wall conditions. The outflow at the flow outlet boundary equals the inflow at the flow inlet for neutral conditions. Under losing and gaining conditions, outflow equals inlet inflow plus/minus gaining/losing flow. The sediment bottom is defined as a no-flow boundary for neutral conditions. For gaining (q G ) and losing (q L ) conditions a 100 cm/day flow is defined in analogy to Fox et al. (2016). This value corresponds to a fixed inlet/outlet volumetric flow rate of 4.85 × 10 −6 m 3 /s in ±z -direction. Pressure is fixed to 0 Pa at the outlet. To adjust the pressure gradient according to boundary fluxes for losing and gaining conditions, a special boundary condition (called fixedFluxPressure in OpenFOAM) is applied to sediment bottom.
In all cases pre-runs of 5 min with initial surface water flow of 0.15 m/s in x -direction are carried out to reach quasi steady state. For losing and gaining cases, q L and q G are applied as well. Results of these preruns are then used as initial conditions and then an inert constant tracer with a concentration of 1 mg/L and a molecular diffusion coefficient of 10 −9 m 2 /s enters the system through the inlet boundary.

Governing Equations and Computational Tool
Simulations were carried out with OpenFOAM (Open-source Operation and Manipulation) software version 2.4.0. The porousInter solver, which is based on the InterFOAM solver, was developed by Oxtoby et al. (2013) and was used here to model the simultaneous flow in surface water and sediment. In this study, the sediment is fully saturated. Being based on an extended version of the Navier-Stokes equations, porousInter, considers porosity, and grain diameter of the sediment as additional resistance terms (see Equations 3-5). The equations for the conservation of mass and momentum are written as ([ ] f is averaged parameter over void region): where is the effective porosity, and α [kg/(m 2 s 2 )] is an additional porous drag term. α accounts for momentum loss by means of fluid friction with the porous medium and flow recirculation within the sediment. Porous drag term α is defined with Equation 3: where d p [m] is the effective grain diameter. Equation 4 defines pressure loss due to fluid friction with sediment after Ergun (1952). In sediment, compared with surface water flow, more momentum is needed to accelerate a given volume of water. This so-called added mass (van Gent 1995) provides extra momentum needed for a larger volume of water to be accelerated. Equation 5 considers this extra momentum due to flow recirculation by porous medium, after van Gent (1995). Where there is only surface water flow, effective porosity is set to 1 and therefore A = B = α = 0. This means that for areas without sediment, the standard Navier-Stokes equations are used. To apply the set of equations to the flow domain, OpenFOAM uses the semi-discrete Finite Volume Method (FVM) discretization scheme and for time it uses Finite Difference Method (FDM). The interface between water and air phase is captured via the Volume of Fluid (VoF) method. In OpenFOAM, velocity is discretized using sparse linear solver (smoothSolver, smoother: symGauss-Seidel). Pressure is calculated by using geometric agglomerated algebraic multigrid preconditioner (GAMG) solver.
Coupling of velocity and pressure is done via the PIM-PLE algorithm. PIMPLE (combined PISO and SIMPLE) is a semi-implicit method for pressure-velocity coupling that accounts for transient flow conditions.
For modeling transport processes, we use the extension of porousInter called porousInterTracer developed by Broecker et al. (2021). They have implemented the following advection-diffusion equation for transport using diffusion coefficient D [m 2 /s]: OpenFOAM, tracer concentration is discretized with BICCG (bioconjugate gradient) solver which is used for asymmetric matrices and has good parallel scaling.
For the present study, we initially used a simple and computationally efficient RANS turbulence model (k -ε; k = turbulent kinetic energy, ε = dissipation rate of turbulent kinetic energy). RANS was not capable of capturing flow velocity and pressure fluctuations in the model. Therefore, we then chose the LES turbulence model (Large Eddy Simulation [LES]; Smagorinsky 1963) which is more robust than the more commonly used RANS model for this study. Modeling the setup of Figure 2 by using 96 parallel processors of the North-German Supercomputing Alliance (HLRN) and applying LES turbulence model takes 24 days (15 days for k -ε) for a 60-min study of tracer propagation. Tracer expansion modeling with this turbulence model led to comparable results as in laboratory experiment serving as a justification for this approach.  Figure 5. An inert tracer with a diffusion coefficient of 10 −9 m 2 /s is continuously injected from the left side of the domain. Figure 5 shows the tracer progression under constant velocity of 0.01 m/s. Simulation results are compared with one-dimensional continuous injection analytical solutions of Kinzelbach (1992). A very good agreement between analytical solution results and simulations can be seen in Figure 5 for different timesteps.

Comparison of Numerical to Experimental Results for Tracer Transport
In this section, for the transport of a tracer into heterogeneous sediment the lab observations of Fox et al. (2016) are compared with the results of the integral model. Simulated tracer propagation is compared with experiment in Figure 6 under the sixth ripple. Numerical results of an alternative model (MIN3P, Mayer et al. 2002;modeled by Fox et al. 2016) after 40 min tracer propagation are also shown in Figure 6 and compared with the integral model in Table 2. Figure 6 shows that the tracer first spreads around the green zone in the upper left with its small hydraulic conductivity (smallest grain size diameter; Figure 6, neutral 10 and 40 min, losing 10 min, gaining 10, 40, and 60 min). Compared with the integral modeling, basically the same behavior is observed in experimental runs. The tracer propagates fastest under losing conditions and slowest under gaining conditions and prefers higher grain size diameter pathways under all conditions.
As suggested by Fox et al. (2016), to model sharp tracer fronts of the experiment, we assumed advection-diffusion dominant tracer propagation. By calculating the Peclet number, we have shown the contribution of advection and diffusion to transport in Figure 7 for tracer propagation at 60 min. We have displayed that the fastest tracer propagation occurs under losing conditions where advective transport is dominant. Tracer propagation is slower under neutral conditions and smallest under gaining conditions. Under these two conditions, diffusive transport is dominant (Figure 7). Figure 8 shows how minor differences in dyed areas between experimental runs and modeled simulations are probably caused by inaccurate implementations of the sediment geometry and pattern in the experiment. Although such inaccuracies can be considered negligible for an experimental setup, in comparison to modeled tracer fronts in heterogeneous sediments they can be important. In Figure 8 it is displayed how the original sediment setup and ripple geometry (as suggested by Fox et al. (2016) and applied for integral model) changed during sediment packing and experiment. Local changes of sediment pattern (e l ) compared with the original setup are illustrated for the investigated zone as an example. Throughout the experiment, these small changes are very important in opening and blocking pathways of flow and tracer propagation. Although compared with the original plan the generally heterogeneous sediment pattern is relatively intact (except e l ), a small bedform migration (e a1−5 ) and deformation (e b1−6 ) of the bedforms can be detected. Presumably, these changes of the sediment surface cause a formation of a thin layer right underneath the bedform. This unidentified (thickness unknown) sediment layer (e u ) is also displayed in Figure 8. As local and bedform geometry changes cannot be quantified exactly, original bedform geometry and sediment composition are used for modeling. The experimental sediment composition and bedform morphometry deviates only slightly from this plan as shown in Figure 8.
To compare results of experimental runs and simulations quantitatively, they are compared through photographs of tracer fronts at 10, 40, and 60 min. Applying "ImageJ" which is used for scientific image analysis, a comparison of the dyed areas of experiment and simulation is shown in Table 1   followed by gaining and neutral conditions. This indicates that less unexpected ripple geometry deformations throughout the experiment under losing conditions result in less discrepancy between simulation and experiment. Considering undefined experimental setup modifications, they fairly match experimental tracer propagation trends both qualitatively ( Figure 6) and quantitatively (Table 1). If experimental setup modifications could be definable in the model or could be eliminated from the experiment, a perfect match between simulation and experiment is to be expected.
Comparing the integral model results at 40 min to an alternative groundwater model (MIN3P; Figure 6; obtained from Fox et al. 2016) in Table 2 shows that the integral solver simulates the tracer propagation for all the conditions very similar to a widely used groundwater model like MIN3P.
By observing a good match between experimental and integral model results (as well as comparison with MIN3P), next, we will describe the modeled pressure distribution and velocity fields in the heterogenous sediment.

Hydrodynamic Pressure
Hydrodynamic pressure distributions under neutral, losing, and gaining conditions on fourth, fifth, and sixth ripples are shown in Figure 9. These ripples are sufficiently far away from the inlet and potential impacts of the flume inlet. Neutral conditions are defined by setting zero velocity to the lower boundary of the sediment bed. For losing conditions an outflux of 4.85 × 10 −6 m 3 /s is set for the bottom of the flume in analogy to Fox et al. (2016). For gaining conditions an influx of 4.85 × 10 −6 m 3 /s through the bottom of the flume is applied equal to the flux used in Fox et al. (2016).
Hydrodynamic pressure results from subtracting total pressure from hydrostatic pressure. In the present study, hydrodynamic pressures are higher on the luv (upstream)  side of the ripple as shown in Figure 9. Compared with losing conditions, gaining, and neutral conditions show lower hydrodynamic pressure. Hydrodynamic pressure values for a horizontal section underneath the ripples (z = 0.195 m) of Figure 9 are displayed in Figure 10. The 60 cm long horizontal cross section (1.575 m < x < 2.175 m, z = 0.195 m) which cuts through several sediment types with different grain diameters is located 2.5 cm under the ripple crests. Looking at Figure 10, hydrodynamic pressure patterns are affected by changes in grain size diameters and porosity. However, a general periodicity of the hydrodynamic pressure distribution due to the presence of bedforms can also be seen in all three flow environments. Hydrodynamic pressures fluctuate under neutral conditions between 0.2 and 0.9 Pa, under losing conditions between 3 and 8 Pa and under gaining conditions between 0.3 and 0.8 Pa. Generally, hydrodynamic pressure peaks under the mid part of the luv sides of the ripples and reaches minimum beneath the troughs between the ripples.
Under neutral, gaining, and losing conditions, differences between hydrodynamic pressures and local pressure oscillations result from different ambient groundwater flow conditions as well as from sediment heterogeneity. To investigate these effects more precisely, flow velocities resulting from hydrodynamic pressures are illustrated for different grain size diameters in the next section.

Flow Patterns
Different hyporheic flow vectors result from the hydrodynamic pressure distribution. Looking at the flow patterns shows the effects of gaining/losing conditions for different porosities and grain size diameters. Figure 11 displays groundwater flow velocities greater than 1 × 10 −3 m/s, 1 × 10 −4 m/s, and 1 × 10 −5 m/s and their directions for neutral, gaining, and losing conditions. These are superimposed on grain diameters map in the sediment underneath ripples 4, 5, and 6. Higher velocity in a zone is an indicator that it is a zone of preferred flow. Unlike neutral and gaining conditions, with v > 1 × 10 −3 m/s, a slight downward flow across ripple surfaces is detected under losing conditions which indicates that under such conditions larger volumes of water enter the sediment compared with gaining and neutral conditions. Looking at v > 1 × 10 −4 m/s, some hyporheic flow occurs under neutral and gaining conditions. Gaining conditions display slight upward flow through the hyporheic sediment. In this velocity class and for losing conditions, not only strong hyporheic flow occurs but also intense downward flow. Through downward pulling of water, hyporheic zone becomes deeper in losing conditions compared with neutral and gaining conditions. Even with velocities as low as 1 × 10 −5 m/s, there is not much flow through the sediment under neutral conditions and the hyporheic flow is limited to the top 5-6 cm sediment. Ambient flow conditions on the other hand cover the entire sediment at these flow velocities.
In Figure 11, the effect of sediment heterogeneity on flow under different conditions can be investigated by tracking the flow over grain size diameters and porosities (see also Figure 3). Under neutral conditions, where no upward/downward ambient groundwater flow exists, water majorly enters the streambed from the luv side of the ripples and leaves the streambed from the lee (downstream) side of the ripples. The extent in which this flow can enter the sediment is controlled by the grain size diameter distribution in the sediment. For v > 1 × 10 −4 m/s, the hyporheic flow of the fourth ripple is allowed to enter deeper into the sediment compared with the fifth and the sixth ripple. This is due to the flow being blocked by the blue stretch of lower grain size diameter (as well as lower porosity) under the fifth and the sixth ripple. At less intense (v > 1 × 10 −5 m/s) flow under neutral conditions, a blue block underneath the fifth ripple and partly the sixth ripple has limited hyporheic flow to the upper 5-6 cm as opposed to the fourth ripple in which flow can enter deeper into the sediment. Under the fourth ripple, hyporheic flow which reaches deeper sediment is deflected by the less porous/smaller grain size diameter areas and instead of flowing back to the surface water, partly continues downward flow. Therefore, as opposed to the conditions in fifth and sixth ripples, for hyporheic NGWA.org V.S. Gollo et al. Groundwater 61, no. 5: 721-732 flow through the fourth ripple longer residence times occur.
As was previously mentioned, hyporheic zone is thicker under losing conditions compared with gaining and neutral conditions. Under losing conditions and at higher flow rates (v > 1 × 10 −4 m/s), deflected hyporheic flow under the fourth ripple is being pulled down by the downwelling groundwater flow. By entering deeper groundwater, the flow meets new heterogeneous zones and is further deflected causing residence times to increase drastically. The hyporheic flow beneath the fifth and sixth ripple is being pulled down left (under fifth ripple) and right (under sixth ripple) of the blue block ( Figure 11). The low grain size diameter/low porosity blue block under the fifth ripple causes the hyporheic flow to partly take a long path of entering the sediment from its left side. A part of the flow takes the shorter path across the blue block. As the result, a portion of the flow resides longer in the sediment compared with the other portion that leaves immediately by flowing towards the lee side of the ripple. Under gaining conditions and at high velocities (v > 1 × 10 −4 m/s), water enters the streambed sediment less deep and behaves similar to neutral conditions. At less intense flow (v > 1 × 10 −5 m/s), upwelling groundwater flow that finds its way up through higher grain size diameters/higher porosity zones, opposes the inflowing hyporheic water in some areas and pushes the outflowing hyporheic flow towards the lee sides of the ripples. A better connectivity of the upwelling flow towards the fourth ripple due to higher grain size diameter/higher porosity underneath this ripple, results in more impact of gaining conditions on the hyporheic flow compared with the fifth ripple. Investigating the flow across rippled heterogeneous streambeds under different ambient groundwater flow conditions shows that a base hyporheic zone generated under neutral conditions becomes deeper by downwelling groundwater flow (losing conditions). Under upwelling groundwater flow (gaining conditions), the depth in which hyporheic flow occurs is decreased. Sediment heterogeneity strongly affects flow under all conditions, that is, neutral, gaining, and losing. Compared with neutral conditions, under losing and gaining conditions the impact is stronger and more complex as the ambient groundwater flow prefers zones of higher grain size diameter/higher porosity. This impact is stronger under losing conditions compared with gaining conditions as downward pulling of flow entraps the hyporheic water in deeper heterogeneous zones. This entrapment which happens by the deflection of flow through complex heterogeneous zones can occur to a lesser extent under neutral and gaining conditions as well. If occurring near the sediment-water interface, flow reflection decreases hyporheic flow path length by redirecting flow towards surface water. Therefore, depending on the base hyporheic flow rates (generated under neutral conditions), sediment heterogeneity, and ambient groundwater flow conditions hyporheic flow residence times can be modified. For investigating such complex phenomenon, the integral solver which can detect flow redirections for heterogeneous sediment under different ambient flow conditions seems very competent.

Conclusion
An integral solver was further developed to model tracer transport through heterogeneous rippled streambeds and thereby gain deep insights into hyporheic flow processes in the sediments. For this purpose, an integral model was compared with the results measured on a physical model, that is, a laboratory flume. Even in the presence of groundwater flow and with complex streambed heterogeneity, the modeled tracer propagation agrees well with the experiment. This confirms the applicability of the model. A look at the resulting flow fields suggests that hyporheic zone thickens under losing conditions more than under neutral and gaining conditions. Sediment heterogeneity can increase the hyporheic flow domain volume. Consequently, with velocity in the stream channel being constant, increased sediment heterogeneity in the porous medium (resulting in thicker hyporheic flow domain) causes higher residence time of the water when comparing flow through sediment underneath ripples. In zones with higher porosity/larger grain size diameter, stronger water flow occurs than in zones with lower porosity/smaller grain size diameter. Different sediment types lead to complex flow patterns and deflect the flow, also towards deeper sediment layers. In part, water is also trapped in low flow zones, which further increases residence times in the hyporheic zone. Small-scale, high-resolution integral modeling of hyporheic flow through heterogeneous sediment provides the necessary information on changes in residence times and flow pathways that are crucial for determining the purification effect and fate of contaminants in the hyporheic zone. In this case, the application of the integral approach should be preferred over the prevailing coupling schemes which neglect the mutual, continuous feedback between groundwater and surface water (Sobhi Gollo et al. 2022).
The integral solver has already been validated to account for streambed morphology  and transport through homogeneous sediment under ambient groundwater conditions . Transport through heterogeneous sediment under ambient groundwater conditions has now been discussed and validated in the present study. The integral approach has so far proven to be very powerful in modeling flow and transport processes across groundwater-surface water interface. Further development of the integral solver by including a wider range of scenarios, transient flow conditions, reactive transport and sediment transport are highly recommended.