A comparative study of integral and coupled approaches for modeling hydraulic exchange processes across a rippled streambed

Although both are crucial parts of the hydrological cycle, groundwater and surface water had traditionally been addressed separately. In recent decades, considering them as a single hydrological continuum in light of their continuous interaction has become well established in the scientific community through the development of numerous measurement and experimental techniques. Nevertheless, numerical models, as necessary tools to study a wide range of scenarios and future event predictions, are still based on outdated concepts that consider groundwater and surface water separately. This study compares these “coupled models”, which result from the successive execution of a surface water model and a groundwater model, to a recently developed “integral model”. The integral model uses a single set of equations to model both groundwater and surface water simultaneously, and can account for the continuous interaction at their interface. For comparison, we investigated small-scale flow across a rippled porous streambed. Although we applied identical model domain details and flow conditions, which resulted in very similar water tables and pressure distributions, comparing the integral and coupled models yielded very dissimilar velocity values across the groundwater–surface water interface. These differences highlight the impact of continuous exchange across the interface in the integral model, which imitates such flow processes more realistically than the coupled model. A few decimeters away from the interface, modeled velocity fields are very similar. Since the integral model and the surface water component of the coupled model are both CFD-based (computational fluid dynamics), they require very similar computational resources, namely access to cluster computers. Unfortunately, replacing the surface water component of the coupled model with the widely used shallow water equations model, which indeed would reduce the computational resources required, produces inaccuracy.


Introduction
On the surface of the planet, groundwater and surface water are in continuous interaction through all surface water types, such as streams, lakes, and wetlands, in many different terrains, from mountains to oceans (Winter et al. 1998;Toran 2017). In the hydrologic continuum, groundwater and surface water are connected through the groundwater-surface water interface (Sophocleous 2002;Bobba 2012). Although groundwater and surface water should be considered a single resource due to their continuous interaction (Winter et al. 1998), they are generally modeled separately due to different characteristics and accessibilities (Brunke 2001;Toran 2017). In recent years, interest in understanding the importance of surface water and groundwater interactions and their integrity as a single continuum has significantly increased (Bobba 2012;Toran 2017;Coluccio and Morgan 2019,). The development and utilization of measurement methods and experiments to determine groundwater-surface water exchange processes has further accentuated their significance (Kalbus et al. 2006;Coluccio and Morgan 2019). Groundwater-surface water interaction happens on the local spatial scale (within the streambed) as well as on large spatial scales (affected by the whole catchment or watershed) (Larkin and Sharp 1992;Brunke and Gonser 1997;Woessner 2000). On the local scale, major morphological features such as meanders (Boano et al. 2006;Revelli et al. 2008;Gomez et al. 2012), sediment bars (Tonina and Buffington 2007;Boano et al. 2010;Marzadri et al. 2010), pool-riffle sequences (Harvey and Bencala 1993;Tonina and Buffington 2009;Trauth et al. 2014), and ripples (Cardenas and Wilson 2007a, 2007b, Bottacin-Busolin and Marion 2010, Fox et al. 2014) all trigger groundwater-surface water interaction by significantly modifying local pressure and shear stresses. These processes occur in the hyporheic zone, where 10% of the groundwater is induced from the surface water (Harvey and Bencala, 1993). Groundwater-surface water interaction influences the fate and transport of contaminants (Conant et al. 2004;Chapman et al. 2007;Kalbus et al. 2007) and impacts the ecosystem of streams (Brunke and Gonser 1997). The water quality of both groundwater and surface water are affected by their interaction (Stanford and Ward 1988;Edwards 1998;Fraser and Williams 1998;Hill et al. 1998;Hayashi and Rosenberry 2002;Fleckenstein et al. 2010;Mojarrad et al. 2019). Various biological activities and the distribution of aerobic and anaerobic conditions within groundwater as well as the fate and transport of waterborne substances such as nutrients and organic carbon (Jones and Holmes 1996;Mulholland et al. 1997;Storey et al. 1999;Stonedahl GEM -International Journal on Geomathematics (2022) 13:16 Page 3 of 27 16 et al. 2012) are also controlled by these interactions (Harvey et al. 2013;Marzadri et al. 2011;Zarnetske et al. 2011). Compared to field and laboratory experiments, numerical models are often less computationally demanding and can provide information for a broader range of scenarios and future predictions. They can also investigate parameters that are hard to measure with physical measurement methods. The first coupling of groundwater-surface water models was presented in the 1990s (Kollet and Maxwell 2006;Jones et al. 2008) to investigate precisely this groundwater-surface water interaction. In a coupled model, the shared boundary between the two models (a groundwater and a surface water model), which is the surface water-groundwater interface, is used to transfer parameter information from one model to another. The most common coupling technique is one-way sequential coupling with pressure as the coupling parameter (Tonina and Buffington 2009;Trauth et al. 2015). This technique is widely used and has been validated for a variety of groundwater-surface water interaction cases, including localscale flow and exchange across streambeds (Saenger et al. 2005;Cardenas and Wilson 2007a;Jin et al. 2010; Bardini et al. 2012;Janssen et al. 2012;Trauth et al. 2013Trauth et al. , 2014Trauth et al. , 2015Chen et al. 2018). Saenger et al. (2005) coupled shallow water equations (for surface water) with Darcy's law (for groundwater) to model various surface water flow rates in a riffle-pool sequence. Their results showed that higher flow rates increase hyporheic exchange and reduce residence times. Cardenas and Wilson (2007a) and Janssen et al. (2012) coupled the RANS equations (Reynolds Averaged Navier-Stokes equations, for surface water) and Darcy's law to investigate the interaction between turbulent flow and groundwater exchange. These models are capable of offering simple flow predictions. Jin et al. (2010) coupled the RANS equations with Darcy's law to study the transport of non-sorbing solutes in a streambed with periodic bedforms. They concluded that these transport processes in the groundwater were advectiondominant. Trauth et al. (2015) investigated the exchange processes in a pool-riffle morphology by coupling the Navier-Stokes equations with the Richards equations (Richards 1931). They stated that it was necessary to use the Navier-Stokes equations to account for turbulence in the surface water. The Navier-Stokes equations are solved using CFD (Computational Fluid Dynamics) modeling software such as OpenFOAM (Open Field Operation and Manipulation;Weller et al. 1998). Open-source software such as OpenFOAM has further facilitated the development of coupling tools such as hyporheicFoam (Li et al. 2004), which allows researchers to couple the RANS equations and Darcy's law.
When studying local groundwater-surface water interaction, the paradigm has shifted towards investigating groundwater-surface water and their interaction as a single resource. In spite of the use of novel measurement and experimental techniques, however, coupled modeling techniques still treat groundwater and surface water as two separate domains, and neglect their continuous spatiotemporal interaction. A decade ago, Oxtoby et al. (2013) presented an alternative approach to the integral modeling of surface water and groundwater flow. The integral approach allows both surface water and groundwater to be modeled using a single set of equations. Oxtoby et al. (2013) generated the porousInter solver (described in Sect. 2.3) of OpenFOAM to manage this approach. porousInter is an extension of the widely used interFoam solver, which uses the Navier-Stokes equations to model multi-phase fluid flow processes. The porous-Inter extension includes porosity and grain diameter as parameters over the region that characterize a porous medium in the model. Due to the high computational requirements of the integral approach, however, its use is limited to cases at the local scale. A local-scale physical setup by Fox et al. (2014), who used a flume experiment to investigate flow and tracer propagation across a homogenous rippled streambed (a similar setup to the rippled domain of this study, which is explained below) was modeled by Broecker et al. (2021) using this integral solver. The results showed a very good agreement between Broecker et al. (2021)'s simulations andFox et al. (2014)'s experiments. In local-scale investigations, the integral approach has proven to be capable of determining high-resolution continuous flow and exchange across the groundwater-surface water interface. Although the ability of this approach to detect continuous interaction across the groundwater-surface water interface appears to be superior to the coupled approach, a systematic comparison between coupled and integral approaches is still lacking. High-resolution modeling via the integral approach utilizes precise mapping of flow and pressure; however, this is what increases its computational requirements. Computational burden is crucial in defining the limitation of the integral approach.
In this study, we therefore aimed to model open channel flow across a rippled porous streambed ( Fig. 1) with integral and coupled approaches to highlight the plausibility of using the integral approach as a step towards modeling groundwater and surface water in a single continuum. Such a study is relevant for local-scale processes, such as hyporheic exchange. For the coupled model (CM), we chose the widely used one-way sequential coupling of groundwater and surface water via pressure. For the groundwater component of the coupled model (GW-CM), we used PCSiWaPro ® (based on the  Guo et al. 2017; PC Sickerwasserprognose (seepage flow prognosis); see Sect. 2.2). For the surface water component of the coupled model (SW-CM), we used the interFoam OpenFOAM solver (based on the Navier-Stokes equations; see Sect. 2.1.1). Throughout many local-scale coupling studies (e.g., Saenger et al. 2005), shallow water equations have been used as the SW-CM. While the use of the Navier-Stokes equations to account for turbulence in the surface water has been recommended (Trauth et al. 2015), up to this point, no comparisons between the use of shallow water equations and the Navier-Stokes equations to determine the impacts of model simplification (see Sect. 2.1.2) on model accuracy and reliability have been made. The selection of the proper surface water modeling equations will determine the computational requirements of the CM. This selection of equations was a major contribution to the primary objective of this study, which was the systematic comparison of the CM and the integral model (IM). In this study, we compared the flow fields and the computational demands of the IM and the CM to answer the scientific questions described in Sect. 1.3.
To demonstrate the capability of the solver for the integral model (porousInter) when modeling groundwater flow, two cases (seepage through dam) have been verified by employing a comparison to PCSiWaPro ® results (used as the groundwater component of the coupled model) as well as several other numerical and analytical solutions (see the "Appendix").
By systematically comparing an integral model and a coupled model in the flow and exchange across a rippled streambed, we aim to answer the following questions: • How do flow fields over the entire domain appear for both models? What consequences does this have for future modeling? • How different are flow fields in the interface region adjacent to surface water/groundwater? What causes these differences? What consequences does this have for future modeling? • How different are the computational requirements of the two models? • In conclusion, which approach is preferable for future modeling of the groundwater-surface water interaction?

Surface water modeling
In this study, we first modeled the surface water flow over of the rippled streambed (SW-CM; Fig. 1) using the interFoam solver of the OpenFOAM software (based on the Navier-Stokes equations; see Sect. 2.1.1). We also investigated the possibility of model simplification by modeling the same domain using shallow water equations (Sect. 2.1.2).

Navier-Stokes equations
We used OpenFOAM version 2.4.0, an open-source computational fluid dynamics program, to simulate free surface flow over the streambed (SW-CM) as well as to simulate the IM (see Sect. 2.3). To model the free surface flow, we used the interFoam solver, which is applicable for multiphase fluid flows. The two-phase (water, air) solver has been widely applied in hydraulic engineering in recent years, with the air phase mainly included to account for the movement of the free water surface (Schulze and Thorenz 2014, Higuera et al. 2014, Schmitt et al. 2015, Bayon et al. 2016 where − → υ [m/s] is the flow velocity vector and t [s] is time. Conservation of mass (Eq. (2)) and momentum (Eq. (3)) are thus written as follows: where ρ [kg/m 3 ] is the density of the fluids, p [Pa] is the pressure, μ [m 2 /s] is the dynamic viscosity (μ = μ ph + μ t ; physical + turbulent viscosity), and g [m/s 2 ] is the vector acceleration of gravity. Besides laminar flow, three types of turbulence models are included in OpenFOAM. One is RANS (Reynolds-averaged Navier-Stokes) models, which include turbulence models such as k-ε, k-ω, and k-ω SST. The others are LES (Large Eddy Simulations) and DNS (Direct Numerical Simulations) models. Compared to RANS models, DNS and LES are more advanced turbulence models. Here, the LES model was applied to predict turbulence within the fluid (see OpenFOAM turbulence guide documentation).
Simulating the free surface shallow water flow using solvers based on the Navier-Stokes equations such as interFoam is computationally demanding. The Navier-Stokes equations can be simplified to shallow water equations in cases where vertical flow as well as deviations from hydrostatic pressure are not important. This simplification reduces the computational time significantly and is explained in the following section.

Shallow water equations
Shallow water equations are widely used for river and stream flow. Shallow water equations are derived from depth-integration of the Navier-Stokes equations. As a consequence, there is no vertical velocity and the pressure distribution in the vertical direction is hydrostatic. In addition to the Navier-Stokes equations, we also used shallow water equations to model the SW-CM.
Conservation of mass and momentum without considering sink or source terms (e.g., precipitation or infiltration) for two-dimensional depth-averaged shallow water equations are written as follows: is the bottom elevation, μ t [m 2 /s] is turbulent kinematic viscosity, and n[s/m 1/3 ] is Manning's roughness coefficient.
To model shallow water, we used the Java-based hms (hydroinformatics modeling system; compare Simons et al. 2014) modeling framework. hms makes it possible to calculate water levels and depth-averaged velocities in consideration of Manning's roughness coefficient.

Groundwater modeling
The groundwater component of the coupled model (GW-CM) was simulated using PCSiWaPro ® . This software is used for unsaturated and saturated groundwater flow modeling. It was developed by IBGW Leipzig (Ingenieurbüro für Grundwasser GmbH) and Technische Universität Dresden. The properties of the sediment were defined according to DIN 4220 (2008). The DIN is the well-established German National Organization for Standardization; DIN 4220 is the German standard for the identification, classification, and derivation of soil parameters. PCSiWaPro ® is based on the Richards equation (Eq. (5)) and uses the wetting and rewetting curves from Luckner et al. (1989) to estimate unsaturated hydraulic properties. In the following equation the extended Darcy's law is introduced in the mass conservation equation: where θ [m 3 /m 3 ] is the volumetric water content, K (θ ) [m/s] is the hydraulic conductivity, and ψ [m] is the pressure head.

Integral model
For the IM, we chose a sub-solver of interFoam called porousInter (Oxtoby et al. 2013). To account for incompressible fluids and porous mediums, the equations for the conservation of mass (Eq. 2) and momentum (Eq. 3) are written as modified Navier-Stokes equations; [] f indicates an averaged parameter over a void region: where ϕ [−] is the effective porosity and D [kg/(m 2 s 2 )] is an additional porous drag term.
where d p [m] is the effective grain size diameter. The term "A" in Eq. (9) describes the pressure loss of the fluid with porous medium due to friction in line with Ergun (1952). In porous medium, compared to free flow, more momentum is needed to accelerate a given volume of water, which is referred to as "added mass"; this is accounted for through "B" ( van Gent 1995). In areas where only free flow exists, effective porosity is set to 1, which results in A = B = D = 0 and allows the use of the original Navier-Stokes equations for free water flow (Eqs. 2 and 3). Similarly, for the SW-CM (see Sect. 2.1.1), for the IM, the LES turbulence model was used.

Comparable hydrogeological conditions assessment
The sediment parameters ϕ and d p are required in porousInter when modeling a flow in a porous medium. PCSiWaPro ® follows the DIN 4220 standard, which classifies sediment according to sediment type, θ s [m 3 /m 3 ] (saturated volumetric water content), θ r [m 3 /m 3 ] (the residual volumetric water content), and K 0 [m/s] (matching hydraulic conductivity at saturation). For consistency with natural streambed sediment (as discussed in the experiment conducted by Fox et al. (2014)), the sediment type derived from DIN 4220 (for the GW-CM) was "pure sand" with a saturated volumetric water content (θ s or ϕ) equal to the porosity chosen in porousInter. When the system is fully saturated, residual volumetric water content is not relevant. For porousInter, a particle size (d p ) of 2 mm and a porosity (ϕ) of 0.25 was chosen. The saturated hydraulic conductivity (K 0 ) used in PCSiWaPro ® was derived from Eq. (11) (following Hazen) using the same particle size chosen for porousInter: where d 10 is a 10% fall-through of the sieve curve [mm].
Since porousInter is limited to a single uniform particle size, d 10 = 2 mm and, with Eq. (11), K 0 = 4.64 × 10 -3 m/s was used in this study.

Computational domain
In Fig. 1, light blue indicates the air phase, dark blue the water phase, and light brown the sediment, including respective boundary conditions. Figure 1a shows the model setup of the surface water flow using hms (shallow water equations model) for single-phase flow over a rippled streambed. We used this model to investigate if model simplification was applicable to the SW-CM (see Fig. 1b, top). The CM is shown in Fig. 1b, where the top depicts the SW-CM (to be modeled using interFoam; Navier-Stokes equations model) and the bottom highlights the GW-CM (to be modeled using PCSiWaPro ® ). Regarding the surface water and air phase, the IM shown in Fig. 1c is identical to that of Fig. 1b, top, and longer than the domain of the GW-CM in its sediment part. Modeling the GW-CM with a groundwater domain of the same length as the one from the IM is not necessary for our purposes. In our study, groundwater flow is only generated due to the presence of the ripple morphology (induced by surface water, which only flows in an x direction). Flow in the flat streambed of the GW-CM ("far" away from the ripples) is therefore negligible. Setting the model boundaries 2 m from both sides of the rippled area is sufficient to ensure full formation of the flow pathways on the studied area of 3 m length with 15 ripples, as shown in Figs. 1c and d.
Model setups of the SW-CM and the IM were created in three dimensions and had a 1 m width along the y axis, 15 m along the x axis, and a maximum of 1.5 m along the z axis. To be able to model turbulence using LES, we set up a 3D domain with a 1 m width on the y axis. Another reason we chose LES was its ability to better capture the formation of small eddies between the ripples in comparison with simpler and faster turbulence models such as k-ω and k-ε. A Smagorinsky sub-grid scale model with a van Driest dumping function was used to reduce the eddy viscosity in the nearwall region. This facilitates the reproduction of the characteristics of direct numerical simulations, which solve the three-dimensional Navier-Stokes equations for all eddies directly, at the near-wall region.
Coupling parameters were transferred to the GW-CM from results of pressure distribution at a y = 0.5 m cross section above the ripples of the SW-CM.

Meshes, initial and boundary conditions, and parameters
This study utilizes unstructured meshes, which were generated using the GMSH threedimensional finite element mesh generator (for the IM and the SW-CM). The SW-CM (Fig. 1b, top) is meshed identically to the surface water part of the IM. The GW-CM was generated using the meshing tool in PCSiWaPro ® using mesh resolutions. Cells along the interface of the surface water and pore water domain coincide for easy transfer of results. The shallow water model (Fig. 1a) is 2D in the x-y plane. With no flow in the y direction, it is practically a 1D illustration of flow in the x direction, consisting of 15,000 cells (cell size = model resolution = 10 -3 m). The IM and the SW-CM consist of more than 430,000 and 116,000 prismatic cells, respectively. The minimum cell size (model resolution; shortest edge size) was 5 × 10 -3 m near ripples compared to a maximum size (longest edge size) of 3.5 × 10 -1 m in the air phase far above the water-air interface for both the SW-CM and the IM. The GW-CM has 15,000 cells (model resolution 5 × 10 -3 m). Model convergence was achieved using these resolutions.
Water flow with a value of 0.5 m 3 /s was set as the boundary condition at the inlet on the left side in Fig. 1 along the x axis for all three models. The top pressure boundary of the GW-CM (yellow line) was defined by taking the pressure distribution from the results of the steady-state run of the SW-CM; left and right boundaries (blue lines) were defined as the linear hydrostatic pressure increase according to depth. Black lines indicate no-flow boundaries with normal velocity being zero, while green lines indicate atmospheric pressure.
No sediment transport was assumed for any cases. We ran both the SW-CM and the IM with steady boundary conditions until a quasi-steady state was achieved. A full steady state would be reached when parameters such as flow velocity and pressure were constant over time. Under turbulent flow conditions, however, such a state can never be fully satisfied. For this reason, we derived initial flow and turbulence conditions after some precomputational time, as soon as the oscillations of the results were relatively small. In this study, we defined quasi-steady state conditions as the point when turbulent eddies between ripples became fully formed (although still moving) and no inconsistent changes (any irregular and unexpected jump in a variable's value) in pressure fields, flow fields, or eddy sizes/shapes occurred. After 60 s of precomputation, these quasi-steady states for the IM and the SW-CM were reached. However, we calculated variable values as averages between 60 and 300 s to account for oscillations within the adjusted range. The shallow water model achieved a full steady state after 300 s. The GW-CM was also modeled starting from a full steady state.
The IM, the SW-CM, and the shallow water model are designed to maintain a0 .5 m water table above the streambed. We defined a 0.5 m water table above the streambed initially for all models so that it was possible to reach this condition faster. We placed a weir in all models to maintain this waterhead. The streambed was initially saturated for the GW-CM and the IM, and all the initial velocities were zero. To account for friction in the shallow water model, we set the Manning coefficient to 0.03 s/m 1/3 (see Sect. 2.1.2), which is the coefficient recommended for clean and straight natural streams.
In the next section, we present the simulation results of the IM and the CM by analyzing the flow velocities in the water above and the sediment underneath the investigated 8th ripple, shown in Fig. 1. Next, we dedicate additional attention to flow fields adjacent to the interface. These were investigated by mapping the hydrodynamic pressures and the flow fields across the interface. sides of theThe next step was to consider the effects of simplifying the models by replacing the SW-CM with a shallow water model.

IM and CM flow velocity fields in the model domain
To investigate the flow velocities in the model domain, we chose a representative slice (8 th ripple in Fig. 1; 735 cm < x < 760 cm, 50 cm < z < 50 cm in Fig. 2) consisting of the water and the sediment parts. The air part was only simulated for precisely capturing the surface water level fluctuations; as this is not relevant here, we excluded it from the representative slice.
To compare the results, we exhibited the simulated flow velocities across several vertical and horizontal cross sections. Vertical cross sections intersect with the luv (x = 748 cm) and the lee (x = 758 cm) sides of the 8th ripple and extend through the sediment and the water parts (Fig. 2). We placed the horizontal cross sections in the water (z = 40 cm), the sediment (z = 40 cm), and the interface-adjacent sediment (z = 2 cm). These cross sections are slightly shifted to the left (to cover the area from the 7th ripple crest to the 8th ripple crest). Figure 3 illustrates the flow velocities of both models (v C M and v I M ) through vertical cross sections at x = 748 cm ( Fig. 3a and b) and x = 758 cm ( Fig. 3c   (z = 2 cm), flow velocity differences between the two models reach as high as 2 × 10 -3 m/s. For deeper groundwater (Fig. 5, z = 40 cm), these differences are almost zero in the z direction (v z ) and peak at 2.5 × 10 -5 in the x direction (v x ).
As displayed in Fig. 6, maximum flow velocity differences between the IM and the CM (for both v x and v z ) for surface water at z = 40 cm are 7 × 10 -3 m/s.
In general, flow velocities in the surface water are higher than in the groundwater. In addition, for the current case, flow in the groundwater is only induced by the surface water flow across a rippled streambed. This causes higher groundwater velocities near the interface compared to the deeper groundwater. A direct comparison of the maximum flow velocity differences between the IM and the CM over the cross sections at z = 2 cm, 40 cm and 40 cm would therefore be inaccurate and misleading. Instead, we divided the averaged (over the horizontal cross section) flow velocity differences   Table 1. D C M and D I M ratios demonstrate the flow velocity differences between the two models considering their deviation from the expected flow velocity simulated by each model (CM and IM). For both v x and v z , across z = 2 cm, these ratios were one to three orders of magnitude greater than those of the z = 40 cm and z = 40 cm cross sections. This means that CM/IM model discrepancy is higher near the interface at the z = 2 cm cross section compared to other horizontal cross sections. D C M and D I M values at z = 40 cm and z = 40 cm are relatively smaller, which means that both models function more similarly a few decimeters away from the groundwater-surface water interface.
We therefore determined that major flow-velocity discrepancies exist between the two models (CM and IM) near the interface. We will further discuss the meaning of these differences in the discussion about flows across interface-adjacent surface water (Sect. 4.2).

IM and CM pressure and velocity fields in the interface-adjacent zone
By maintaining a 0.5 m water level over the sediment in the whole domain throughout the entire simulation period, the hydrostatic pressures of both models remain equal. The surface water flow was also defined identically for the surface water component of the coupled model (SW-CM) and the surface water part of the integral model. The simulated hydrodynamic pressures (p_rgh, Fig. 7) of the two models on the 8th ripple (see Fig. 1 above), however, differ from each other on the ripple surface. These differences are displayed in Fig. 7b. Another indicator of the differences between the two models near the interface is shown in Fig. 8. Higher flow velocities and bigger turbulent eddies can be seen within the troughs between the ripples of the CM (Fig. 8b) compared to those of the IM (Fig. 8a). In Fig. 8a, flow vectors are formed along and through the interface. These paths of flow display a continuous exchange between the surface water and the groundwater, and they show how surface water is transported through the sediment matrix from one side of the ripple to the other. Further information about the groundwater can be gained by looking at the groundwater component separately (Fig. 9).  The patterns and hotspots of the flow differ when comparing the interface-adjacent flow in the groundwater (Fig. 9). We also realized that a zone on the crest of the ripple in the IM (the area highlighted green in Fig. 9) has a velocity of 0.02 m/s < v I M < 0.2 m/s, which for this case is outside of the area where Darcy's law would apply (Reynolds number (−) > 10).
The results emphasize that near the groundwater-surface water interface, the IM and the CM behave very differently. The fundamentals of these differences and their significance in selecting a plausible small-scale groundwater-surface water modeling approach is further discussed in Sect. 4.

Shallow water simplification
Here we examine the plausibility of using a simpler set of equations to model the SW-CM, which could drastically decrease the computational requirements. For this purpose, we have chosen the hms solver, which solves shallow water equations (see Sect. 2.1.2). Using shallow water equations, v z is zero and v x values along a vertical cross section are constant. Nevertheless, we generated two parabolic velocity profiles  Fig. 3. In shallow water equations, pressure is assumed to be hydrostatic, which conflicts with the hydrodynamic pressure results displayed in Fig. 7 for the SW-CM. Following the discussions of the IM and the CM flow and exchange in the next section, we discuss the reduction of computational requirements by model simplification using the shallow water equations model.

Required computational resources
Modeling the SW-CM via hms requires only 30 min of runtime on a PC (8 core 3.4 Ghz AMD processors). However, to simulate the SW-CM based on the Navier-Stokes equations, access to a computer cluster is required. Using 40 processors at the highperformance computing cluster of Technische Universität Berlin, it takes about 24 h for the SW-CM to complete the simulation. Using the same resources, it takes about 70 h for the IM to run completely. For the GW-CM, about one hour of computing with the aforementioned PC is sufficient.

Discussion
In the introduction (Sect. 1.2), we described the coupled model studies that are most relevant to the domain of our study (flow across a morphologically modified streambed). In these studies, groundwater was modeled with Darcy's law (or the Richards equations-i.e., relatively saturated Darcy's law equations). The surface water component of the coupled model was either simulated with the Navier-Stokes equations or shallow water equations. By modeling the GW-CM with Darcy's law (Richards equations in saturated sediment are equivalent to Darcy's law) and the SW-CM with both the Navier-Stokes equations and shallow water equations, our coupled model is representative of the development of previous models. Here, we discuss this representative coupled model alongside the recently developed integral modeling approach. For the discussion, we divide the comparison of the flow fields of the IM and the CM into two parts: • The area a few decimeters away from the interface.
• The area near the interface.
We use the Navier-Stokes-based SW-CM for our initial discussions. In the next step, we discuss simplifying the SW-CM approach by using shallow water equations. Then, we compare the computational requirements of the CM to the IM.

Flow comparison a few decimeters away from the interface
Darcy's law, and the Richards equation based on it (as used for the GW-CM via PCSiWaPro ® ), is the standard mathematical approach to determining flows in a porous medium. The IM investigated here, however, takes a different approach to determining flow in this zone (see Sect. 2.3). The appendix contains results showing the successful verification of the integral solver against several analytical and numerical solutions (including PCSiWaPro®) for two cases simulating seepage through a dam. Incidentally, this solver has also been validated for the case of flow triggered by bioturbation in benthic zones in comparison to a model based on Darcy's law ). Nevertheless, for flow across a morphologically modified (rippled) streambed, the current study is a further step in validating the integral approach. In Fig. 3, we can see that when starting from the interface (z~0), flow velocities of the models in both flow directions (v x and v z ) converge as depth increases. It also shows that flow velocities correspond in deep groundwater (z = 40 cm).
For the surface water part of the domain, both models use the same set of flow simulation equations (Navier-Stokes equations) and the meshes are identical. However, the models' dissimilarities starting from the interface are propagated upwards. By around z = 40 cm, flow velocity results indicate that both models behave similarly.
Both models behave similarly a few decimeters from the groundwater-surface water interface. This shows that the mathematical solver of the IM is capable of plausibly simulating the flow processes, especially in the groundwater, where its equations are distinct from Darcy's law. Nevertheless, the notable source of the difference between the two models-the interface-adjacent zone-requires further discussion.

Flow comparison near the interface
Here, we separately discuss the groundwater and surface water flow near the groundwater-surface water interface.

Interface-adjacent surface water
Looking at Fig. 8, it is clear that velocity vectors meet the interface (body of the ripple on the luv side) at different angles. In the IM, velocity vectors are bent upon the incident, and depending on the impact angle, they either enter the groundwater or are deflected back into the surface water. The velocity of the water that enters the sediment is dampened. The deflected flow vector (if deflected towards the foot of the lee side) creates a turbulent eddy in the trough between the ripples. In the CM, there is no path through the sediment as the interface is a no-flow boundary. As a result, the only option of the surface water flow vectors that meet the interface is to be deflected back into the surface water. Here, all the flow vectors that will be deflected and trapped in the trough between the ripples join together and create larger eddies compared to the ones in the IM. Accumulation of these larger eddies in the CM result in higher flow velocities in the interface-adjacent surface water compared to the IM. The CM (SW-CM) velocity distributions in the interface-adjacent area define the hydrodynamic pressures on the interface. In Fig. 7, the hydrodynamic pressure of the CM is therefore greater than that of the IM. Fox et al. (2014) and (2016) physically modeled a very similar flow across a rippled streambed setup. In these experiments, no sudden jump of flow as simulated by the SW-CM could be detected. Roche et al. (2018) studied wave propagation from the surface water into the groundwater. The absence of high-velocity zones near the interface was also confirmed by the integral model of this study.
Our results show that the assumption of a no-flow boundary at the interface leads to an overestimation of the flow at the interface adjacent to the surface water. In turn, by allowing the surface water to continuously interact with the groundwater, the IM yields more plausible results of the flow processes near the groundwater-surface water interface.

Interface-adjacent groundwater
Pressure (hydrostatic and hydrodynamic pressure) as the coupling parameter for the GW-CM was acquired from the above-mentioned SW-CM velocity distribution across the interface. The interface in the IM is included in the model domain, and therefore pressure values are constantly readjusted in light of the flow. The velocity distribution of the models across the ripple in Fig. 9 is therefore different. The formation of one large eddy in the trough between ripples (see Fig. 8) causes maximum hydrodynamic pressure on the foot of the ripple, where the maximum groundwater flow velocity of the GW-CM is detectable. In contrast to this, allowing continuous flow movement from the surface water into the groundwater via the IM results in an accumulation of flow hotspots in various areas across the body of the ripple. In addition, flow velocities outside of the boundaries of Darcy's law can be seen in the IM. Due to the relatively small grain size considered in this study, these flows are limited to the top sediment (distance from the interface < 1 cm). However, for larger grain diameters, many studies (e.g., Blois et al. 2014;Roche et al. 2018) have shown that in a groundwater-surface water continuum, surface water turbulence can be transported into the groundwater and create areas of water velocities approaching those of surface water. Such phenomena cannot be simulated via a groundwater model based on Darcy's law (such as the GW-CM). Distinctive numerical solvers and discretization methods of the groundwater component of the CM and the IM could potentially impact the results. However, we think that this impact is negligible in this study, as we have checked the grid convergence for both models.
Our findings demonstrate that in the area near the surface water-groundwater interface, the IM yields more plausible results compared to the CM. In the area near the interface, where the flow fields of the two models are distinctive, the advantages of the IM are particularly stark compared to the CM. Due to the similarity of the models further from this zone, however, the computational requirements should also be discussed when choosing the appropriate model.

SW-CM simplification using a shallow water model
When modeling the surface water using the shallow water equations model (hms), flow in the vertical direction (v z ) is assumed to be zero. Nevertheless, we have shown that in the presence of actual streambed morphology, multidirectional flow (by accumulation of eddies in the troughs between ripples) is generated. In addition to generating velocity in the vertical direction (v z ), eddies modify horizontal velocities (v x ) and pressure fields. A model simplification via the shallow water equations for the SW-CM is inappropriate for the following reasons: • Exclusion of the vertical velocities; • Oversimplification of the horizontal velocities (see Fig. 10); • Assumption of hydrostatic pressure.
It was essential to use a model based on the Navier-Stokes equations to model the SW-CM of this study. This corroborates the statement by Trauth et al. (2015), who claimed that due to the presence of turbulent eddies across morphologically modified streambeds, shallow water equations could not be used to model surface water flow.

Computational-demand comparison
Both the CM and the IM require access to computer clusters. Although the runtimes of the IM and the SW-CM differed by around 46 h in our case, we state that the necessity to employ significant computational capacity, which is a key parameter in CFD modeling, is the same for both models.
When comparing the GW-CM to the groundwater part of the IM, it should be noted that very large domains will increase the computational requirements of the IM substantially. We have demonstrated how both models function similarly for flows far from the interface. Therefore, where the effect of the interface is negligible, a Darcy's law groundwater model, which is computationally much less demanding than the IM, is sufficient.

Conclusions
Groundwater-surface water modeling approaches should be adjusted to reflect the current paradigm shift that conceptualizes groundwater and surface water as a single resource. Specifically, unlike current modeling approaches, the surface watergroundwater interface, which is the hotspot of biogeochemical interaction, should be accurately represented in the simulated domain. To pursue this, in our study we exhibited the plausibility of using a single set of equations (modified Navier-Stokes equations) governing small-scale flow processes in both groundwater and surface water. We achieved this by comparing our "integral approach" to the well-established approach of one-way sequential coupling via pressure in the case of flow across a rippled streambed. As a result of the integral model's ability to reflect the continuous interaction between groundwater and surface water, it proved to represent flows close to the water-sediment interface more plausibly than the coupled model could.
Our study advocates for the use of the integral model for questions requiring a very exact simulation of local (small)-scale processes in the interface domain within a few decimeters of the interface into the surface water and the groundwater. The model domain of this study presented such an interface domain, where surface watergroundwater interaction was induced by ripples; this is one of the main drivers of hyporheic zone processes (Lewandowski et al. 2020). Furthermore, we determined that high-resolution modeling of "deep" (more than a few decimeters) groundwater is computationally too demanding when using the integral approach, and a standard groundwater model would be sufficient for these cases.
The integral approach can be applied to a wide range of scenarios, such as transient flow conditions, non-Darcy flows in porous media, and combined free and porous media flow around breakwaters and dikes. article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.