Document [original]

Integral modeling approach for flow and transport

at surface water-groundwater interfaces

Master of Science

Vahid Sobhi Gollo

ORCID: 0000-0002-5970-2475

an der Fakultät VI - Planen Bauen Umwelt -

der Technischen Universität Berlin

zur Erlangung des akademischen Grades

Doktor der Ingenieurwissenschaften

- Dr.-Ing. -

genehmigte Dissertation

Promotionsausschuss:

Vorsitzender: Prof. Dr.-Ing. Matthias Barjenbruch

Gutachter: Prof. Dr.-Ing. Reinhard Hinkelmann

Gutachter: Dr. Jörg Lewandowski

Gutachter: Prof. Dr.-Ing. Ilhan Özgen-Xian

Gutachterin: Prof. Dr.-Ing. Nicole Saenger

Tag der wissenschaftlichen Aussprache: 13. Januar 2023

Berlin 2023

Preface

I would like to give special thanks to my supervisor Prof. Dr. Reinhard (Phillip) Hinkelmann,

for giving me a chance to be part of the Urban Water Interfaces Research Training Group and

for accepting me as his PhD student. I would like to thank him dearly for his constructive

feedbacks, patience and scientific guidance throughout my PhD candidacy. I would like to

thank my supervisors Dr. Jörg Lewandowski and Prof. Gunnar Nützmann for the fruitful

feedbacks and discussions. I express my gratitude to Dr. Tabea Broecker for her help and

productive collaboration on the research topic. I would like to thank all my former colleagues

from the Chair of Water Resources Management and Modeling of Hydrosystems at TU Berlin

as well as all members of the UWI Research Training Group for all the good memories and the

positive impact they had on my scientific career as well as personal growth.

Here I would like to appreciate Ralf Duda, a colleague and a friend with unmatching patience

and good mood who has always supported me through technical issues. I am likewise very

thankful for the administrational support and the encouragement by Dr. Gwendolin Porst and

Tosca Piotrowski.

My biggest gratitude goes to my wife Aisan for all the love and support during this time.

Berlin, January, 2023

Acknowledgments

This contribution was elaborated within project H2 of the interdisciplinary Research Training

Group “Urban Water Interfaces” (UWI, RTG 2032/1+2) funded by the Deutsche

Forschungsgemeinschaft (DFG) located at the Technische Universität Berlin (TUB) and the

Leibniz-Institute of Freshwater Ecology and Inland Fisheries (IGB).

For modeling in sections 2 and 3, high performance computers (HPC) of TU Berlin were used.

The origin and the foundation of the software PSCiWaPro® used in section 2 is a result of the

BMBF research project “Sickerwasserprognose“. The use of the software was permitted with

an agreement between Technische Universität Dresden and Technische Universität Berlin.

For modeling in section 4, North-German Supercomputing Alliance (HLRN) were used. The

experimental data in section 4 were provided by Prof. Arnon from the Ben-Gurion University

of the Negev.

iii

Publications of cumulative doctoral thesis

First author publications:

1. Sobhi Gollo, V., T. Broecker, J. Lewandowski, G. Nützmann, and R. Hinkelmann.

2021. An integral approach to simulate three-dimensional flow in and around a

ventilated U-shaped chironomid dwelled burrow, Journal of Ecohydraulics,

https://doi.org/10.1080/24705357.2021.1938258.

2. Sobhi Gollo, V., Broecker, T., Marx, C., Lewandowski, J., Nützmann, G. &

Hinkelmann R. 2022. A comparative study of integral and coupled approaches for

modeling hydraulic exchange processes across a rippled streambed, International

Journal on Geomathematics 13, no 16: 1-27. https://doi.org/10.1007/s13137-022-

00206-5.

3. Sobhi Gollo, V., Broecker, T., Lewandowski, J., Nützmann, G. & Hinkelmann R. 2022.

Flow and transport modeling in heterogeneous sediments using an integral approach.

Groundwater. https://doi.org/10.1111/gwat.13275.

Conference contributions:

1. Sobhi Gollo, V., T. Broecker, J. Lewandowski, G. Nützmann, and R. Hinkelmann. 2019.

Comparison of an integral solver to coupled modelling approach for hydraulic exchange at

surface water-groundwater interface along rippled streambed, Vol. 21, EGU (European

Geosciences Union) 2019-7759-1.

2. Sobhi Gollo, V., T. Broecker, J. Lewandowski, G. Nützmann, and R. Hinkelmann. 2019.

Integral flow and transport modelling approach for surface water-groundwater interface

domain, GQ 19 (Groundwater Quality), Liege, Belgium.

3. Sobhi Gollo, V., T. Broecker, and R. Hinkelmann. 2021. High-resolution integral modelling

approach for flow and transport in groundwater - surface water interaction space. SIAM

Conference on Mathematical and Computational Issues in the Geosciences (Online)

Supplementary contributions

Co-authored publications:

1. Broecker, T., K. Teuber, V. Sobhi Gollo, G. Nützmann, J. Lewandowski, and R.

Hinkelmann. 2019. Integral Flow Modelling Approach for Surface Water-Groundwater

Interactions along a Rippled Streambed. Water, 11(7), 1517, DOI: 10.3390/w11071517.

2. Broecker, T., V. Sobhi Gollo, A. Fox, J. Lewandowski, G. Nützmann, S. Arnon, and R.

Hinkelmann. 2021. High Resolution Integrated Transport Model for Studying Surface

Water–Groundwater Interaction. Groundwater, DOI: 10.1111/gwat.13071.

An overview of all supplementary scientific work is given in section 5.

The thesis was carried out as project H2 of the Research Training Group ‘Urban Water

Interfaces‘(GRK 2032/1+2) and was funded by the German Research Foundation (DFG).

Abstract

Groundwater and surface waters are crucial parts of hydrological cycle which traditionally have

been treated separately. Continuous interaction between these two water bodies have attracted

requisite attention towards considering them as a single hydrogeochemical continuum in recent

decades. Many studies have addressed these interactions on large scales such as catchments as

well as local scale hyporheic zone interactions. Influenced by factors such as streambed

morphology, ambient groundwater conditions, sediment heterogeneity and bioturbation,

hyporheic zone is the hotspot of retention, transformation and attenuation of solutes as well as

habitat of a variety of aquatic organisms. Through development of novel measurement methods

and experimental techniques, investigating groundwater and surface water as a single

hydrologic unit is now very well established in the scientific community. Nevertheless,

numerical models as necessary tools to study wide range of scenarios and future event

predictions are still based on concepts that consider groundwater and surface water separately.

These so called “coupled models” result from successive execution of a surface water model

and a groundwater model, oftentimes without any feedback effects. Continuous feedback across

groundwater-surface water interface is a key component in successfully investigating small

(local) scale intense exchange processes such as hyporheic flow and bioturbation. As opposed

to coupled modeling approach, in the current study small-scale groundwater-surface water

processes using a novel integral approach is discussed which allows the continuous feedback

across groundwater-surface water interface. Due to high computational effort of this approach,

only small-scale high-resolution scenarios in the close vicinity of the interface such as

hyporheic flow through rippled streambed and bioturbation (pumping activity of tube-dwelling

organisms) are investigated.

In the first step, for a case of flow and exchange across a rippled streambed the integral approach

is compared to the widely used one-way sequential coupling approach. For integral modeling,

porousInter solver of OpenFOAM computational fluid dynamics model is applied solving

Navier-Stokes equations which are extended by resistance terms including porosity and grain

size diameter parameters to account for flow in porous medium. For coupling, next to the

groundwater model “PCSiWaPro®”, OpenFOAM is used for surface water modeling. The

necessity of using Navier-Stokes equations-based OpenFOAM for surface water modeling is

discussed next to a simpler shallow water equations model called “hms”. Integral to coupled

model comparison shows that under turbulent surface water flow in the vicinity of the ripples,

continuous feedback of flow between surface water and hyporheic water is simulated with the

integral model leading to reduced turbulence compared to the coupled model and flow and

pressure fields near ripples are affected. Results of the two models are very similar far away

from the interface in deeper sediment (> 40 cm), however, significantly differ above the

interface in the surface water.

In a further step the area of application of the integral approach is expanded towards another

important small-scale groundwater-surface water exchange phenomenon. Here, using the

integral approach, the pumping activity of tube-dwelling macroinvertebrates in a U-shaped

burrow and its impact on hydrodynamic exchange processes between porewater and overlying

(surface) water are modeled and compared to experiments and coupled model (Brand et al.

2013). It was realized that unlike coupled model, by using integral model continuous feedback

between burrow water and surrounding porewater resulted in modification of pressure and

velocity fields in the vicinity of the burrow. As experimental approaches are not capable of

direct measurement of flow and exchange processes in the burrow and its surrounding

porewater, the integral approach offers a novel possibility about quantifying continuous

porewater-burrow. In the first two steps, integral approach was compared to widely used

coupled models. Comparison results were plausible and showed that not only integral approach

can model a variety of small-scale high-resolution groundwater-surface water interaction

processes, but also by allowing continuous feedback between two water bodies, could offer

valuable insight towards expanding the understanding of groundwater-surface water

interactions.

In the final step, the integral solver “porousInter”, which was formerly tested for transport

through homogenous sediment, is further applied for tracer transport through heterogenous

sediment. The solver is verified by comparing the modeled transport through heterogenous

sediment to analytical solutions. It is then validated by simulating the flow and tracer transport

across rippled heterogenous sediment experiments of Fox et al. (2016). Simulated tracer

propagation agreed well with the experiments. The simultaneous effect of ambient groundwater

conditions and sediment heterogeneity on hyporheic flow and residence times is discussed.

The main outcomes of this thesis are validation of the integral approach, discussion of its

advantages in comparison to other groundwater-surface water interface modeling approaches,

extension of its area of applicability in modeling small-scale groundwater-surface water

interaction induced by pumping activity of tube-dwelling macroinvertebrates in a U-shaped

burrow and further application for modeling tracer transport through heterogenous sediment

under ambient groundwater conditions. Integral approach can be used for investigations of

engineered hyporheic zones and more general for river and lake sediment ecological questions

vii

as well as a component in integrated water resources management. This approach can further

be applied to free surface and porous media flow interaction processes such as flow through

porous breakwaters and the simultaneous overtopping of and seepage through dikes and dams.

viii

Zusammenfassung

Grundwasser und Oberflächenwasser sind wichtige Komponenten des hydrologischen

Kreislaufs, die bisher getrennt voneinander behandelt worden sind. Die kontinuierliche

Interaktion zwischen diesen beiden Wasserkörpern hat in den letzten Jahren die notwendige

Aufmerksamkeit geweckt, um sie als ein einzelnes hydrogeochemisches Kontinuum zu

betrachten. Viele Studien haben sich mit Wechselwirkungen auf großer Skala (z.B.

Einzugsgebiete) sowie mit den Wechselwirkungen in der hyporheischen Zone auf lokaler Skala

befasst. Beeinflusst von Faktoren wie der Morphologie des Flussbettes, den Bedingungen des

umgebenden Grundwassers, der Heterogenität der Sedimente und der Bioturbation, ist die

hyporheische Zone der Hotspot für die Retention, die Transformation und den Abbau von

Stoffen sowie der Lebensraum für eine Vielzahl von Wasserorganismen. Durch die

Entwicklung neuartiger Messmethoden und experimenteller Techniken ist die Untersuchung

von Grund- und Oberflächenwasser als ein einziges hydrologisches Kontinuum in der

wissenschaftlichen Gemeinschaft inzwischen sehr etabliert. Allerdings basieren die

numerischen Modelle, die für die Untersuchung einer Vielzahl von Szenarien und die

Vorhersage künftiger Ereignisse erforderlich sind, immer noch auf den Konzepten, die Grund-

und Oberflächenwasser getrennt betrachten. Diese sogenannten "gekoppelten Modelle"

ergeben sich aus der aufeinanderfolgenden Ausführung eines Oberflächenwassermodells und

eines Grundwassermodells, oft ohne Rückkopplungseffekte. Eine kontinuierliche

Rückkopplung an der Schnittstelle zwischen Grund- und Oberflächenwasser ist eine

Schlüsselkomponente für die erfolgreiche Untersuchung intensiver Austauschprozesse auf

kleiner (lokaler) Skala wie hyporheische Strömung und Bioturbation. Im Gegensatz zu

gekoppelten Modellierungsansätzen werden in der vorliegenden Studie kleinskalige

Grundwasser-Oberflächenwasser-Prozesse mit einem neuartigen integralen Ansatz betrachtet,

der eine kontinuierliche Rückkopplung über die Grundwasser-Oberflächenwasser-Grenzfläche

ermöglicht. Aufgrund des hohen Rechenaufwands dieses Ansatzes werden nur kleinskalige,

hochaufgelöste Szenarien in unmittelbarer Nähe der Grenzfläche untersucht, wie z.B.

hyporheische Strömungen durch ein gerippeltes Flussbett und Bioturbation (Pumpaktivität von

Wurmgang bewohnenden Organismen).

In einem ersten Schritt wird für den Fall der Strömung und des Austauschs durch ein gerippeltes

Flussbett der integrale Ansatz mit dem häufig verwendeten einseitigen sequenziellen

Kopplungsansatz verglichen. Für die integrale Modellierung wird der porousInter-Löser des

OpenFOAM-Modells verwendet, der die Navier-Stokes-Gleichungen löst, die um

Widerstandsterme erweitert werden, die von der Porosität und dem Korngrößendurchmesser

abhängen, um die Strömung in einem porösen Medium zu berücksichtigen. Zur Kopplung wird

neben dem Grundwassermodell "PCSiWaPro®" auch OpenFOAM für die Modellierung von

Oberflächenwasser verwendet. Die Notwendigkeit der Verwendung des auf den Navier-Stokes-

Gleichungen basierenden OpenFOAM für die Oberflächenwassermodellierung wird neben

einem einfacheren Flachwassergleichungsmodell "hms" diskutiert. Der Vergleich zwischen

dem integralen und dem gekoppelten Modell zeigt, dass bei einer turbulenten

Oberflächenwasserströmung in der Nähe der Rippel eine kontinuierliche Rückkopplung

zwischen dem Oberflächenwasser und dem hyporheischen Wasser mit dem integralen Modell

simuliert wird, was im Vergleich zum gekoppelten Modell zu einer geringeren Turbulenz führt

und die Strömungs- und Druckfelder in der Nähe der Rippel beeinflusst. Die Ergebnisse der

beiden Modelle ähneln sich weit entfernt von der Grenzfläche im tieferen Sediment (> 40 cm),

unterscheiden sich aber deutlich oberhalb der Grenzfläche im Oberflächenwasser.

In einem weiteren Schritt wird der Anwendungsbereich des integralen Ansatzes auf ein weiteres

wichtiges kleinskaliges Grundwasser-Oberflächenwasser-Austauschphänomen erweitert. Hier

wird mit Hilfe des integralen Ansatzes die Pumpaktivität von Wurmgang bewohnenden

Makroinvertebraten in einem U-förmigen Wurmgang und deren Auswirkung auf

hydrodynamische Austauschprozesse zwischen Porenwasser und dem darüber liegendem

(Oberflächen-)Wasser modelliert und mit Experimenten und gekoppelten Modellen (Brand et

al. 2013) verglichen. Es wurde festgestellt, dass im Gegensatz zum gekoppelten Modell bei der

Verwendung des integralen Modells eine kontinuierliche Rückkopplung zwischen dem

Grabenwasser und dem umgebenden Porenwasser zu einer Veränderung der Druck- und

Geschwindigkeitsfelder in der Umgebung des Wurmgangs führt. Da experimentelle Ansätze

nicht in der Lage sind, die Strömungs- und Austauschprozesse im Wurmgang und dem sie

umgebenden Porenwasser direkt zu messen, bietet der integrale Ansatz eine neuartige

Möglichkeit zur Quantifizierung des kontinuierlichen Porenwassers im Wurmgang. In den

ersten beiden Schritten wurde der integrale Ansatz mit häufig verwendeten gekoppelten

Modellen verglichen. Die Vergleichsergebnisse waren plausibel und zeigten, dass der integrale

Ansatz nicht nur eine Vielzahl von kleinskaligen, hochaufgelösten Grundwasser-

Oberflächenwasser-Wechselwirkungsprozessen modellieren kann, sondern auch durch die

Ermöglichung einer kontinuierlichen Rückkopplung zwischen zwei Wasserkörpern wertvolle

Erkenntnisse zur Erweiterung des Verständnisses von Grundwasser-Oberflächenwasser-

Wechselwirkungen liefern kann.

Imletzten Schritt wird der integrale Löser porousInter, der zuvor für den Tracer-Transport durch

homogenes Sediment getestet wurde, auch für den Tracer-Transport durch heterogenes

Sediment angewendet. Der Löser wird durch den Vergleich des modellierten Transports durch

heterogenes Sediment mit analytischen Lösungen überprüft. Anschließend wird er durch die

Simulation von Strömung und Tracer-Transports durch gerippeltes heterogenes Sediment und

Vergleich mit den Experimenten von Fox et al. (2016) validiert. Die simulierte Tracer-

Ausbreitung stimmte gut mit den Experimenten überein. Dabei wird das Zusammenwirkender

Umgebungsbedingungen des Grundwassers und der Heterogenität des Sediments auf den

hyporheischen Fluss und Aufenthaltszeiten diskutiert.

Die wichtigsten Ergebnisse dieser Arbeit sind die Validierung des integralen Ansatzes, die

Diskussion seiner Vorteile im Vergleich zu anderen Ansätzen zur Modellierung der

Grundwasser-Oberflächenwasser-Grenzfläche, die Erweiterung seines Anwendungsbereichs

zur Modellierung kleinskaliger Grundwasser-Oberflächenwasser-Wechselwirkungen, die

durch die Pumpaktivität von Wurmgang bewohnenden Makroinvertebraten in einem U-

förmigen Wurmgang erzeugt werden, und die weitere Anwendung zur Modellierung des

Tracer-Transports durch heterogenes Sediment im Grundwasser. Der integrale Ansatz kann für

Untersuchungen von ingenieurtechnischgebauten hyporheischen Zonen und allgemeiner für

ökologische Fragen im Zusammenhang mit Fluss- und Seesedimenten sowie als Komponente

der integrierten Wasserbewirtschaftung verwendet werden. Dieser Ansatz kann auch für

Interaktionsprozesse zwischen freien Oberflächenwasser und Wasser in porösen Medien

angewandt werden, wie z.B. die Strömung durch poröse Wellenbrecher und das gleichzeitge

Überlaufen und Durchsickern von Deichen und Dämmen.

Contents

Preface i

Acknowledgments ii

Publications of cumulative doctoral thesis iii

Abstract v

Zusammenfassung viii

1. Introduction 1

1.1 Background 1

1.2 Groundwater-surface water interaction measurement 1

1.3 Spatial scales of groundwater-surface water interaction 3

1.4 Hyporheic zone 6

1.5 Concepts of modeling groundwater-surface water interactions 9

1.6 Flow and transport modeling with OpenFOAM 11

1.6.1 Modeling surface water flow and transport 12

1.6.2 Turbulence models 13

1.6.3 Integral modeling of groundwater-surface water interactions 14

1.7 Coupled modeling of groundwater-surface water flow 15

1.8 Scope of this thesis 16

1.8.1 Current research deficits 16

1.8.2 Aim of this thesis 17

2. Comparison of integral approach to coupled approach for flow across rippled

streambed 19

2.1 Abstract 20

2.2 Introduction 20

2.3 Materials and methods 24

2.3.1 Surface water modeling 24

2.3.2 Groundwater modeling 26

2.3.3 Integral model 27

2.3.4 Comparable hydrogeological conditions assessment 27

2.3.5 Computational domain 28

2.3.6 Meshes, initial and boundary conditions and parameters 29

2.4 Results 31

2.4.1 IM and CM flow velocity fields in the model domain 31

2.4.2 IM and CM pressure and velocity fields in the interface-adjacent zone 37

2.4.3 Shallow water simplification 39

2.4.4 Required computational resources 40

2.5 Discussion 40

2.5.1 Flow comparison a few decimeters away from the interface 41

2.5.2 Flow comparison near the interface 41

2.5.3 SW-CM simplification using a shallow water model 43

2.5.4 Computational-demand comparison 43

2.6 Conclusions 44

2.7 Appendix 45

xii

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow

using integral approach 46

3.1 Abstract 46

3.2 Introduction 47

3.3 Materials and methods 50

3.3.1 Computational domain 50

3.3.2 Governing equations 51

3.3.3 Grid and boundary conditions 53

3.3.4 System parameters 54

3.3.5 Modeling cases 55

3.4 Results and discussion 56

3.4.1 Comparison of integral to coupled approach 56

3.4.2 Comparison of flow patterns to experiment results 62

3.5 Conclusion 64

4. Flow and transport modeling in heterogenous sediments using integral approach 65

4.1 Abstract 66

4.2 Introduction 66

4.3 Materials and methods 68

4.3.1 Computational domain 68

4.3.2 Boundary and initial conditions 69

4.3.3 Governing equations and computational tool 70

4.4 Results and discussion 72

4.4.1 Validation of flow and transport in heterogenous sediment with a 1D

analytical solution 72

4.4.2 Comparison of numerical to experimental results for tracer transport 73

4.4.3 Pressure and flow distribution in heterogenous sediment 77

4.5 Conclusion 82

5. Supplementary contributions 83

5.1 Integral Flow simulations for groundwater-surface water interactions at rippled

streambeds 83

5.2 Transport simulations for groundwater-surface water interactions under neutral,

losing and gaining flow conditions for homogenous streambed 84

6. Synthesis 86

6.1 Summary and conclusions 86

6.1.1 General outcomes 86

6.1.2 Outcomes of comparing integral to coupled approach for flow across rippled

streambed 88

6.1.3 Outcomes of comparing integral to coupled approach for ventilation through

U-shaped burrow 89

6.1.4 Outcomes of transport modeling for heterogenous groundwater-surface

water- interactions 90

6.1.5 Recommendations for modeling and data processing 91

6.2 Outlook 92

xiii

Appendix A: Validation of porousInter for dispersive transport 94

Appendix B: Test case overview 98

B1 Seepage through a dike modeled with porousInter and PCSiWaPro® 98

B2 Seepage through a rectangular dam modeled with porousInter and

PCSiWaPro® 101

B3 One-dimensional one-phase modeling of surface water flow over rippled

streambed with hms (shallow water equations) 104

B4 Three-dimensional two-phase modeling of surface water flow over rippled

streambed with interFoam 105

B5 Two-dimensional modeling of groundwater flow in sediment under rippled

streambed with PCSiWaPro® 107

B6 Three-dimensional two-phase integral modeling of flow across rippled

streambed with porousInter 108

B7 Three-dimensional flow simulation in and around a ventilated U-shaped

burrow using porousInter – case 1 110

B8 Three-dimensional flow simulation in and around a ventilated U-shaped

burrow using porousInter – case 2 112

B9 One-dimensional tracer transport in heterogenous groundwater with

porousInterTracer 114

B10 Transport simulations for groundwater-surface water interactions at

heterogenous rippled streambeds 116

Bibliography 119

xiv

List of Figures

Figure 1-1: Direct-real time oxygen concentration measurement setup in hyporheic pore

water placed at the side channel of the river Erpe, Berlin/Brandenburg;

measurement data is collected for the project H1 of the DFG Research Training

Group Urban Water Interfaces; photo by Birgit Maria Müller ......................................... 2

Figure 1-2: Groundwater-surface water exchange types: a) gaining condition, b) losing

condition, c) neutral condition, after Winter et al. (1998) ................................................ 4

Figure 1-3: Hyporheic exchange through (1) meander, (2) dry bar, (3) pool-riffle, (4)

ripple and (5) bioturbation, cross section “a” across flow direction ................................. 5

Figure 1-4: Major hyporheic drivers and processes after Lewandowski et al. (2019): (1)

water-sediment interface, (2) synthetic or organic matter (dissolved or particulate)

from urban areas and industry, (3) local and large scale groundwater-surface water

exchange, (4) dynamic biofilm, (5) bioturbation and microbial activity, (6)

macroinvertebrates, (7) protozoa, (8) labile carbon gradients, (9) anoxic

microzones, (10) solute and particle transport, (11) sorption sites, (12) temperature

gradients, (13) surface flow forcing, (14) hyporheic exchange flux, (15)

heterogenous streambed conductivity and morphology, (16) infiltration zones and

(17) exfiltration zones ....................................................................................................... 8

Figure 1-5: a) Sediment components, b) surface and groundwater flow at interface, c)

single domain and d) two domain groundwater-surface water coupling schemes;

modified after Mosthaf et al. (2011) ............................................................................... 10

Figure 2-1: Model setup and boundary conditions of a) the shallow water equations

surface water model; b) the coupled model (CM), including the surface water

component of the coupled model (SW-CM, b, top) and the groundwater component

of the coupled model (GW-CM, b, bottom); c) the integral model (IM); and d) the

ripple geometry. .............................................................................................................. 28

Figure 2-2: Graphic visualization of the ripples at the sediment–water interface and the

horizontal (blue dotted lines) and the vertical (green dotted lines) cross sections

used to compare flow results of the integral model (IM) and the coupled model

(CM). ............................................................................................................................... 31

Figure 2-3: A comparison of flow velocities vx and vz in the x (a, c) and z (b, d) directions

for vertical cross sections at x = 748 cm (a, b) and x = 758 cm (c, d) for the integral

model (IM) and the coupled model (CM); absolute values of subtracting flow

velocities have been calculated using the CM and the IM (|vCM- vIM|) in the x (e)

and z (f) directions for both vertical cross sections. ....................................................... 32

Figure 2-4: A comparison of flow velocities vx and vz in the x (a) and z (b) directions for

the horizontal cross section at z = -2 cm along 735 cm < x < 755 cm for the integral

model (IM) and the coupled model (CM); absolute values of subtracting flow

velocities have been calculated using the CM and the IM (|vCM- vIM|) in the x (c)

and z (d) directions. ........................................................................................................ 33

Figure 2-5: A comparison of flow velocities vx and vz in x (a) and z (b) directions for the

horizontal cross section at z = -40 cm along 735 cm < x < 755 cm for the integral

model (IM) and the coupled model (CM); absolute values of subtracting flow

velocities have been calculated using the CM and the IM (|vCM- vIM|) in the x (c)

and z (d) directions. ........................................................................................................ 34

Figure 2-6: A comparison of flow velocities vx and vz in the x (a) and z (b) directions for

the horizontal cross section at z = 40 cm along 735 cm < x < 755 cm for the integral

model (IM) and the coupled model (CM); absolute values of subtracting flow

velocities have been calculated using the CM and the IM (vCM- vIM|) in the x (c)

and z (d) directions. ........................................................................................................ 35

Figure 2-7: a) Hydrodynamic pressure values at the surface water–groundwater interface

on the investigated ripple (representative cross section, 8th ripple) for the IM and

the CM; b) subtraction of the hydrodynamic pressures of the CM from that of the

IM on the investigated ripple. ......................................................................................... 37

Figure 2-8: Flow velocities and directions over and through the investigated ripple (8th

ripple, Figure 1) for a) the integral model and b) the coupled model. Note: red

dashed line indicates the sediment–water interface; black arrows only indicate

directions, and are not scaled regarding magnitude. Magnitude is shown as a color

scale. ............................................................................................................................... 38

Figure 2-9: Velocity and Reynolds number (in log10 scale) distributions through the

investigated ripple (8th ripple, Figure 2-1) for a) the integral model (IM) and b) the

coupled model (CM). ...................................................................................................... 38

Figure 2-10: A comparison of the horizontal surface water flow (vx) for vertical cross

sections at a) x = 748 cm and b) x = 758 cm for the surface water component of the

coupled model (SW-CM) and the shallow water equations in the surface water

model (SW-SM). ............................................................................................................. 39

Figure 2-11: Model geometry and results of different models and analytical solutions for

seepage through a rectangular dam (left) and a dike (right). ........................................ 45

Figure 3-1: Model domain (lengths are not to scale) and boundaries B1-B4 reported in

Table 1; B1 is atmospheric pressure at cylinder top; B2 are cylinder walls and

cylinder bottom; B3 is symmetry plane to reduce computational effort; B4

represents the pumping chironomid larvae. The left panel is a detail enlargement of

the burrow shown in the right panel ............................................................................... 50

Figure 3-2: Pressure (left) and velocity (right) fields of Case 2 (top) and Case 1 (bottom)

with pumping velocity = 0.009 m/s and dp (effective grain size) = 0.25 mm. Green

and yellow lines: see Figure 3.1 ...................................................................................... 56

Figure 3-3: Pressure in the pumping zone (orange = left branch of the burrow, blue = right

branch of the burrow) for different water depths ............................................................ 57

Figure 3-4: Maximum/minimum pressure at the pumping zone for 10 different cases of

U/dp; U is the maximum velocity at the pumping zone and dp is the effective grain

size diameter; blue stands for the right side of the baffle where pressure is positive

and orange stands for the left side of the baffle where pressure is negative ................... 58

Figure 3-5: Velocity and pressure profiles at z = -0.05 m through inlet and outlet branches

of the burrow for dp (effective grain size) = 0.25 mm, Case 1 ....................................... 59

Figure 3-6: Velocity magnitudes (colors) and flow directions (arrows) a) in the entire

domain, b) around pumping zone and c) between inlet and outlet branches .................. 60

Figure 3-7: Dominant flow directions (white arrows) on z = -0.05 m plane, looking from

a) top b, c) lateral sides; blue plane is the cross section of a half cylinder .................... 61

Figure 3-8: Comparison of the modeled flow pattern (left) above inlet (top row) and outlet

(bottom row) with results of Morad et al. (2010) (right); dashed lines in modeled

xvi

images mark the sediment surface; note: values differ due to unidentical burrow

lengths and pumping velocities. Images based on the same measurements/raw data

as images shown in Morad et al. (2010)/Roskosch et al. (2010) .................................... 62

Figure 3-9: Velocity at pumping zone to velocity at inlet/outlet (measured at outlet) for 5

different scenarios; dp (effective grain size) = 0.25 mm ................................................ 63

Figure 4-1: Schematic view of hyporheic zone. qh – hyporheic flow from the main channel

into the streambed sediment and returning to the stream, qG – gaining flux from

groundwater to surface water, qL – losing flux from surface water to groundwater. ..... 67

Figure 4-2: Model setup and geometry for modeling flow and transport in heterogenous

sediment based on Fox et al. (2016); sediment heterogeneity is explained in Fig. 4-

3. ..................................................................................................................................... 69

Figure 4-3: (top) Grain diameters and (bottom) porosities. ..................................................... 70

Figure 4-4: Model boundary conditions. .................................................................................. 70

Figure 4-5: Comparison of integral model to analytical solutions for a continuous tracer

injection in a simple 1D model domain: Model domain and porosities (left top),

tracer fronts modelled with integral model after t=400 and 600 s (left bottom) and

comparison of tracer concentrations for different timesteps for the integral model

and the analytical solution. ............................................................................................. 73

Figure 4-6: Tracer propagation for the investigated 6th ripple after 10, 40 and 60 min.

Results shown for both the lab experiment (a) and the integral model (b; purple).

For 40 min, tracer fronts simulated via MIN3P groundwater model (from Fox et al.

2016) are as well illustrated. The ripples of the lab experiments deviated slightly

from the shape at the start of the experiment because of minor sediment movement

during previous experiments. .......................................................................................... 73

Figure 4-7: Contribution of advection and diffusion to transport under different conditions

at 60 min tracer propagation for the 6th ripple. Peclet number = Advective transport

rate/diffusive transport rate. In dark red areas (Peclet number > 1) advection is

predominant. Other areas (Peclet number < 1) are low-turbulent zones where

diffusion is dominant. ..................................................................................................... 74

Figure 4-8: Differences of the real sediment setup in the experiment of Fox et al. (2016)

(top panel) and the originally planned sediment setup which was used in the integral

model (middle panel). The hyporheic zone beneath the 6th ripple, i.e. the discussed

ripple in Fig. 4-9 is magnified (lower panel) to display local sediment packing

differences (el). Note: compared to Fig. 4-2, the first ripple from left is excluded

from the model setup so it matches the image provided by Fox et al. (2016). ............... 75

Figure 4-9: Hydrodynamic pressure distribution under neutral, losing and gaining

conditions and the location of a horizontal cross section (z = 0.195 m, white line)

as well as grain diameter maps at 4th, 5th and 6th ripples................................................. 78

Figure 4-10: Hydrodynamic pressures under neutral, losing and gaining conditions for a

horizontal cross section through the sediment bed in the depth marked in Fig. 4-9

as white line. ................................................................................................................... 78

Figure 4-11: Flow vectors (dimensionless) mapped over grain size diameters of the

sediment underneath 4th, 5th and 6th ripples for v > 1 × 10-3 m/s, v > 1 × 10-4 m/s and

v > 1 × 10-5 m/s for neutral, losing and gaining conditions. ........................................... 80

xvii

Figure A-1: Model domain of 2D dispersive transport validation. Small red square

centered at x = 1 m and y = 5 m shows the location of initial pulse injection with C

= 1 mg/l. Flow is 1D in x-direction. ............................................................................... 95

Figure A-2: Model of dispersive pulse tracer transport progression using porousInter after

1 s, 50 s and 100 s. .......................................................................................................... 95

Figure A-3: Simulated (left) and analytical (right) tracer concentrations after 1, 50 and

100 s over the model domain. Tracer concentration is displayed in the vertical axis.

........................................................................................................................................ 96

Figure A-4: Comparison of the simulated and analytical tracer propagation after 1, 50 and

100 s for cross sections in 𝑥𝑥 and 𝑦𝑦 directions. These cross-sections are highlighted

as red dotted lines in Fig. A-3. ........................................................................................ 97

xviii

List of Tables

|𝑣𝑣𝑣𝑣𝑣𝑣| and |𝑣𝑣𝑣𝑣𝑣𝑣| velocities for z = -2 cm, -40 cm and 40 cm in the 𝑥𝑥 and 𝑧𝑧

directions. ........................................................................................................................ 36

Table 3-1: Boundary locations and conditions for geometry of Figure 1 and parameter

values .............................................................................................................................. 53

Table 4-1: Tracer propagation comparison with time between integral model and

experiment under neutral, losing and gaining conditions. .............................................. 76

Table 4-2: Tracer propagation comparison after 40 minutes between the MIN3P and the

integral model (porousInter) under neutral, losing and gaining conditions. ................... 77

Table B 1: Model setup for test case B1 (seepage through a dike modeled with porousInter

and PCSiWaPro®) .......................................................................................................... 98

Table B 2: Boundary conditions for test case B1 - porousInter ............................................... 99

Table B 3: Boundary conditions for test case B1 – PCSiWaPro®. ....................................... 100

Table B 4: Model setup for test case B2 (Seepage through a rectangular dam modeled

with porousInter and PCSiWaPro®). ........................................................................... 101

Table B 5: Boundary conditions for test case B2 - porousInter ............................................. 102

Table B 6: Boundary conditions for test case B2 – PCSiWaPro® ........................................ 103

Table B 7: Model setup for test case B3 (surface water flow over ripples using hms) .......... 104

Table B 8: Boundary conditions for test case B3 – PCSiWaPro®. ....................................... 104

Table B 9: Model setup for test case B4 (surface water flow over ripples using interFoam)

...................................................................................................................................... 105

Table B 10: Boundary conditions for test case B4 ................................................................. 106

Table B 11: Model setup for test case B5 (surface water flow over ripples using

PCSiWaPro®) ............................................................................................................... 107

Table B 12: Boundary conditions for test case B5 ................................................................. 107

Table B 13: Model setup for test case B6 (Integral modeling of groundwater-surface water

interactions across rippled streambeds using porousInter) ........................................... 108

Table B 14: Boundary conditions for test case B6 ................................................................. 109

Table B 15: Model setup for test case B7 (Integral modeling of U-shaped burrow

ventilation using porousInter - case 1) .......................................................................... 110

Table B 16: Boundary conditions for test case B7 ................................................................. 111

Table B 17: Model setup for test case B8 (Integral modeling of U-shaped burrow

ventilation using porousInter - case 2) .......................................................................... 112

Table B 18: Boundary conditions for test case B8 ................................................................. 113

Table B 19: Model setup for test case B9 (One-dimensional tracer transport in

groundwater with two porosities) ................................................................................. 114

Table B 20: Boundary conditions for test case B9. ................................................................ 115

xix

Table B 21: Model setup for test case B10 (Transport simulations for groundwater-surface

water interactions at heterogenous rippled streambeds). .............................................. 116

Table B 22: Boundary conditions for test case B10. .............................................................. 117

1. Introduction

1.1 Background

At the earth’s scale, hydraulic groundwater and surface water connection is a crucial part of

hydrologic cycle. On the surface of the planet these two major water bodies are in continuous

interaction through all surface water types such as streams, lakes and wetlands in many different

terrains from the mountains to the oceans (Winter et al. 1998). In inland, major groundwater-

surface water interaction areas are wetlands, karst, glaciers, streams, lakes, coastal areas and

polar ice caps (Khan and Khan 2019). In the hydrologic continuum, groundwater and surface

water are connected through groundwater-surface water interface (Sophocleous 2002, Bobba

2012). Surface water-groundwater interactions happen by groundwater lateral flow through

unsaturated sediment (Sophocleous 2002) and by infiltration into or exfiltration from the

saturated zones (Brunke 2001, Sophocleous 2002) hence affecting the quantity and/or quality

of both systems (Bobba 2012). Although due to continuous interaction groundwater and surface

water should be considered a single resource (Winter et al. 1998, Lewandowski et al. 2020),

due to different characteristics and accessibilities they are generally studied separately (Brunke

2001, Lewandowski et al. 2020). In recent years the interest in understanding the importance

of surface water and groundwater interactions and their integrity as a single continuum has

significantly increased (Bobba 2012, Lewandowski et al. 2020). Legislative attention is as well

gained towards integral consideration of groundwater and surface water bodies. Directive

2000/60/EC (34) of the European parliament states that to pursue “the purposes of

environmental protection, there is a need for a greater integration of qualitative and

quantitative aspects of both surface waters and groundwaters, taking into account the natural

flow conditions of water within the hydrological cycle”.

1.2 Groundwater-surface water interaction measurement

Development and utilization of measurements methods and models (see section 1.5) to

determine groundwater-surface water exchange processes has accentuated the significance of

these processes. Next to hydraulic exchange measurements, a major focus has been given to

biogeochemical processes and their significance for aquatic ecology and water quality

(Rosenberry and La Baugh 2008). Many methods are being used to measure several parameters

that are crucial for groundwater-surface water interaction determination. Based on the chosen

method, parameters can be directly or indirectly acquired (Kalbus et al. 2006). Direct parameter

measurements are conducted by performing for instance seepage meter; and temperature profile

1. Introduction

tests to determine seepage flux; and temperature gradients, respectively. Direct oxygen

concentration measurement of hyporheic water (see section 1.4) is illustrated in Figure 1-1.

Small water sampling tubes (small green tubes inside the water) installed in the sediment are

connected to an oxygen meter (blue device on the bridge) that directly measures oxygen

concentration of the sampled water and transmits data to the laptop. Water is pumped through

the tubes using syringe pumps (red devices) connected to the syringes. Indirect parameter

measurements are conducted either using experimental and field methods to determine

hydraulic head, hydraulic conductivity, porosity and groundwater velocity based on Darcy’s

law or using tracer tests, monitoring wells and pumping tests to determine groundwater

components and contaminant concentration based on mass balance approaches (Kalbus et al.

2006, Khan and Khan 2019). Spatial measurement scales of the investigated parameters with

different methods vary from a single point measurement to > 1000 m (Fleckenstein et al. 2010,

González-Pinzón et al. 2015, Magliozzi et al. 2018).

Figure 1-1: Direct-real time oxygen concentration measurement setup in hyporheic pore water

placed at the side channel of the river Erpe, Berlin/Brandenburg; measurement data is collected

for the project H1 of the DFG Research Training Group Urban Water Interfaces; photo by Birgit

Maria Müller

While direct measurements only determine parameters such as seepage flux for a 10-1 m spatial

measuring scale, Darcy’s law-based pumping tests as well as mass balance-based tracer tests

can measure up to 102 – 103 m sections (Kalbus et al. 2006). Seepage meter direct measurements

tools follow a simple concept and are inexpensive but in practice their placement could block

the streamflow and they cannot detect excessive shallow hyporheic flow (Khan and Khan

1. Introduction

2019). Several enhancement tools such as slow pulling of hyporheic water without disturbing

the surrounding sediment using syringe pumps (Figure 1-1) were developed.

Indirect measurement methods although give very good parameter characterization for

measured points, however, can become challenging when dealing with large areas with high

hydrogeologic complexity. Using indirect measurement for parameter characterization, a trade-

off between resolution of heterogeneities and sampled subsurface volume is inevitable (Kalbus

et al. 2006). Although implementing advanced measurements techniques like PIV (Particle

Image Velocimetry, Blois et al. (2014), Roche et al. (2018)) has recently improved parameter

characterization for laboratory experiments, generalization of acquired parameters through

laboratory and especially field measurement campaigns remains challenging. Implementing

different scenarios under varying hydrologic conditions to predict future events and gaining a

deeper understanding of the complex dynamics of the groundwater-surface water exchange is

possibly through modeling studies (Boano et al. 2009, Stonedahl et al. 2012, Brunner et al.

2017, Broecker et al. 2019). The current contribution examines the applicability of the advanced

numerical models to study the groundwater-surface water exchange processes (namely the

integral approach; see section 1.6). In section 1.5, various approaches to model groundwater-

surface water exchange are presented. Next to the methods, it is crucial to discuss the scale of

the groundwater-surface water processes.

1.3 Spatial scales of groundwater-surface water interaction

Groundwater-surface water interaction is a multiscale processes, where local scale effects

aggregate and then manifest across all spatial scales. The spatial scale of a selected method has

a considerable effect on the results of parameter measurement (Kalbus et al. 2006, Khan and

Khan 2019). Groundwater-surface water interactions occur in a variety of spatial scales. A scale

categorization done by Larkin and Sharp (1992), Brunke and Gonser (1997) and Woessner

(2000) is as follows:

1) Large-scale interactions: interaction is affected by whole catchment or watershed

2) Local-scale interactions: within the streambed controlled by e.g. hyporheic flow (see

section 1.4) properties

1. Introduction

For each spatial scale, landscape diversifies physiographic and climatic settings of the area

(Winter et al. 1998). Surface water types such as oceans, streams, lakes and wetlands unify

perspective of the interaction of groundwater and surface water in different landscapes (Winter

et al. 1998).

Figure 1-2: Groundwater-surface water exchange types: a) gaining condition, b) losing

condition, c) neutral condition, after Winter et al. (1998)

For different landscapes and spatial scales, higher/lower groundwater table compared to surface

water table can cause two distinctive flow patterns (Cardenas 2009, Fox et al. 2014, Fox et al.

2016). With a higher groundwater table compared to surface water table, an upwelling

groundwater movement results in groundwater infiltration into the surface water. This state is

referred to as “gaining condition” (Figure 1.2a). Conversely, with a lower groundwater table

compared to surface water, hydraulic gradient is reversed and water exfiltrates from the surface

water into the groundwater. This creates a so called “losing condition” (Figure 1.2b). In a

“neutral condition” water table at the ground- and surface water are in equilibrium (Figure 1.2c;

Silliman & Booth 1993, Woessner 2000).

1. Introduction

Considering the impact of all spatial scales simultaneously is crucial in investigating regional

groundwater-surface water interactions (Ellis et al. 2007; Krause et al. 2007, Mojarrad et al.

2019). Therefore, next to interactions that are triggered by groundwater-surface water table

imbalances, local scale interactions like hyporheic exchange (Merill and Tonjes 2014, Mojarrad

et al. 2019,) and bioturbation (Roskosch et al. 2012, Hupfer et al. 2019, Shrivastva et al. 2021)

should be highlighted.

Figure 1-3: Hyporheic exchange through (1) meander, (2) dry bar, (3) pool-riffle, (4) ripple and

(5) bioturbation, cross section “a” across flow direction

In local scale, large morphological features such as meanders (Figure 1-3 (1); Boano et al. 2006,

Revelli et al. 2008, Gomez et al. 2012), sediment bars (Figure 1-3 (2); Tonina and Buffington

2007, Boano et al. 2010, Marzadri et al. 2010), pool-riffle (Figure 1-3 (3); Harvey & Bencala,

1993; Tonina and Buffington, 2009, Trauth et al. 2014) and ripples (Figure 1-3 (4); Cardenas

& Wilson 2007a, 2007b, Bottacin‐Busolin & Marion, 2010, Fox et al. 2014) trigger hyporheic

flow by significantly modifying local pressure and shear stresses. Hyporheic zone and complex

hydrological process occuring in this zone are further discussed in section 1.4.

Bioturbation is a generic term which is referred to a variety of activities of benthic (Kristensen

et al. 2012) and hyporheic (Shrivastva et al. 2021) organisms that impact streambed/lakebed

sediments. This impact is through particle reworking and burrow ventilation. By mixing the

sediment, freshy settled easily degradable material is transported into deeper sediment and

inorganic material from deeper sediment is introduced to the oxic sediment zones

(Lewandowski et al. 2005a). Burrow ventilation of tube dwelling macroinvertebrates (Figure

1-3 (5)) such as Chironomus Plumosus (Stief et al. 2010b, Hupfer et al. 2019) results in diffusive

and advective aeration of the sediment called bioirrigation (Kristensen et al. 2012) and

significant loss of fixed nitrogen from shallow aquatic ecosystems (Jeppesen et al. 1998,

Sondergaard et al. 2007).

1. Introduction

Groundwater-surface water interaction spatial heterogeneities influence the fate and transport

of contaminants (Conant et al. 2004, Chapman et al. 2007, Kalbus et al. 2007) and impact the

ecosystem of the streams (Brunke & Gonser 1997). Water quality of both groundwater and

surface water are impacted 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). Diverse 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 et al. 2012) are controlled by this interaction (Harvey et al. 2013, Marzadri et

al. 2011, Zarnetske et al. 2011). In dealing with stream water quality issues such as spreading

of industrial, agricultural and sewage contaminants, hyporheic zone, due to its self-purification

effects, is an attractive hydrogeological feature. Its self-purification effect includes nutrient

turnover, degradation of contaminants and the removal of trace organic compounds (Biksey

and Gross 2001, Lewandowski et al. 2019).

1.4 Hyporheic zone

In terms of ecology, geochemistry and hydrogeology “Hyporheic zone” can be defined

differently (Merill and Tonjes 2014, Cardenas 2009). In ecology it refers to biologically active

zone near the stream where groundwater and surface water are interacting (Harvey and Bencala

1993, Harvey and Wagner 2000). Key ecosystem processes such as temperature, primary

productivity, mixing and transport of dissolved oxygen, nutrient cycling, microbe and

invertebrate communities and fish diversity are influenced by these interactions (Brunke and

Gonser 1997, Boulton et al. 2010, Krause et al. 2011). Many freshwater organisms dwell

permanently (Gibert and Deharveng 2002), temporarily, or under certain conditions in

hyporheic zone. Hyporheic zone is a crucial component of the life cycle of many biota such as

fishes, macroinvertebrates and amphibians (Lopez-Rodriguez et al. 2009, Resh and Rosenberg

2010, Williams et al. 2010). In presence of predators or under unfavourable surface conditions

such as floods, droughts or contaminated surface water (Stubbington et al. 2011), some of

organisms (hyporhoes, Boulton et al. 1998, Wood et al. 2010) use hyporheic zone as a

temporary refuge. In areas where excessive deposition of fine sediment as a by-product of

human development deranges the natural hyporheic zone sediment composition and hyporheos

community is negatively affected (Paul and Meyer 2001, Coleman et al. 2011). Under such

conditions, clogging of interstitial space between natural streambed grain blocks stream flow

from downwelling into hyporheic zone (Arnon et al. 2010, Song et al. 2010) which interrupts

1. Introduction

the dissolved oxygen supply and nutrient cycling. Fortunately, presence of (macro)invertebrates

and other organisms can restore permeability in the sediment through bioturbation (Nogaro et

al. 2010).

By increase of solute residence time and solute contact with substrates, hyporheic zone

influences the biogeochemistry of stream ecosystems (Bencala 2000). This impacts stream

water quality by attenuation of contaminants (Biksey and Gross 2001, Conant et al. 2004,

Fischer et al. 2005, Chapman et al. 2007, Kalbus et al. 2007) and through removal of nitrate

and phosphate (Hill and Cardaci 2004, Doyle 2005, Fischer et al. 2005), heavy metals (Feris et

al. 2003a, 2003b), microplastics (Klein et al. 2015), pesticides (Boutron et al. 2011), organic

stormwater contaminants (Peter et al 2019) and pharmaceuticals (Kunkel and Radke 2008, Riml

et al. 2013, Schaper et al. 2018, Schaper et al. 2019).

In hydrogeology, hyporheic zone is defined as sediment surrounding a river in which recent

water from surface flow mixes with groundwater and soon returns to surface water (Harvey and

Bencala 1993, Boano et al. 2014). Next to spatial scales heterogeneity, temporal dynamics of

hyporheic zone are defined with relatively fast transport processes in surface water (Ibisch et

al. 2009, Kennedy et al. 2009) and slower transport in groundwater with velocities within lower

orders of magnitude (Müller 2014). This causes the dynamic exchange of water, substances and

energy and circulations in the hyporheic zone (Krause et al. 2011). Water exchange at the

hyporheic - surface water interface which is caused by local pressure differences and turbulence

(Cardenas and Wilson 2007a) can lead to water level variations and thus pressure gradients

(Elliot and Brooks 1997a, 1997b, Packman et al. 2004). Presence of features like pool-riffle or

ripples in streambed (Figure 1-3) create high pressures on luv side and lower pressure on lee

side of the feature thus allowing surface water infiltration from luv side and exfiltration from

lee side (Elliot and Brooks 1997a, 1997b, Broecker et al. 2021). Next to that, aquifer properties

(Brunke and Gonser 1997), ambient groundwater (Cardenas and Wilson 2007a) and turbulence

in surface water and its propagation (Roche et al. 2018) into the sediment influence hyporheic

exchange. Heterogeneity which is caused by transport and settling of sediment particles over a

variety of flow conditions (Powell 1998), is the main driver in flow adjustment by redirecting

flow through areas with higher hydraulic conductivities (Vaux 1968). Many studies (e.g. Freeze

and Cherry 1979, Pryschlak et al. 2015) have shown that due to its major impact on hyporheic

exchange processes in streams, considering heterogeneity of the streambed has as much of a

priority as stream morphology and groundwater flow conditions. Therefore, the best

1. Introduction

hydrogeological representation of hyporheic exchange processes is achieved by considering

heterogeneity alongside morphology and groundwater conditions (Figure 1-4 Hydrology).

Figure 1-4: Major hyporheic drivers and processes after Lewandowski et al. (2019): (1) water-

sediment interface, (2) synthetic or organic matter (dissolved or particulate) from urban areas

and industry, (3) local and large scale groundwater-surface water exchange, (4) dynamic

biofilm, (5) bioturbation and microbial activity, (6) macroinvertebrates, (7) protozoa, (8) labile

carbon gradients, (9) anoxic microzones, (10) solute and particle transport, (11) sorption sites,

(12) temperature gradients, (13) surface flow forcing, (14) hyporheic exchange flux, (15)

heterogenous streambed conductivity and morphology, (16) infiltration zones and (17)

exfiltration zones

Intensive material turnover and dynamic water interaction in hyporheic zone drives hydro-

biogeochemical processes from many different disciplines (Figure 1-4, Lewandowski et al.

2019). An improved awareness of processes which is achieved though interdisciplinary

investigations enables these processes to be considered more holistically (Lawrence et al. 2013).

A combination of spatially and temporally resolved field and laboratory data with physically

based numerical models is crucial in understanding dynamic processes at the groundwater-

surface water interface (Tonina and Buffington 2007, Cardenas 2010, Cuthbert et al. 2010, Derx

et al. 2010, Westhoff et al. 2010). In the following modeling concepts and tools to describe

surface water-groundwater interactions are discussed.

1. Introduction

1.5 Concepts of modeling groundwater-surface water interactions

For surface water and groundwater investigations, numerical models are widely applied

techniques to study hydrodynamic and hydrogeochemical processes by computational

simulation of different scenarios. Understanding of these complex flow, transport and

transformation processes is achieved using mathematical models that are designed to describe

physical laws and equations. Compared to field and laboratory experiments, they are often less

expensive and can provide information for a broader range of scenarios and future predictions

and investigate parameters that are hard to measure with physical measurement methods

(section 1.2). Groundwater-surface water interaction concepts were first incorporated into

groundwater models in the early 1980s (e.g. Prudic 1989). Early applications were limited

mostly to large scale water quantity management which resulted in coarse model resolutions in

time and space (Gorelick 1986). Surface water flow processes were simplified and aquifer

homogeneity was assumed (Fleckenstein et al. 2010). Later efforts to include ecological

functions (Stanford and Ward 1988) and biogeochemical processes (Hill 1990) were made by

presentation of transient storage models. These advances have led to the development of various

advanced hydrological models that operate at the catchment (large) scale (e.g. Parflow by

Ashby and Falgout 1996, Amanzi/ATS by Moulton et al. 2012, OpenGeoSys by Kolditz et al.

2012). With more focus being given to the local-scale processes (Elliot and Brooks 1997a,

Packman and Brooks 2001), existing groundwater models were linked to more complex surface

water models (Swain and Wexler 1993) and coupled groundwater-surface water models were

presented (Kollet and Maxwell 2006, Jones et al. 2008).

Coupling of flow and transport between groundwater and surface water is achieved by either

modeling of flow through pore spaces or single/two-domain strategies. Pore network models

are based on modeling the flow using e.g. Navier-Stokes equations with considering the

sediment matrix as well as water and air components (Figure 1-5 a). These models are used to

improve the understanding of phenomena such as impact of pore geometry on microbial

communities (Young and Ritz 2000), relation of pore water network to residence times (Bijeljic

and Blunt 2006) and water-air interface dynamics (Hassanizadeh et al. 2002). Detailed mapping

of pore spaces is often required for pore network modeling. Although recent method

developments such as x-ray tomography can help in determining pore networks, they still lack

the ability to differentiate between air and water phase (Van Loo et al. 2014). High

computational cost along with expensive and inaccurate pore network measurements methods

are the main challenges of this modeling approach.

1. Introduction

Using superposed groundwater-surface water equations to deal with both groundwater and

surface water flow is known as single-domain approach. Superposed equations are applied on

a transition zone between surface and groundwater and a gradual gradient (Figure 1-5 c) of key

flow parameters between them is defined. Using this approach, defining the thickness of the

transition zone as well as the parameters and their gradients are very challenging (Goyeau et al.

2003, Rosenzweig and Shavit 2007).

Figure 1-5: a) Sediment components, b) surface and groundwater flow at interface, c) single

domain and d) two domain groundwater-surface water coupling schemes; modified after

Mosthaf et al. (2011)

In two-domain approach, groundwater flow is modelled with groundwater equations and

surface water flow with surface water flow equations (Figure 1-5 d). The shared boundary

between two models which is the surface water-groundwater interface is then used to transfer

parameter information from one model to other (and back). Using this approach, there are

several coupling techniques who couple surface flow velocity to groundwater flow velocity

with adopting a slip coefficient. A relatively precise coupling outcome is expected through

coupling with feedback effects (two-way coupling, Nützmann and May 2007) where the

selected coupling parameter values from the surface water model run is fed to the groundwater

model. Resulting values from groundwater model run are returned to the surface water model

in an iteration step and surface water model is run again with the new value. Iteration steps

continue until desired parameter convergence is achieved. Flow conditions like velocity in

groundwater are orders of magnitude smaller than surface water which makes temporal

coupling interval selection of both models very challenging. The most common two-domain

coupling technique is one-way sequential coupling with pressure as the coupling parameter

(Tonina and Buffington et al. 2009, Trauth et al. 2015). This technique is widely used and

1. Introduction

validated for a variety of groundwater-surface water interaction cases including hyporheic zone

and bioturbation processes (Saenger et al. 2005, Cardenas and Wilson 2007a, Bardini et al.

2012, Brand et al. 2013, Trauth et al. 2013, Trauth et al. 2014, Chen et al. 2015 and 2018).

In previous sections, it was explained how for local scale processes such as hyporheic flow and

bioturbation, the surface water-groundwater interface itself was hosting crucial ecohydrological

processes. Via coupling, this interface is considered as a side boundary and its continous

transition is replaced by a discontinuity . Next to other advanced numerical models such as

DUMUX (Weishaupt et al. 2019) who applied the direct numerical solution of flow through a

random pore-network, recently, a high-resolution integral modeling approach (Oxtoby et al.

2013) was applied to hyporheic flow and transport processes through homogenous sediment

across rippled streambed by Broecker et al. (2019, 2021). This approach uses the same

conceptual approach for both groundwater and surface water without adding any additional

coupling parameters for the interface. To further develop this approach and its area of

applicability for local groundwater-surface water exchange modeling, it is necessary to validate

it in comparison to widely used coupled approaches and experiments. Local-scale high

resolution modeling via this approach is done through computational fluid dynamic model

OpenFOAM which is introduced in the following.

1.6 Flow and transport modeling with OpenFOAM

OpenFOAM (Open Field Operation And Manipulation) is the most widely used free, open-

source computational fluid dynamics software with a vast user base across many areas of

science both commercially and academically. It includes a wide range of features (libraries,

solvers, utilities) for complex fluid mechanics modeling as well as fields like acoustics, solid

mechanics and electromagnetics. Numerical fluid mechanics problems can be modeled via

OpenFOAM in compressible/incompressible, single-/multi-phase, 1-/2-/3-dimensions forms. It

is mostly based on the FVM (Finite Volume Method) in space and the FDM (Finite Difference

Method) in time (Schulze and Thorenz 2014).

In pre-processing phase, structured numerical mesh is generated using the built-in utility

‘blockMesh’. Additionally, unstructured meshes from other meshing tools like ‘GMSH’

(Geuzaine and Remacle 2009) and Salome (Ribes and Caremoli 2007) can be imported.

Boundaries as well as initial conditions should be selected accordantly to each other for model

stability and good convergence.

The partial differential equations governing flow and transport in OpenFOAM cannot be solved

directly and thus are spatially and temporally discretized as mentioned above and are transferred

1. Introduction

to algebraic equations. Based on computational capacity, desired outcome precision and

necessary model convergence, spatio-temporal discretization, margin of error and interpolation

of parameters can be defined by user; with considering limiting factors like Courant number

(Schulze and Thorenz 2014).

To keep computational capacity limits of the models presented in the current study small,

supercomputer cluster HPC (High Performance Computing) of the Technische Universität

Berlin is used for modeling. Mathematical equations, algorithms and modeling techniques of

OpenFOAM used here for modeling surface water flow as well as integral modeling of flow

and exchange between groundwater and surface water are presented in the following.

1.6.1 Modeling surface water flow and transport

Most surface water flow scenarios of the present study are modeled with InterFoam solver, a

component of OpenFOAM software. This widely used solver (Higuera et al. 2014, Schmitt et

al. 2015, Bayon et al. 2016) allows modeling of three-dimensional Navier-Stokes equations for

two (here water and air, multiphase in general) incompressible, immiscible fluids. In the current

study, it is often essential to consider water level fluctuations and their effects on pressure

distribution in the model as also was the case in Broecker et al. (2019). To capture the sharp

water-air interfaces, the VOF (Volume of Fluid) interface locating method implemented in

interFoam is applied. Fractions of air and water in partially saturated cells are defined with

altering fluid properties for each fraction. Volume of fractions are detected based on their

fraction indicator values (0 [air] < 𝛼𝛼< 1 [water]). Fluid properties include density and viscosity

which must be weighted according to their fractions (𝛼𝛼 [-]) and are calculated as follows:

𝜇𝜇=𝜇𝜇𝑤𝑤𝛼𝛼+𝜇𝜇𝑎𝑎(1 −𝛼𝛼) Eq. (1.1)

𝜌𝜌=𝜌𝜌𝑤𝑤𝛼𝛼+𝜌𝜌𝑎𝑎(1 −𝛼𝛼) Eq. (1.2)

where 𝜌𝜌 [kg/m3] is the density of the fluids (𝜌𝜌𝑤𝑤 = water density, 𝜌𝜌𝑎𝑎= air density) 𝑎𝑎𝑎𝑎𝑎𝑎 𝜇𝜇

[kg/(m s)] is the dynamic viscosity (𝜇𝜇= 𝜇𝜇𝑝𝑝 + 𝜇𝜇𝑡𝑡; physical + turbulent dynamic viscosity, 𝜇𝜇𝑤𝑤=

water dynamic viscosity, 𝜇𝜇𝑎𝑎= air dynamic viscosity).

The interface convection is described by the following equation:

𝜕𝜕α

𝜕𝜕t+𝛻𝛻∙(α𝑣𝑣)= 0 Eq. (1.3)

where 𝑣𝑣 [m/s] is the flow velocity vector and 𝑡𝑡 [s] is time.

Conservation of mass (Eq. (1.4)) and momentum (Eq. (1.5)) are written as follows:

1. Introduction

𝛻𝛻𝑣𝑣 = 0 Eq. (1.4)

�𝜕𝜕(𝜌𝜌𝑣𝑣

�



)

𝜕𝜕𝑡𝑡 +𝑣𝑣𝛻𝛻(𝜌𝜌𝑣𝑣)�=−𝛻𝛻𝛻𝛻+𝜇𝜇𝛻𝛻2𝑣𝑣+𝜌𝜌𝜌𝜌 Eq. (1.5)

where p [Pa] is the pressure and g [m/s2] is the gravity vector acceleration.

The above equations calculate velocity in three dimensions and pressure, while pressure-

velocity coupling in interFoam is done with PIMPLE algorithm which is based on the PISO

(Pressure Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for

Pressure Linked Equations) algorithm (Caretto et al. 1973, Issa 1986, Greenshields 2010).

To describe the transport of a conservative tracer with a concentration C [kg/m³], Broecker et

al. (2018) have implemented the following advection-diffusion equation into interFoam:

𝜕𝜕𝐶𝐶

𝜕𝜕𝑡𝑡+𝛻𝛻 ∙(𝑣𝑣𝑣𝑣)+𝛻𝛻∙(𝐷𝐷𝛻𝛻𝑣𝑣)= 0 Eq. (1.6)

D [m2/s] is diffusion coefficient which is divided into Dph [m²/s] for the physical diffusivity and

Dt [m²/s] for the turbulent diffusivity. Turbulent Schmidt number relates the turbulent

diffusivity to the turbulent viscosity as follows:

𝑆𝑆𝑆𝑆𝑡𝑡=𝜇𝜇𝑡𝑡/𝜌𝜌

𝐷𝐷𝑡𝑡 Eq. (1.7)

The physical diffusivity is a fluid property, while the turbulent Schmidt number Sct [-] can be

used for model calibration.

1.6.2 Turbulence models

In fluid dynamics, a fluid can either flow in parallel layers with no disruption between those

layers (laminar flow) or undergo chaotic changes in pressure and flow velocity both in time and

space (turbulent flow). The extent of turbulence in a flow is often characterized using Reynolds

number (Re [-]) which describes the ratio of inertia to viscous forces acting on the fluid.

Reynolds number is used to predict flow patterns in different flow conditions. In fully formed

open channel and pipe flow with Re < 2300 flow is considered laminar and with Re > 2300 it

is considered turbulent. While some modeling cases of this study have relatively low Reynolds

numbers and are considered laminar (section 3), others need to be modeled with a proper

turbulence model for precise high-resolution modeling outcomes (section 2 and 4).

OpenFOAM offers a wide range of turbulence models. The main three types of numerical

turbulence simulation methods are RANS, LES and DNS.

1. Introduction

A simpler turbulence model compared to others in terms of computational effort is RANS

(Reynolds-averaged Navier-Stokes). It consists of a variety of models such as k-ε, k-ω and k-ω

SST. k-ε is the most common RANS turbulence model to simulate mean flow characteristics

for turbulent flow conditions. It is a two-equation model that consists of two transport equations

to account for turbulent kinetic energy (k) and turbulent dissipation rate (ε). It is applicable for

free-shear flows and small pressure gradients but not good for complex flows with strong

curvatures. k-ω (ω is the frequency of energy dissipating eddies) is a RANS model with strength

close to walls. It is a two-equation model for turbulent energy (k) and frequency (ω) and is used

mostly for low Reynolds numbers. k-ω SST (k-ω Shear Stress Transport) is a combination of

the formerly mentioned RANS models and uses k-ε for inside of the model domain and k-ω for

the boundaries. Although computationally suitable, RANS models sometimes cannot properly

capture all flow velocity fluctuations and anisotropy of eddies

Large turbulent eddies are directly resolved with LES (Large Eddies Simulation) on relatively

fine grids. Applying a user-defined low-pass filtering, the smallest length scale is put aside and

its effects are modelled via subgrid algebraic models. In comparison to RANS, LES is much

more computationally burdening. In DNS (Direct Numerical Solutions) all spatial and temporal

eddy scales are resolved leading to extremely fine grids and extremely high computational

effort. Compared to DNS, LES is computationally reasonable. In this thesis next to laminar

modeling, for turbulence modeling LES is used.

1.6.3 Integral modeling of groundwater-surface water interactions

A modified version of interFoam to account for flow in porous medium was presented by

Oxtoby et al. (2013). For incompressible fluids the equations averaged over region ([ ]f) for the

conservation of mass (Eq. (1.4)) and momentum (Eq. (1.5)) are written as:

φ𝛻𝛻[𝑣𝑣]𝑓𝑓= 0 Eq. (1.8)

φ�𝜕𝜕[𝜌𝜌]𝑓𝑓[𝑣𝑣

�



]𝑓𝑓

𝜕𝜕𝑡𝑡 +𝛻𝛻[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓[𝑣𝑣]𝑓𝑓�=−φ𝛻𝛻[𝛻𝛻]𝑓𝑓+φ[𝜇𝜇]𝑓𝑓𝛻𝛻2[𝑣𝑣]𝑓𝑓+φ[𝜌𝜌]𝑓𝑓∙g + 𝐷𝐷 Eq. (1.9)

where φ [-] is the effective porosity and 𝐷𝐷 [kg/(m2s2)] is an additional porous drag term.

𝐷𝐷 = 𝐴𝐴 + 𝐵𝐵 Eq. (1.10)

𝐴𝐴 = −�150 (1−𝜑𝜑)

𝑑𝑑𝑝𝑝𝜑𝜑[𝜇𝜇]𝑓𝑓+ 1.75[𝜌𝜌]𝑓𝑓�[𝑣𝑣]𝑓𝑓��(1−𝜑𝜑)

𝑑𝑑𝑝𝑝[𝑣𝑣]𝑓𝑓 Eq. (1.11)

𝐵𝐵 = −0.34 (1−𝜑𝜑)

𝜑𝜑𝜕𝜕[𝑣𝑣]𝑓𝑓[𝜌𝜌]𝑓𝑓

𝜕𝜕𝑡𝑡 Eq. (1.12)

1. Introduction

where 𝑎𝑎𝑝𝑝 [m] is the effective grain size diameter.

In Eq. (1.11), the term “𝐴𝐴” describes pressure loss due to friction of the fluid with porous

medium after Ergun (1952). In porous medium, compared to free flow, more momentum is

needed to accelerate a given volume of water which is called added mass. Added mass is

accounted for in Eq. (1.12) (van Gent 1995).

In areas where only free flow exists, effective porosity is set to 1 which results in 𝐴𝐴=𝐵𝐵=𝐷𝐷=

0 and ensures the use of the original Navier-Stokes equations for free water flow (Eq. (1.4) and

Eq. (1.5)).

Former discussions about turbulence models apply here as well. Similar to surface water

transport modeling Broecker et al. (2021) have presented “porousInterTracer” for integral

modeling of conservative transport through porous zones:

𝜕𝜕[𝐶𝐶]𝑓𝑓

𝜕𝜕𝑡𝑡 +𝛻𝛻 ∙��𝑣𝑣𝑝𝑝

�



�𝑓𝑓[𝑣𝑣]𝑓𝑓�+𝛻𝛻∙([𝐷𝐷]𝑓𝑓𝛻𝛻[𝑣𝑣]𝑓𝑓)= 0 Eq. (1.13)

Where 𝑣𝑣𝑝𝑝

�



[m/s] is pore velocity.

In dealing with transport in porous medium, next to advection-diffusion based equations

provided in “porousInterTracer”, one could consider dispersion as well. Although for the cases

of this study, only diffusion was accounted for, verification of advective-diffusive-dispersive

version of porousInterTracer through comparison with analytical solutions is presented in the

appendix A.

1.7 Coupled modeling of groundwater-surface water flow

As mentioned in section 1.5 coupled models are the prevalent tools to investigate groundwater-

surface water interactions and are widely applied and validated. A necessary step to discuss the

advantages of the integral approach and address its area of application is to compare it to

coupled approach. Due to its high validity, wide application and prominence, one-way

sequential coupling with pressure as the coupling parameter is the best choice for this

comparison. In this thesis, first an identical coupled model is made for comparison to integral

model for the case of a rippled riverbed to inspect the full procedure of pre-processing,

simulating and post-processing of both models and second, integral model is compared to an

existing small-scale coupled model for the case of bioturbation. For surface water part of the

coupled model, interFoam solver of the OpenFOAM is used (discussed in 1.6.1). Thereafter, to

discuss the necessity of choosing a computationally expensive CFD surface water model for

1. Introduction

coupling, it is compared to a simpler shallow water equations modelShallow water equations

are applicable, when the wavelength is at least 20 times larger than the water depth leading to

hydrostatic pressure. These equations are derived from depth-integration of the Navier-Stokes

equations and in consequence they do not include vertical components. 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:

𝜕𝜕ℎ

𝜕𝜕𝑡𝑡 + 𝛻𝛻(𝑣𝑣ℎ)= 0 Eq. (1.14)

𝜕𝜕(𝑣𝑣

�



ℎ)

𝜕𝜕𝑡𝑡 +𝛻𝛻�𝑣𝑣·𝑣𝑣ℎ+1

2 𝜌𝜌ℎ2�− 𝜇𝜇𝑡𝑡𝛻𝛻2(𝑣𝑣ℎ) = −𝜌𝜌ℎ𝛻𝛻𝑧𝑧𝐵𝐵+𝑛𝑛2𝜌𝜌𝜌𝜌

ℎ1

3𝑣𝑣∙|𝑣𝑣| Eq. (1.15)

where ℎ[m] is the water depth, 𝑧𝑧𝐵𝐵[m] is the bottom elevation, 𝜇𝜇𝑡𝑡 [m2/s] is turbulent kinematic

viscosity and 𝑎𝑎 [s/m1/3] is Manning’s roughness coefficient.

Java-based modelling framework hms (hydroinformatics modeling system, Simons et al. 2014)

which uses the Finite Volume Method is chosen for shallow water flow modeling.

Coupling is achieved by transferring the modeled pressures results of the shared boundary from

surface water (interFoam) to groundwater model. These are then used as boundary conditions

of the PCSiWaPro® groundwater model. PCSiWaPro® is a 2D Finite Element-based

groundwater modeling software developed by IBGW Leipzig (Ingenieurbüro für Grundwasser

GmbH Leipzig) and TU Dresden. It uses the the extended Darcy’s law (Richards equations)

which is introduced in the mass conservation equation as follows:

𝜕𝜕𝜕𝜕

𝜕𝜕𝑡𝑡 =𝛻𝛻[𝐾𝐾(𝜃𝜃)(𝛻𝛻𝜓𝜓−1)] Eq. (1.16)

where 𝜃𝜃 [m3/m3] is the volumetric water content, 𝐾𝐾(𝜃𝜃) [m/s] is the hydraulic conductivity, 𝜓𝜓

[m] is the pressure head. Further information can be found on PCSiWaPro® user guide

(PCSiWaPro® Benutzerhandbuch).

1.8 Scope of this thesis

1.8.1 Current research deficits

Primarily, Broecker et al. (2019) have applied the integral approach to model flow and exchange

across a rippled streambed domain. Nevertheless, a systematic comparison of this approach to

prevalent coupling approach which is famously applied in groundwater-surface water

interaction investigations is missing. Extensive model comparison is necessary to prove the

1. Introduction

flexibility and applicability of the integral approach and its superiority to the coupled approach

in modeling small-scale high-resolution groundwater-surface water exchange processes.

In investigating cases such as flow and ventilation in and around invertebrate-dwelled burrows,

current experimental approaches only observe the overlying water flow processes and lack the

ability to measure flow and exchange inside the burrow. Existing coupled models that

investigate flow inside and around the burrow cannot show the connected flow system of the

overlying water, burrow and the surrounding sediment. Integral approach however is potentially

capable of modeling connected flow and exchange over, through and around the burrow.

Validation of the integral approach for tracking the tracer concentration along and across the

groundwater-surface water interface, is the key for its application to water quality assessment

and groundwater contamination modeling. Broecker et al. (2021) have validated the approach

for tracer transport in homogenous sediment. However, sediment heterogeneity which strongly

impacts contamination propagation by changing solute residence times in the sediment has not

yet been validated for the integral approach.

1.8.2 Aim of this thesis

An integral modeling approach to describe high-resolution flow and transport processes

occurring at the interface of groundwater and surface water is investigated in this thesis for

local (small) scale groundwater-surface water interactions. A conclusive comparison of the

approach next to the coupled approach is done for two distinctive cases and advantages of using

integral approach to model local (small) scale flow near surface water-groundwater interface

are discussed and area of application is determined. Thereafter, to extend the area of application

to modeling tracer transport in heterogenous sediments, flow patterns and tracer propagation at

heterogenous rippled streambeds with ambient groundwater flow are investigated and validated

with the help of experiments.

The thesis is structured in six sections consisting of the current introduction, three peer-

reviewed journal articles (one accepted, two submitted), further supplementary contributions

and a synthesis. Additionally, further solver validation and complementary information are

provided in two appendices.

Section 2 describes flow over and through a porous and rippled streambed. This is done through

modeling flow and exchange processes with integral and coupled approaches. Integral approach

uses the CFD software OpenFOAM to solve modified Navier-Stokes equations. A Richards

equation-based groundwater model (PCSiWaPro®) is coupled with an interFOAM surface

1. Introduction

water model and the necessity of using interFOAM as the surface water model is investigated

against a shallow water equations-based model (hms) for the same setup. Due to integral

approach being computationally slightly more expensive than coupled approach, area of

application meaning distance from the groundwater-surface water interface where coupled and

integral models lead to different results is analysed.

In section 3 flow induced by bioturbation activity of the macroinvertebrate Chironomus

plumosus worms in a U-shaped burrow in lake sediments is modelled via integral approach and

compared to similar setups modeled with coupled approach and experiments. Importance of

modeling continuous surface water-groundwater flow exchanges in these burrows as hotspots

of bioirrigation is emphasized in comparison to coupled approach which neglects feedback

effects and to experiments which lack representing flow through sediment and are mostly

conducted only for tracking flow over burrows.

By comparing the integral to coupled approach in section 2 and 3, advantages of the integral

approach are revealed and its superiority to model flow and exchange in the groundwater-

surface interaction space is assessed for a variety of cases.

In section 4 flow and tracer transport processes for heterogenous sediment were studied under

neutral, losing and gaining groundwater flow conditions. For this purpose, the extended integral

solver for modeling tracer propagation in surface water and sediment is used and further

validated by performing numerical modeling of flow and tracer propagation similar to

experiments of flow and tracer transport over a heterogenous rippled streambed.

An overview of the supplementary scientific work is presented in section 5. This section

includes a brief introduction of two co-authored journal articles concerning flow and exchange

across rippled streambed under the influence of ripple geometries, surface hydraulics and grain

sizes as well as flow and transport through homogenous surface water-groundwater interface.

In section 6 the main findings of this study are conjunctively summarized and their common

goal of applying, validating and extending the integral approach as a suitable modeling

instrument for groundwater-surface water interface is discussed; limitations of the approach are

addressed and further steps are described.

In appendix A, flow and transport integral model is validated for modeling advective-diffusive-

dispersive tracer propagation with the help of analytical solutions. In appendix B, an overview

of the modeled cases is given.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

2. Comparison of integral approach to coupled approach for

flow across rippled streambed

This study was published in international Journal on Geomathematics as:

________________________________________________________________________

Sobhi Gollo, V., Broecker, T., Marx, C., Lewandowski, J., Nützmann, G. & Hinkelmann R.

2022. A comparative study of integral and coupled approaches for modeling hydraulic

exchange processes across a rippled streambed, International Journal on Geomathematics 13,

no 16: 1-27. https://doi.org/10.1007/s13137-022-00206-5.

This article is licensed under a Creative Commons Attribution 4.0 International License, which

permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source, provide a link to the

Creative Commons license, and indicate if changes were made. The images or other third-party

material in this article are included in the article’s Creative Commons license, unless indicated

otherwise in a credit line to the material. If material is not included in the article’s Creative

Commons license 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 license, visit http://creativecommons.org/licenses/ by/4.0/.

___________________________________________________________________________

This is the postprint version of the article.

The test cases setups are listed in Appendix B (B1 Seepage through a dike modelled with

porousInter and PCSiWaPro®, B2 Seepage through a rectangular dam modelled with

porousInter and PCSiWaPro®, B3 one-dimensional, one-phase modeling of surface water flow

over rippled streambed with hms (shallow water equations), B4 Three-dimensional, two-phase

modeling of surface water flow over rippled streambed with InterFoam, B5 Two-dimensional,

modeling of groundwater flow in the sediment under the rippled streambed with PCSiWaPro®,

B6 Three-dimensional, two-phase integral modeling of flow across rippled streambed with

porousInter.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

2.1 Abstract

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” (CM),

which result from the successive execution of a surface water model and a groundwater model,

to a recently developed “integral model” (IM). 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.

2.2 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). 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). In recent years,

2. Comparison of integral approach to coupled approach for flow across rippled streambed

interest in understanding the importance of surface water and groundwater interactions and their

integrity as a single continuum has significantly increased (Bobba 2012). 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).

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 & Bencala, 1993; Tonina and Buffington, 2009, Trauth et al.

2014), and ripples (Cardenas & Wilson 2007a, 2007b, Bottacin‐Busolin & 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 & 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 &

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

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

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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 local-scale 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. 2013, Trauth et al. 2014, Trauth et al. 2015,

Chen 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 advection-

dominant. 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. 2020), 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 Section 2.3.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 porousInter extension includes porosity and grain diameter as parameters over

the region that characterize a porous medium in the model. Due to the high computational

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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 and Fox 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

(Figure 2-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 Richards equations; see Guo et al. 2017;

PC Sickerwasserprognose (seepage flow prognosis); see Section 2.3.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 Section 2.3.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 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 below.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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?

2.3 Materials and methods

2.3.1 Surface water modeling

In this study, we first modeled the surface water flow over of the rippled streambed (SW-CM;

Figure 2-1) using the interFoam solver of the OpenFOAM software (based on the Navier-Stokes

equations). We also investigated the possibility of model simplification by modeling the same

domain using shallow water equations.

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. 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).

Cardenas and Wilson (2007a), Tonina and Buffington (2009), and Janssen et al. (2012)

recommended a Navier-Stokes equations solver when running simulations of relatively

2. Comparison of integral approach to coupled approach for flow across rippled streambed

complex geometries such as those in this study, because flow and pressure distributions are

more realistically captured compared to simpler approaches.

Eq. (2-1) depicts the indicator fraction (0 [air] < 𝛼𝛼< 1 [water]) for the interface convection

equation:

𝜕𝜕α

𝜕𝜕t+𝛻𝛻∙(α𝑣𝑣)= 0 Eq. (2-1)

where 𝑣𝑣 [m/s] is the flow velocity vector and 𝑡𝑡 [s] is time.

Conservation of mass (Eq. (2-2)) and momentum (Eq. (2-3)) are thus written as follows:

𝛻𝛻𝑣𝑣 = 0 Eq. (2-2)

�𝜕𝜕(𝜌𝜌𝑣𝑣

�



)

𝜕𝜕𝑡𝑡 +𝑣𝑣𝛻𝛻(𝜌𝜌𝑣𝑣)�=−𝛻𝛻𝛻𝛻+𝜇𝜇𝛻𝛻2𝑣𝑣+𝜌𝜌𝜌𝜌 Eq. (2-3)

where 𝜌𝜌 [kg/m3] is the density of the fluids, p [Pa] is the pressure, 𝜇𝜇 [m2/s] is the dynamic

viscosity (𝜇𝜇= 𝜇𝜇𝑝𝑝ℎ + 𝜇𝜇𝑡𝑡; physical + turbulent viscosity), and g [m/s2] 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

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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:

𝜕𝜕ℎ

𝜕𝜕𝑡𝑡 + 𝛻𝛻(𝑣𝑣ℎ)= 0

𝜕𝜕(𝑣𝑣ℎ)

𝜕𝜕𝑡𝑡 +𝛻𝛻�𝑣𝑣·𝑣𝑣ℎ+1

2 𝜌𝜌ℎ2�− 𝜇𝜇𝑡𝑡𝛻𝛻2(𝑣𝑣ℎ) = −𝜌𝜌ℎ𝛻𝛻𝑧𝑧𝐵𝐵+𝑎𝑎2𝜌𝜌𝜌𝜌

ℎ1

3𝑣𝑣∙|𝑣𝑣|

where ℎ[m] is the water depth, 𝑧𝑧𝐵𝐵[m] is the bottom elevation, 𝜇𝜇𝑡𝑡 [m2/s] is turbulent kinematic

viscosity, and 𝑎𝑎 [s/m1/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.

2.3.2 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. (2-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:

𝜕𝜕𝜕𝜕

𝜕𝜕𝑡𝑡 =𝛻𝛻[𝐾𝐾(𝜃𝜃)(𝛻𝛻𝜓𝜓−1)] Eq. (2-5)

where 𝜃𝜃 [m3/m3] is the volumetric water content, 𝐾𝐾(𝜃𝜃) [m/s] is the hydraulic conductivity, and

𝜓𝜓 [m] is the pressure head.

Eq. (2-4)

2. Comparison of integral approach to coupled approach for flow across rippled streambed

2.3.3 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-2)) and momentum (Eq. (2-3)) are written as modified Navier-Stokes equations;

[ ]f indicates an averaged parameter over a void region:

φ𝛻𝛻[𝑣𝑣]𝑓𝑓= 0 Eq. (2-6)

φ�𝜕𝜕[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓

𝜕𝜕𝑡𝑡 +𝛻𝛻[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓[𝑣𝑣]𝑓𝑓�=−φ𝛻𝛻[𝛻𝛻]𝑓𝑓+φ[𝜇𝜇]𝑓𝑓𝛻𝛻2[𝑣𝑣]𝑓𝑓+φ[𝜌𝜌]𝑓𝑓∙g + 𝐷𝐷 Eq. (2-7)

where φ [-] is the effective porosity and 𝐷𝐷 [kg/(m2s2)] is an additional porous drag term.

𝐷𝐷 = 𝐴𝐴 + 𝐵𝐵 Eq. (2-8)

𝐴𝐴 = −�150 (1−𝜑𝜑)

𝑑𝑑𝑝𝑝𝜑𝜑[𝜇𝜇]𝑓𝑓+ 1.75[𝜌𝜌]𝑓𝑓�[𝑣𝑣]𝑓𝑓��(1−𝜑𝜑)

𝑑𝑑𝑝𝑝[𝑣𝑣]𝑓𝑓 Eq. (2-9)

𝐵𝐵 = −0.34 (1−𝜑𝜑)

𝜑𝜑𝜕𝜕[𝑣𝑣]𝑓𝑓[𝜌𝜌]𝑓𝑓

𝜕𝜕𝑡𝑡 Eq. (2-10)

where 𝑎𝑎𝑝𝑝 [m] is the effective grain size diameter.

The term “𝐴𝐴” in Eq. (2-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 𝐴𝐴=𝐵𝐵=𝐷𝐷= 0 and allows the use of the original Navier-

Stokes equations for free water flow (Eq. (2-2) and Eq. (2-3)). Similarly, for the SW-CM, for

the IM, the LES turbulence model was used.

2.3.4 Comparable hydrogeological conditions assessment

The sediment parameters φ and 𝑎𝑎𝑝𝑝 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, 𝜃𝜃𝑠𝑠 [m3/m3] (saturated volumetric water content), 𝜃𝜃𝑟𝑟 [m3/m3] (the residual

volumetric water content), and 𝐾𝐾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 (𝜃𝜃𝑠𝑠 or φ) equal to the porosity chosen in porousInter.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

When the system is fully saturated, residual volumetric water content is not relevant. For

porousInter, a particle size (𝑎𝑎𝑝𝑝) of 2 mm and a porosity (φ) of 0.25 was chosen. The saturated

hydraulic conductivity (𝐾𝐾0) used in PCSiWaPro® was derived from Eq. (2-11) (following

Hazen) using the same particle size chosen for porousInter:

𝐾𝐾0= 0.00116 𝑎𝑎10

2 Eq. (2-11)

where 𝑎𝑎10 is a 10% fall-through of the sieve curve [mm].

Since porousInter is limited to a single uniform particle size, 𝑎𝑎10 = 2 mm and, with Eq. (2-11),

𝐾𝐾0 = 4.64 × 10-3 m/s was used in this study.

2.3.5 Computational domain

Figure 2-1: Model setup and boundary conditions of a) the shallow water equations surface

water model; b) the coupled model (CM), including the surface water component of the coupled

model (SW-CM, b, top) and the groundwater component of the coupled model (GW-CM, b,

bottom); c) the integral model (IM); and d) the ripple geometry.

In Figure 2-1, light blue indicates the air phase, dark blue the water phase, and light brown the

sediment, including respective boundary conditions. Figure 2-1a shows the model setup of the

surface water flow using hms (shallow water equations model) for single-phase flow over a

2. Comparison of integral approach to coupled approach for flow across rippled streambed

rippled streambed. We used this model to investigate if model simplification was applicable to

the SW-CM (see Figure 2-1b, top). The CM is shown in Figure 2-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 Figure 2-1c is identical to that of Figure 2-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 Figures 2-1c and 2-1d.

Model setups of the SW-CM and the IM were created in three dimensions and had a 1m 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 near-wall 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.

2.3.6 Meshes, initial and boundary conditions and parameters

This study utilizes unstructured meshes, which were generated using the GMSH three-

dimensional finite element mesh generator (for the IM and the SW-CM). The SW-CM (Figure

2-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 (Figure 2-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

2. Comparison of integral approach to coupled approach for flow across rippled streambed

= 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³/s was set as the boundary condition at the inlet on the left

side in Figure 2-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 60s and 300s 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 a ~0.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/m1/3, which is the coefficient recommended for clean and

straight natural streams.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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 Figure 2-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. The next step was to consider the effects of simplifying the models by

replacing the SW-CM with a shallow water model.

2.4 Results

2.4.1 IM and CM flow velocity fields in the model domain

To investigate the flow velocities in the model domain, we chose a representative slice (8th

ripple in Figure 2-1; 735 cm < 𝑥𝑥 < 760 cm, -50 cm < 𝑧𝑧 < 50 cm in Figure 2-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.

Figure 2-2: Graphic visualization of the ripples at the sediment–water interface and the

horizontal (blue dotted lines) and the vertical (green dotted lines) cross sections used to compare

flow results of the integral model (IM) and the coupled model (CM).

To compare the results, we exhibited the simulated flow velocities across several vertical and

horizontal cross sections. Vertical cross sections intersect with the luv (𝑥𝑥 = 748 cm) and the lee

(𝑥𝑥 = 758 cm) sides of the 8th ripple and extend through the sediment and the water parts (Figure

2-2). We placed the horizontal cross sections in the water (𝑧𝑧 = 40 cm), the sediment (𝑧𝑧 = -40

2. Comparison of integral approach to coupled approach for flow across rippled streambed

cm), and the interface-adjacent sediment (𝑧𝑧 = -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 2-3: A comparison of flow velocities 𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧 in the x (a, c) and z (b, d) directions for

vertical cross sections at x = 748 cm (a, b) and x = 758 cm (c, d) for the integral model (IM)

and the coupled model (CM); absolute values of subtracting flow velocities have been

calculated using the CM and the IM (|𝑣𝑣𝐶𝐶𝐶𝐶 − 𝑣𝑣𝐼𝐼𝐶𝐶|) in the x (e) and z (f) directions for both

vertical cross sections.

-0.23 0.00 0.23 0.46 0.69 0.92 1.15

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

x = 748 cm

z (m)

(m/s)

vIM

vCM

x = 748 cm

x = 758 cm x = 758 cm

vIM

vCM

-0.22 -0.11 0.00 0.11 0.22 0.33 0.44 0.55 0.66 0.77

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

z (m)

(m/s)

-0.23 0.00 0.23 0.46 0.69 0.92 1.15

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

z (m)

(m/s)

- v

||v

- v

-0.22 -0.11 0.00 0.11 0.22 0.33 0.44 0.55 0.66 0.77

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

z (m)

(m/s)

vIM

vCM

vIM

vCM

-0.23 0.00 0.23 0.46 0.69 0.92 1.15

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

z (m)

(m/s)

x = 748 cm

x = 758 cm

e x = 748 cm

x = 758 cm

-0.22 -0.11 0.00 0.11 0.22 0.33 0.44 0.55 0.66 0.77

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

z (m)

(m/s)

2. Comparison of integral approach to coupled approach for flow across rippled streambed

Figure 2-3 illustrates the flow velocities of both models (𝑣𝑣𝐶𝐶𝐶𝐶 and 𝑣𝑣𝐼𝐼𝐶𝐶) through vertical cross

sections at 𝑥𝑥 = 748 cm (Figures 2-3a and 2-3b) and 𝑥𝑥 = 758 cm (Figures 2-3c and 2-3d). These

are shown separately for the velocities in the x direction (𝑣𝑣𝑥𝑥; Figures 2-3a and 2-3c) and the z

direction (𝑣𝑣𝑧𝑧; Figures 2-3b and 2-3d). CM velocities were subtracted from IM velocities for

each vertical cross section and in both directions. The absolute values of these subtractions are

displayed in Figures 2-3e and 2-3f. Major differences between the flow velocities of the two

models in both directions (𝑥𝑥 and 𝑧𝑧) were detected in the area close to the interface (around 𝑧𝑧 =

0 cm). A few decimeters away from the interface (both in groundwater and surface water),

however, these differences fade and reach values close to zero.

Figure 2-4: A comparison of flow velocities 𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧 in the x (a) and z (b) directions for the

horizontal cross section at z = -2 cm along 735 cm < x < 755 cm for the integral model (IM)

and the coupled model (CM); absolute values of subtracting flow velocities have been

calculated using the CM and the IM (|vCM − vIM|) in the x (c) and z (d) directions.

7.35 7.40 7.45 7.50 7.55

-2.0×10

-3

-1.5×10

-3

-1.0×10

-3

-5.0×10

-4

0.0

5.0×10

-4

1.0×10

-3

1.5×10

-3

2.0×10

-3

Vx (m/s)

x (m)

VIM

VCM

7.35 7.40 7.45 7.50 7.55

-2.0×10

-3

-1.5×10

-3

-1.0×10

-3

-5.0×10

-4

0.0

5.0×10

-4

1.0×10

-3

1.5×10

-3

2.0×10

-3

Vz (m/s)

x (m)

VIM

VCM

7.35 7.40 7.45 7.50 7.55

0.0

5.0×10

-4

1.0×10

-3

1.5×10

-3

2.0×10

-3

Vx (m/s)

x (m)

|VCM- VIM|

7.35 7.40 7.45 7.50 7.55

0.0

5.0×10

-4

1.0×10

-3

1.5×10

-3

2.0×10

-3

Vz (m/s)

x (m)

|VCM- VIM|

2. Comparison of integral approach to coupled approach for flow across rippled streambed

To further substantiate the fact that unlike areas near the interface, differences between the flow

velocities farther from the interface are negligible in the two models, we present the flow

velocities across horizontal cross sections near the interface (z = -2 cm, Figure 2-4) and away

from it (z = -40 cm, Figure 2-5; z = 40 cm, Figure 2-6).

Figure 2-5: A comparison of flow velocities 𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧 in x (a) and z (b) directions for the

horizontal cross section at z = -40 cm along 735 cm < x < 755 cm for the integral model (IM)

and the coupled model (CM); absolute values of subtracting flow velocities have been

calculated using the CM and the IM (|vCM − vIM|) in the x (c) and z (d) directions.

Figures 2-4, 2-5, and 2-6 show the absolute value of the differences between the simulated

velocities of the CM and the IM in both directions (𝑥𝑥 and 𝑧𝑧) as well. Looking at Figure 2-4 (z

= -2 cm), flow velocity differences between the two models reach as high as 2 × 10-3 m/s. For

deeper groundwater (Figure 2-5, z = - 40 cm), these differences are almost zero in the z direction

(𝑣𝑣𝑧𝑧) and peak at 2.5 × 10-5 in the x direction (𝑣𝑣𝑥𝑥).

7.35 7.40 7.45 7.50 7.55

1×10-5

2×10-5

3×10-5

4×10-5

5×10-5

6×10-5

7×10-5

8×10-5

9×10-5

1×10-4

Vx (m/s)

x (m)

VIM

VCM

7.35 7.40 7.45 7.50 7.55

1×10-5

2×10-5

3×10-5

4×10-5

5×10-5

6×10-5

7×10-5

8×10-5

9×10-5

1×10-4

Vz (m/s)

x (m)

VIM

VCM

7.35 7.40 7.45 7.50 7.55

0.00

5.20×10-6

1.04×10-5

1.56×10-5

2.08×10-5

2.60×10-5

Vx (m/s)

x (m)

|VCM - VIM|

7.35 7.40 7.45 7.50 7.55

0.00

5.20×10-6

1.04×10-5

1.56×10-5

2.08×10-5

2.60×10-5

Vz (m/s)

x (m)

|VCM - VIM|

2. Comparison of integral approach to coupled approach for flow across rippled streambed

Figure 2-6: A comparison of flow velocities 𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧 in the x (a) and z (b) directions for the

horizontal cross section at z = 40 cm along 735 cm < x < 755 cm for the integral model (IM)

and the coupled model (CM); absolute values of subtracting flow velocities have been

calculated using the CM and the IM (vCM − vIM|) in the x (c) and z (d) directions.

As displayed in Figure 2-6, maximum flow velocity differences between the IM and the CM

(for both 𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧) 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 2-1) by the averaged flow velocities of the IM and the CM (|𝑣𝑣𝐶𝐶𝐶𝐶

�

| and

7.35 7.40 7.45 7.50 7.55

1.062

1.064

1.066

1.068

1.070

1.072

1.074

1.076

1.078

1.080

(m/s)

x (m)

7.35 7.40 7.45 7.50 7.55

6.0×10

-3

8.0×10

-3

1.0×10

-2

1.2×10

-2

1.4×10

-2

1.6×10

-2

1.8×10

-2

(m/s)

x (m)

7.35 7.40 7.45 7.50 7.55

0.0

1.0×10

-3

2.0×10

-3

3.0×10

-3

4.0×10

-3

5.0×10

-3

6.0×10

-3

7.0×10

-3

8.0×10

-3

(m/s)

x (m)

- v

7.35 7.40 7.45 7.50 7.55

0.0

1.0×10

-3

2.0×10

-3

3.0×10

-3

4.0×10

-3

5.0×10

-3

6.0×10

-3

7.0×10

-3

8.0×10

-3

(m/s)

x (m)

- v

2. Comparison of integral approach to coupled approach for flow across rippled streambed

|𝑣𝑣𝐼𝐼𝐶𝐶

�

|, Table 2-1). The results (𝐷𝐷𝐶𝐶𝐶𝐶,𝐷𝐷𝐼𝐼𝐶𝐶) are shown for all the horizontal cross sections and in

both flow directions (𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧) in Table 2-1.

Table 2-1: |𝑣𝑣𝐶𝐶𝐶𝐶 − 𝑣𝑣𝐼𝐼𝐶𝐶

�

|, |𝑣𝑣𝐼𝐼𝐶𝐶

�

| and |𝑣𝑣𝐶𝐶𝐶𝐶

�

| velocities and |𝑣𝑣𝐶𝐶𝐶𝐶 − 𝑣𝑣𝐼𝐼𝐶𝐶

�

|velocities over |𝑣𝑣𝐼𝐼𝐶𝐶

�

| and

|𝑣𝑣𝐶𝐶𝐶𝐶

�

| velocities for z = -2 cm, -40 cm and 40 cm in the 𝑥𝑥 and 𝑧𝑧 directions.

𝐷𝐷𝐶𝐶𝐶𝐶 and 𝐷𝐷𝐼𝐼𝐶𝐶 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 𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧, 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. 𝐷𝐷𝐶𝐶𝐶𝐶 and 𝐷𝐷𝐼𝐼𝐶𝐶 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.

𝑣𝑣

𝑥𝑥

(m/s)

𝑣𝑣

𝑧𝑧

(m/s)

z = -2

z = -40 cm

z = 40

z = -2 cm

z = -40

Z = 40

|𝑣𝑣

𝐶𝐶𝐶𝐶

− 𝑣𝑣

𝐼𝐼𝐶𝐶

�

|(m/s)

4.9×10-4

1.95×10-5

2.60×10-3

4.59×10-4

2.65×10-7

3.33×10-3

|𝑣𝑣

𝐶𝐶𝐶𝐶

�

| (m/s)

5.75×10-5

6.02×10-5

1.071

8.17×10-5

1.41×10-5

1.35×10-2

|𝑣𝑣

𝐼𝐼𝐶𝐶

�

| (m/s)

2.56×10-5

7.97×10-5

1.073

6.83×10-5

1.39×10-5

1.02×10-2

𝐷𝐷𝐶𝐶𝐶𝐶 =

|𝑣𝑣

𝐶𝐶𝐶𝐶

− 𝑣𝑣

𝐼𝐼𝐶𝐶

�

|𝑣𝑣𝐶𝐶𝐶𝐶

�

| (-)

8.60

3.20×10-1

2.43×10-3

5.61

1.87×10-2

2.46×10-1

𝐷𝐷𝐼𝐼𝐶𝐶 =

|𝑣𝑣

𝐶𝐶𝐶𝐶

− 𝑣𝑣

𝐼𝐼𝐶𝐶

�

|𝑣𝑣𝐼𝐼𝐶𝐶

�

| (-)

19.28

2.40×10-1

2.42×10-3

6.71

1.90×10-2

3.27×10-1

2. Comparison of integral approach to coupled approach for flow across rippled streambed

2.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, Figure 2-7) of the two models on the 8th ripple (see Figure 2-1 above), however, differ

from each other on the ripple surface. These differences are displayed in Figure 2-7b.

Figure 2-7: a) Hydrodynamic pressure values at the surface water–groundwater interface on the

investigated ripple (representative cross section, 8th ripple) for the IM and the CM; b)

subtraction of the hydrodynamic pressures of the CM from that of the IM on the investigated

ripple.

Another indicator of the differences between the two models near the interface is shown in

Figure 2-8. Higher flow velocities and bigger turbulent eddies can be seen within the troughs

between the ripples of the CM (Figure 2-8b) compared to those of the IM (Figure 2-8a). In

Figure 2-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 (Figure 2-9).

7.40 7.45 7.50 7.55 7.60

4900

5000

5100

5200

5300

5400

p_rgh (Pa)

x (m)

7.40 7.45 7.50 7.55 7.60

100

150

200

250

300

350

400

450

p_rgh (Pa)

x (m)

CM - IM

2. Comparison of integral approach to coupled approach for flow across rippled streambed

Figure 2-8: Flow velocities and directions over and through the investigated ripple (8th ripple,

Figure 1) for a) the integral model and b) the coupled model. Note: red dashed line indicates

the sediment–water interface; black arrows only indicate directions, and are not scaled

regarding magnitude. Magnitude is shown as a color scale.

Figure 2-9: Velocity and Reynolds number (in log10 scale) distributions through the investigated

ripple (8th ripple, Figure 2-1) for a) the integral model (IM) and b) the coupled model (CM).

The patterns and hotspots of the flow differ when comparing the interface-adjacent flow in the

groundwater (Figure 2-9). We also realized that a zone on the crest of the ripple in the IM (the

2. Comparison of integral approach to coupled approach for flow across rippled streambed

area highlighted green in Figure 2-9) has a velocity of 0.02 m/s < 𝑣𝑣𝐼𝐼𝐶𝐶 < 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 Section 2.5.

2.4.3 Shallow water simplification

Figure 2-10: A comparison of the horizontal surface water flow (𝑣𝑣𝑥𝑥) for vertical cross sections

at a) x = 748 cm and b) x = 758 cm for the surface water component of the coupled model (SW-

CM) and the shallow water equations in the surface water model (SW-SM).

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. Using shallow water equations,

𝑣𝑣𝑧𝑧 is zero and 𝑣𝑣𝑥𝑥 values along a vertical cross section are constant. Nevertheless, we generated

two parabolic velocity profiles (according to Bartam and Balance 1996) for vertical cross

sections at x = 748 cm and x = 758 cm based on the average velocities of the shallow water

equations model (𝑣𝑣𝑥𝑥

�

748𝑐𝑐𝑐𝑐 = 0.803 𝑐𝑐

𝑠𝑠,𝑣𝑣𝑥𝑥

�

758𝑐𝑐𝑐𝑐 = 0.807 𝑐𝑐

𝑠𝑠). As shown in Figure 2-10,

generated 𝑣𝑣𝑥𝑥 flow profiles of the shallow water equations model (𝑣𝑣𝑆𝑆𝑆𝑆−𝑆𝑆𝐶𝐶) are very different

compared to the SW-CM (which uses the Navier-Stokes equations). Furthermore, the 𝑣𝑣𝑧𝑧 = 0

assumption in shallow water equations contradicts the 𝑣𝑣𝑧𝑧 profiles across the vertical cross

sections displayed in Figure 2-3. In shallow water equations, pressure is assumed to be

-0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2

0.0

0.1

0.2

0.3

0.4

0.5

z (m)

vx (m/s)

SW-CM

SW-SM

-0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2

0.0

0.1

0.2

0.3

0.4

0.5

x = 758 cm

z (m)

vx (m/s)

x = 748 cm

vSW-CM

vSW-SM

2. Comparison of integral approach to coupled approach for flow across rippled streambed

hydrostatic, which conflicts with the hydrodynamic pressure results displayed in Figure 2-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.

2.4.4 Required computational resources

Modeling the SW-CM via hms requires only 30 minutes 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 high-performance

computing cluster of Technische Universität Berlin, it takes about 24 hours for the SW-CM to

complete the simulation. Using the same resources, it takes about 70 hours for the IM to run

completely. For the GW-CM, about one hour of computing with the aforementioned PC is

sufficient.

2.5 Discussion

In the introduction, 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.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

2.5.1 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. 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

(Sobhi Gollo et al. 2021). Nevertheless, for flow across a morphologically modified (rippled)

streambed, the current study is a further step in validating the integral approach. In Figure 2-3,

we can see that when starting from the interface (𝑧𝑧 ~ 0), flow velocities of the models in both

flow directions (𝑣𝑣𝑥𝑥 and 𝑣𝑣𝑧𝑧) 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.

2.5.2 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 Figure 2-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

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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 Figure 2-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 Figure

2-9 is therefore different. The formation of one large eddy in the trough between ripples (see

Figure 2-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.,

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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.

2.5.3 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 (𝑣𝑣𝑧𝑧) 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

(𝑣𝑣𝑧𝑧), eddies modify horizontal velocities (𝑣𝑣𝑥𝑥) 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 Figure 2-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.

2.5.4 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 hours 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.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

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.

2.6 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 water–groundwater 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 water–groundwater 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.

2. Comparison of integral approach to coupled approach for flow across rippled streambed

2.7 Appendix

porousInter verification for seepage flow

We chose the seepage of water through a simplified rectangular dam and a dike as test cases;

the results of the integral solver (porousInter), PCSiWaPro®, analytical solutions (Kobus and

Keim 2001 (1D), Di Nucci 2015 (2D), Casagrande 1937 [extended Kozeny] (2D), Lattermann

2010 [Kozeny] (2D)), and other numerical solutions (Aitchison and Coulson 1972, Westbrook

1985) are compared in Figure 2-11. A particle size (𝑎𝑎𝑝𝑝) of 2 mm (𝐾𝐾0 = 4.64 × 10-3 m/s) and

porosity of 0.25 were used. The results in Figure 2-11 demonstrate very good agreement

between most analytical solutions and model results, including PCSiWaPro® and the integral

solver (porousInter). Only the 1D analytical solution deviates from the rest, since 1D tended to

oversimplify the solution here.

Figure 2-11: Model geometry and results of different models and analytical solutions for

seepage through a rectangular dam (left) and a dike (right).

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

3. Flow simulations in and around a ventilated U-shaped

chironomid dwelled burrow using integral approach

This study was published in Journal of Ecohydraulics as:

___________________________________________________________________________

Sobhi Gollo, V., T. Broecker, J. Lewandowski, G. Nützmann, and R. Hinkelmann. 2021. An

integral approach to simulate three-dimensional flow in and around a ventilated U-shaped

chironomid dwelled burrow, Journal of Ecohydraulics,

https://doi.org/10.1080/24705357.2021.1938258.

This is an Accepted Manuscript of an article published by Taylor & Francis in Journal of

Ecohydraulics on 02.07.2021, available online:

http://www.tandfonline.com/10.1080/24705357.2021.1938258

& Francis Group

___________________________________________________________________________

This is the postprint version of the article.

Flow simulations using porousInter are listed in Appendix B7 and B8 for the test cases 1 and 2

respectively.

3.1 Abstract

Tube dwelling of chironomids often dominate benthic communities in freshwater ecosystems

with high population density and pumping rates which strongly enhance exchange across the

sediment-water interface and impact biogeochemical processes. Such processes are

investigated by tracking the flow initiated by chironomid’s pumping through and around

burrows using laboratory and computer models. We used modeling and experimental results of

other authors considering U-shaped burrows embedded in the sediment to improve process

understanding and prove the plausibility of an integral modeling approach. In contrast to

coupled models of pipe (burrow), free surface (overlying water column) and groundwater flow

(surrounding sediment), we present a novel high-resolution integral formulation for the porous

medium-free water domain (called porousInter as part of OpenFOAM (Open Field Operation

and Manipulation)) solving extended versions of the Navier-Stokes equations which allows us

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

to simultaneously simulate flow in the burrow, the overlying water column and the surrounding

sediment to better account for feedback effects between the sediment and free water. Using

similar model setup as of a coupled approach, we performed scenarios of flow through burrow

and sediment triggered by pumping in the center of the burrow. Plausible agreement of our

integral model with results of a coupled model and with experimental results was obtained when

comparing flow patterns around the burrows, between two burrow branches and at burrow inlet

and outlet.

3.2 Introduction

Charles Darwin realized the relevance of biological reworking of soils and sediments (Darwin

1881, Meysman et al. 2006). Nowadays bioturbation is defined as “all transport processes

carried out by animals that directly or indirectly affect sediment matrices. These processes

include both particle reworking and burrow ventilation” (Kristensen et al. 2012). Bioturbators

are considered ecosystem engineers (Meysman et al. 2006, Kristensen et al. 2012). Some

bioturbators live in marine and freshwater sediment where they build open-ended or blind-

ended burrows and actively pump overlying water through their burrows (ventilation) to supply

themselves with oxygen and/or food. Burrow ventilation causes advective and diffusive

exchange of solutes between overlying water, water in the burrows and sediment pore water

called bioirrigation (Kristensen et al. 2012).

The pumping activities (ventilation) of macroinvertebrates such as Chironomus plumosus have

severe impacts on hydrodynamic and biogeochemical processes in lake sediments as well as on

exchange processes between pore water and the overlying water body (Hupfer et al. 2019;

Hoelker et al. 2015, Brand et al. 2013, Roskosch et al. 2010 and 2011 and 2012, Morad et al.

2010, Lewandowski et al. 2007). Applying the 2D pore water sampler introduced by

Lewandowski et al. (2002), Lewandowski et al. (2005a, 2005b, 2005c) used laboratory and

field studies that systematically investigated the influence of macrozoobenthos (including C.

plumosus) on the spatial distribution of pore water phosphorus concentrations. The results of

the laboratory experiments showed that there is a change in retention of phosphorus

(Lewandowski et al. 2005a, 2005b), ammonium and ferrous iron (Laskov et al. 2007,

Lewandowski et al. 2007) in the sediment in the presence of macrozoobenthos. Field

experiments supported the hypothesis that macrozoobenthos causes small-scale horizontal

pore-water heterogeneity and regulates phosphorus exchange (Lewandowski et al. 2005c). In

fact, for chironomid-dwelled burrows, there is a positive diffusive iron (II) and ammonium flux

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

from the pore water to the burrows and a reduced phosphorus release across the sediment-water

interface (Lewandowski et al. 2007).

Chironomus plumosus larvae are known to build single U-tubes down to a sediment depth of

10 to 20 cm (Hölker et al. 2015). Roskosch et al. (2010) reported that C. plumosus can be found

in amounts up to 10,000 individuals 1/m2 in lake sediments. In Lake Müggelsee (Berlin,

Germany) Roskosch et al. (2011) measured a density of 745 4th instar larvae 1/m2 and reported

a population pumping rate of 1.09 to 1.45 m3 m-2 d-1 depending on the assumed burrow diameter

(1.6 or 2.0 mm). This value shows the impressive impact of chironomids on hydrodynamics

and solute exchange in lake sediments which impacts the above-mentioned biogeochemical

processes (Roskosch et al. 2010).

Although to date, there have not been any in situ flow measurements inside burrows dwelled

by Chironomus plumosus larvae due to technical difficulties associated with the measurement

in fragile, 1.6 mm diameter burrows in soft, muddy sediment (Morad et al. 2010), there are

several measurements from mesocosm laboratory experiments of Roskosch et al. (2010) and

Morad et al. (2010) as well as modeling approaches (e.g. Brand et al. 2013). Brand et al. (2013)

used during larval pumping activity a flow velocity of 0.009 m/s in chironomid burrows in their

coupled approach to model the hydrodynamic exchanges between sediment and burrow. That

value is rather small comparted to Morad et al. (2010) who measured a mean flow velocity of

0.0137 to 0.0154 m/s during pumping periods, Roskosch et al. (2010) determined 0.0150 m/s

and Roskosch et al. (2011) reported 0.0149 m/s. There is an alternation of pumping and resting

phases every few minutes. The flow velocity of pumping and resting phases can be averaged

using 33 min pumping time per hour which was determined by Roskosch et al. (2011) resulting

in an average flow velocity of 0.0082 m/s. They also realized that pressures initiated by water

fluxes depend on burrow shape, sediment properties and ventilation activities of the

macrozoobenthos. Brand et al. (2013) used an average (pumping + resting phases) flow velocity

of 0.04 m/s.

Modeling of pumping activities of macrozoobenthos in U-shaped burrows started with the

assumption of radially symmetric tubes. Aller (1980) set up a marine tube irrigation model

assuming that the water in the burrow has the same composition as the overlying water and that

a diffusive transfer mechanism transports solutes in the surrounding sediment (bioirrigation).

(Chemical) mass transfer, pore water transfer, impacts of sediment type and depth of the burrow

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

have been discussed in later modeling studies (e.g. Furukawa et al. 2001, Koretsky et al. 2002,

Meysman et al. 2006, Meile et al. 2005, Stief et al. 2010a).

Brand et al. (2013) presented a 3D coupled model (free fluid flow in burrow and porous medium

flow in sediment) of the pumping activity of a Chironomus plumosus larvae around a U-shaped

burrow. They simulated the hydrodynamics of the flow in the porous medium and applied it to

the transport of an inert tracer as well as of oxygen using different sediment permeabilities. The

model of Brand et al. (2013) gave insights into flow fields and flow directions as well as

pressure distributions in the sediment around the burrow which could hardly be assessed using

laboratory experiments. However, we think that this model like many other coupling schemes,

lacks crucial aspects. In Brand et al. (2013) the overlying water above the sediment had not

been considered. Their approach can be considered as one-way coupling without feedback. This

limitation is due to the inability of the coupled approach to simultaneously model flow in the

sediment as well as in the burrow and in the overlying water.

In this contribution we present an alternative to coupled models, an integral model for the

investigation of a U-shaped burrow under the influence of pumping activities. This approach

includes the groundwater-surface water interface as a crucial exchange zone directly in the

model concept. It was initially developed by Oxtoby et al. (2013) who extended the InterFoam

solver of OpenFOAM, which is an Open-Source solver for the Navier-Stokes equations, by a

porosity and a special resistance term in the porous domain (Broecker et al. 2019). This

extension called porousInter is an alternative to prevalent coupled approaches such as the one

from Brand et al. (2013). PorousInter has never been applied to small scale pore water-surface

water exchange processes such as invertebrate-dwelled burrows in lake sediments where there

is the interaction of free flow in the burrow with the surrounding sediment and the overlying

water column. Two separate models for sediment (using groundwater flow model) and free

water (using e.g. Navier-Stokes equations) are coupled through the sediment-water interface as

a common closed boundary. The integral approach includes the interface as well as both pore

water and surface water in a single model using a modified (porosity included) version of

Navier-Stokes equations resulting in a model which is very much likely to capture continuous

flow exchange on the interface between surface water and groundwater in cases where it is

crucial. In the present work, the same setup as of Brand et al. (2013) is used for two reasons:

On the one hand as a validation case and on the other hand to discuss the differences between

the two approaches. We aim to present more realistic flow and pressure fields of U-shaped

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

burrows under the influence of the pumping activity of Chironomus plumosus larvae and pupae.

In addition, we include the overlying water above the sediment for plausibility tests according

to the experiments of Roskosch et al. (2010) and Morad et al. (2010). A wider range of flow

velocities and sediment characteristics compared to Brand et al. (2013) are employed in this

study in order to satisfy the discussions on how a connected system of sediment and surface

water will react to pressure/flow fields on the interface compared to a constant pressure

boundary coupling scheme and to relate to results of physical experiments.

3.3 Materials and methods

3.3.1 Computational domain

Figure 3-1: Model domain (lengths are not to scale) and boundaries B1-B4 reported in Table 1;

B1 is atmospheric pressure at cylinder top; B2 are cylinder walls and cylinder bottom; B3 is

symmetry plane to reduce computational effort; B4 represents the pumping chironomid larvae.

The left panel is a detail enlargement of the burrow shown in the right panel

Figure 3.1 illustrates the geometry of the computational domain. This geometry has been

selected to match the one of Brand et al. (2013). The only altered lengths for some scenarios

are the radius and the height of the sediment cylinder which are both 0.25 m in our study while

they were both 1 m in Brand et al. (2013). The impacts of this reduction have been checked by

running some test scenarios and comparing the results in the areas close to the burrow where

no considerable changes have been observed. Therefore, the geometry in Figure 3.1 acts nearly

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

identical to the one of Brand et al. (2013). This step was necessary to reduce the considerable

computational effort. Using the same number of processors, the down sized model requires only

10 % of the CPU time when compared to the large model.

The water table marked yellow in Figure 3.1 shows 0.1 m overlying column of water on top of

the saturated sediment. In comparison to no overlying water shown with a green line this is an

alternative scenario later used for comparison to other modeling approaches and experimental

results. The burrow has a tube-like shape with a diameter of 0.002 m. It goes 0.150 m into the

sediment starting from the inlet, including a curved part in the lower area and is symmetrical.

The pumping zone of the Chironomus plumosus is located at the bottom of the burrow (B4,

marked with a short vertical red line).

3.3.2 Governing equations

The Open-Source computational fluid dynamics software OpenFOAM (Open Field Operation

and Manipulation) version 2.4.0 is chosen to simulate flow processes in and around the U-

shaped burrow using the interFoam solver which is applicable for multiphase fluid flows. It is

based on the Navier-Stokes equations and is widely used in hydraulic engineering in recent

years (Schulze and Thorenz 2014). In detail we have applied and extended porousInter (Oxtoby

et al. 2013), a sub-branch of interFoam. Our solver has been developed for two-phase (water,

air) flow in the context of groundwater-surface water interactions (Broecker et al. 2019). In the

present contribution, only single phase (water) flow is investigated. In simulations of relatively

complex geometries as in our case, flow and pressure distributions are realistically captured

using Navier-Stokes equations (Cardenas and Wilson 2007a, Tonina and Buffington 2009,

Janssen et al. 2012). For incompressible fluids and accounting for the porous medium, the

equations for the conservation of mass (Eq. (3.1)) and momentum (Eq. (3.2)) are written as ([ ]f

is averaged parameter over void region):

φ𝛻𝛻[𝑣𝑣]𝑓𝑓= 0 Eq. (3.1)

φ�𝜕𝜕[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓

𝜕𝜕𝑡𝑡 +𝛻𝛻[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓[𝑣𝑣]𝑓𝑓�=−φ𝛻𝛻[𝛻𝛻]𝑓𝑓+φ[𝜇𝜇]𝑓𝑓𝛻𝛻2[𝑣𝑣]𝑓𝑓+φ[𝜌𝜌]𝑓𝑓∙g + 𝐷𝐷 Eq. (3.2)

where φ [-] is the effective porosity, 𝑣𝑣 [m/s] is the flow velocity, 𝛻𝛻 [Pa] is pressure, 𝜌𝜌 [kg/m3]

is the density of the fluids, 𝑡𝑡 [s] is time, 𝜇𝜇 [kg/(m s)] is the dynamic viscosity (𝜇𝜇= 𝜇𝜇𝑝𝑝ℎ𝑦𝑦𝑠𝑠 +

𝜇𝜇𝑡𝑡𝑡𝑡𝑟𝑟𝑡𝑡; physical + turbulent viscosity), g [m/s2] is the gravitational acceleration and 𝐷𝐷

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

[kg/(m2s2)] 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 Eq. (3.3):

𝐷𝐷 = 𝐴𝐴 + 𝐵𝐵 Eq. (3.3)

𝐴𝐴 = −�150 (1−φ)

dpφ[μ]f+ 1.75[ρ]f�[𝑣𝑣]f��(1−φ)

dp[𝑣𝑣]f Eq. (3.4)

𝐵𝐵 = −0.34 (1−φ)

φ𝜕𝜕[𝑣𝑣]f[ρ]f

𝜕𝜕t Eq. (3.5)

where 𝑎𝑎𝑝𝑝 [m] is the effective grain size diameter.

The term “𝐴𝐴” in Eq. (3.4) considers pressure loss due to friction of the fluid with porous medium

after Ergun (1952). In sediment, compared to free flow, more momentum is needed to accelerate

a given volume of water. This is called added mass (van Gent 1995) because the extra

momentum needed suggests that a larger volume of water has to be accelerated. “𝐵𝐵” in Eq. (3.5)

acts to add this added mass of the fluid due to flow recirculation caused by porous medium after

van Gent (1995). In areas where only free flow exists, effective porosity is set to 1 which results

in 𝐴𝐴=𝐵𝐵=𝐷𝐷= 0 and ensures the use of the original Navier-Stokes equations for free water

flow.

The modeling tool uses the Finite-Volume-Method in space and the Finite-Difference-Method

in time and allows massively parallel computations. In order to couple pressure and velocity, a

combination of the SIMPLE (Semi-Implicit Method for Pressure Linked Equations, outer-

correction loop) and the PISO (Pressure-Implicit with Splitting of Operators, inner corrector

loop) algorithms called PIMPLE (merged PISO and SIMPLE) for robust modeling is

implemented (Rodrigues et al. 2011).

According to the definition of the Reynolds number (Eq. (3.6)) (Reynolds 1883):

𝑅𝑅𝑅𝑅= 𝜌𝜌𝑣𝑣𝜌𝜌

𝜇𝜇 Eq. (3.6)

where 𝐿𝐿 [m] is a characteristic length.

and considering the highest velocity in the burrow in this study (0.097 m/s), the density of water

(1000 kg/m3), the dynamic viscosity of the water (~0.001 kg/(m s)), the burrow diameter 𝐿𝐿

(0.002 m), the critical Reynolds number 2300 for the transition between laminar and turbulent

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

flow, the highest Reynolds number here is ~200 which means that flow is always laminar and

therefore, turbulence is not considered in the model.

3.3.3 Grid and boundary conditions

To set up a grid which is comparable with both, the one of Brand et al. (2013) and the

experiments of Roskosch et al. (2010) and Morad et al. (2010) for the domain shown in Figure

3.1, we generated a model with 0.5 million cells which are distributed with higher resolution

near and inside the burrow (cell length < 0.5 mm) and coarser resolution (cell length ~ 1 cm) in

the other regions of the soil using SALOME 9.3.0. With 0.5 million cells grid convergence was

obtained, i.e. the results do not change when the mesh is further refined. For this study we

assumed a rigid porous medium. Movement of sediment including erosion and deposition is not

considered. As initial condition, the system was fully-filled with water. We set boundaries as

shown in Table 3.1. We defined the top boundary with atmospheric pressure (B1, Table 3.1).

The top boundary is the sediment surface (Figure 3.1, green line) or the water surface of the 10

cm overlying water column (Figure 3.1, yellow line).

Table 3-1: Boundary locations and conditions for geometry of Figure 1 and parameter values

Boundary Location Pressure (Pa) Velocity (m/s)

B1 cylinder top atmospheric pressure pressure defined velocity

(pressureInletOutletVelo

city)

B2 cylinder sides and bottom velocity

(fixedFluxPressure)

B3 planar surface symmetry plane symmetry plane

circular patch at bottom

of burrow

velocity

(fixedFluxPressure)

0.009 -0.023

Effective grain size diameter [mm] 0.125 -0.250

Effective porosity [-

]

0.39

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

Case 1 is set up with comparable boundary conditions as the coupling scheme of Brand et al.

(2013) while Case 2 is designed to be able to track the movement of water into the burrow inlet

and out of the burrow outlet. Thus, there is the opportunity of a comparison of the model with

the experiments of Roskosch et al. (2010) and Morad et al. (2010).

The system is a cylindrical column with closed sides. Therefore we define no flow conditions

by setting velocity to zero on the side walls and on the bottom (B2, Table 3.1). In order to

reduce computational efforts, only half of the cylinder is considered (Figure 3.1). Consequently,

a planar surface (B3, Table 3.1) acts as a symmetry plane to mirror the other half of the cylinder.

We then defined a baffle in the bottom of the burrow (B4, Table 3.1) which is a feature in

OpenFOAM allowing the user to turn internal faces into boundary faces. The use of this baffle

is to mimick the pumping activity of Chironomus plumosus larvae (or pupae) in the bottom of

the burrow. Like any other boundaries in OpenFOAM, velocity or/and pressure conditions can

be applied to such a baffle. The defined baffle has in fact a half-circular surface perpendicular

to the tubed shaped burrow. Flow is free to enter and exit a baffle from both sides. Equal to the

maximum velocity value used in Brand et al. (2013), we have initially chosen 0.009 m/s as the

internal pumping velocity in the baffle. Thereafter, we investigated the impact of higher

pumping velocities on pressure and flow recirculation patterns in the sediment. It is applied in

0.15 m depth below the sediment surface in the center of the burrow. The difference to the

boundaries of Brand et al. (2013) is that we need no inner pressure boundary between burrow

and sediment as it is required and common in coupled approaches. For achieving steady state,

we run each instationary model for 40 to 100 s which took about 2-3 days CPU time on 40

processors of a high performance computer (TU Berlin HPC).

3.3.4 System parameters

Previous investigations including Brand et al. (2013), Morad et al. (2010), Hupfer et al. (2019),

and Roskosch et al. (2010, 2011, 2012) have all used grain size distributions found in Lake

Müggelsee (Berlin, Germany) to provide insight into pumping activities of Chironomus

plumosus. Brand et al. (2013) have set the range of permeabilities from 3 × 10-12 to 1 × 10-14 m2

to cover the range of fine sandy to unconsolidated loamy fine sandy sediments.

PorousInter includes sediment in the model by using effective porosity and an additional drag

term (see Eq. 3.3-3.5) which depends among other terms on the effective grain size diameter

(here uniform grain size distribution) (Oxtoby et al. 2013). The relation between the intrinsic

permeability 𝐾𝐾 [m²] and hydraulic conductivity 𝐾𝐾𝑓𝑓 [m/s] is given by Eq. (3.7):

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

𝐾𝐾𝑓𝑓=𝐾𝐾 𝜌𝜌𝜌𝜌

𝜇𝜇 Eq. (3.7)

The relation can be simplified to (Hinkelmann 2005):

𝐾𝐾𝑓𝑓=𝐾𝐾 ∙107

𝑐𝑐.𝑠𝑠 Eq. (3.8)

if 𝐾𝐾𝑓𝑓 is given in [m/s] and 𝐾𝐾 in [m2].

Therefore, the range of hydraulic conductivities used by Brand et al. (2013) is 3 × 10-5 to 1 ×

10-7 m/s. We used the following equation to convert hydraulic conductivity to effective grain

size (after Hazen):

𝐾𝐾𝑓𝑓= 0.00116 𝑎𝑎10

2 Eq. (3.9)

where 𝑎𝑎10 is 10 % fall through of the sieve curve [mm].

Using Eq. (3.8) and Eq. (3.9), 𝑎𝑎10 ranges for the permeabilities used by Brand et al. (2013)

between 0.160 and 0.009 mm. Brand et al. (2013) have chosen 𝑎𝑎50 = 0.25 mm. Assuming a

single uniform grain size (𝑎𝑎10 = 𝑎𝑎50) and since the main comparison with Brand et al. (2013) is

based on cases with 𝐾𝐾 = 3 × 10-12 m2, the effective grain size is set to 0.25 mm. Scenarios with

a smaller grain size diameter down to 0.125 mm are also calculated. These values are in

accordance with experiments of Roskosch et al. (2010) and Morad et al. (2010).

According to DIN 4220:2008-11 (2008), the effective porosity is set to 0.39 corresponding to

the porosities of different sandy soils, from pure to middle loamy sand (similar material as used

by Brand et al. 2013, Roskosch et al. 2010 and Morad et al. 2010) (Table 3.1).

3.3.5 Modeling cases

Two different cases with small adjustments in model geometry are performed.

Case 1: In order to compare the integral approach to the coupled approach of Brand et al. (2013),

the geometry of Figure 3.1 is used and the green line is taken as the top boundary. The advective

flow was observed for uniform grain size diameter of 0.25 mm and the pumping velocity of

0.009 m/s (identical to Brand et al. 2013). The flow in and around inlet, outlet and pumping

zone (B4, Figure 3.1) as well as recirculation patterns through the sediment are analysed.

Furthermore, impacts of different pumping velocities and effective grain size diameters were

investigated. Case 1 allows the replication of the pumping activity of chironomid larvae in the

burrow systems using a standard case of a U-shaped burrow with focus on exchange processes

at the interface between pore water and surface water. Added to the comparison to the coupled

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

approach, this case gives insight towards the effect of pumping activity of Chironomus

plomosus on flow characteristics by mimicking the pumping activity of the larvae using

pumping velocities at B4 and a flow though a structure designed to replicate the natural upward

burrow shape.

Case 2: By adding an overlying water column of 10 cm on top of the saturated system (yellow

line as the top boundary in Figure 3.1, extended version: 1 m height and radius), the aim was to

investigate flow patterns over burrow inlet and outlet with the intention to compare these flow

patterns to experimental results of Morad et al. (2010) and Roskosch et al. (2010) and to derive

a relation between velocity at inlet/outlet to pumping velocity inside the burrow triggered by

Chironomus plumosus larvae.

3.4 Results and discussion

3.4.1 Comparison of integral to coupled approach

Internal pumping velocity and pressure differences at the pumping zone

Figure 3-2: Pressure (left) and velocity (right) fields of Case 2 (top) and Case 1 (bottom) with

pumping velocity = 0.009 m/s and 𝑎𝑎𝑝𝑝 (effective grain size) = 0.25 mm. Green and yellow lines:

see Figure 3.1

Using the geometry of Figure 3.1 and sediment characteristics of Table 3.1, the goal was to test

the integral approach for feasibility in ventilated U-shaped burrow case. By applying a 2400

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

N/m3 force, Brand et al. (2013) created a velocity equal to 0.009 m/s during pumping periods.

In Brand et al. (2013), this corresponded to pressure reduction from -6 Pa right behind the

pumping zone to 0 at the inlet, followed by a 6 Pa pressure right in front of the pumping zone

which then reduces to 0 at the outlet. In the coupled approach of Brand et al. (2013), the burrow

is assumed to be a pipe where the pressure in the burrow is imposed as boundary condition in

the groundwater flow model. No feedback of the surrounding sediment on flow in the burrow

is taken into account. With the integrated model this setup was modelled as Case 1 which is

shown as lower panel of Figure 3.2. The lower left panel of Figure 3.2 shows the pressure

distribution results for the same setup.

In Case 1, the resulting pressure difference was not only caused by the flow in the burrow

(unlike the assumption of one-way pressure head boundary from Brand et al. (2013)) but also

by the flow in the sediment. Therefore a range of -2.9 to +2.9 Pa was observed. In fact, with the

integral approach, -6 to +6 Pa were only achieved right and left of the baffle by applying 0.023

m/s as internal boundary velocity for the pumping zone. In general, symmetric pressure and

velocity fields can be seen. A velocity field similar to the one of the coupled approach is found

around the burrow (Figure 3.2, bottom right, Case 1). For Case 1, absolute pressures on the left

and right of the pumping zone are equal. However, they differ in sign.

Figure 3-3: Pressure in the pumping zone (orange = left branch of the burrow, blue = right

branch of the burrow) for different water depths

-50000

50000

100000

150000

200000

0 5 10 15 20

Total Pressure (Pa)

Water column (m)

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

Figure 3.2 top shows a similar setup for Case 2. Compared to Case 1, under the same velocity

conditions Case 2 experiences a larger total pressure (Ptotal = Ppumping + Patm + Pρgh) in the pumping

zone with Ppumping - pumping pressure, Patm - atmospheric pressure and Pρgh - hydrostatic pressure

which increases linearly with depth below the water surface. due to added pressure (Pρgh) of the

overlying water column (Figure 3.2). To investigate the impact of overlying water height on total

pressure at pumping zone, we ran case 1 with an altering overlying column of water as an

addition to the top boundary (top boundary B1 = atmospheric pressure + pressure caused by a

column of water). Figure 3.3 shows the relation between lake bed depth to pressure at the

pumping zone for a range of lake bed depths with U = 0.009 m/s and dp= 0.25 mm. A linear

trend is observed.

The pumping pressure distribution in the sediment depends on the effective grain size. With

growing effective grain size, water eases more to flow through sediment and thus the pressure

drops. With smaller grain sizes, the system will look more like pipe flow. Flow is more

concentrated in the burrow and less distributed in the sediment. The relation between pressure

P at the pumping zone and the ratio of the maximum velocity U at the pumping zone to the

effective grain size diameter 𝑎𝑎𝑝𝑝 is shown in Figure 3.4. Using 10 different cases with varying

velocities and grain sizes a linear approximation between pressure and U/𝑎𝑎𝑝𝑝 ratio is obtained.

It shows that higher pumping velocities (U) and smaller grain sizes (𝑎𝑎𝑝𝑝) result in a higher

absolute value of the pressure at the pumping zone.

Figure 3-4: Maximum/minimum pressure at the pumping zone for 10 different cases of U/𝑎𝑎𝑝𝑝;

U is the maximum velocity at the pumping zone and 𝑎𝑎𝑝𝑝 is the effective grain size diameter; blue

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

stands for the right side of the baffle where pressure is positive and orange stands for the left

side of the baffle where pressure is negative

Pressure and flow fields

The velocity in a porous medium is generally orders of magnitude smaller than in free flow

(surface water or pipe flow). In a burrow embedded in fine sediment, the flow is similar as in a

pipe where a fully developed semi-circular velocity profile occurs between the walls. To see if

OpenFOAM can handle smooth transitions of the maximum flow in the center of the pipe to its

‘walls’ where velocity is nearly zero, the flow in the burrow was investigated. Such a

plausibility test has not been published so far to our knowledge. Figure 3.5 shows the flow

profile inside the inlet and outlet branches of the burrow in a cross section at z = -0.05 m in

which a semi-circular flow profile can be seen. The formation of this semi-circular flow profile

provides a good basis for investigations of interactions between free and porewater flow. The

flow in the sediment is not only influenced by the effective porosity and grain size, but also by

the velocity and pressure field in the interacing burrow. In Brand et al. (2013) the mass

exchange between burrow and sediment is driven by a pressure boundary condition, while in

our approach such a forcing is not required as the pressure along the burrow wall, the mass

exchange and feedback effects between burrow and sediment are directly included in the model

concept. Therefore, we think our approach is more realistic when compared to Brand et al.

(2013).

Figure 3-5: Velocity and pressure profiles at z = -0.05 m through inlet and outlet branches of

the burrow for 𝑎𝑎𝑝𝑝 (effective grain size) = 0.25 mm, Case 1

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

Another noticeable aspect is the distribution of pressure between two branches of the burrow.

The pressure difference of these two branches is caused by the pumping activity at the bottom

of the burrow. Figure 3.5 shows the pressure distribution along the cross section at z = -0.05 m

which is point-symmetrical with zero pressure at x = 0 m. The highest and the lowest pressures

occur in the burrows. This induces (1) flow in the burrows and (2) flow through the sediment

from the burrow branch with high pressure (outlet) to the burrow branch with low pressure

(inlet). The pressure differences are zero at the sediment surface and gradually increase with

depth.

Figure 3.6 shows flow directions in the domain (Figure 3.6a), close to pumping zone (Figure

3.6b) and between two branches (Figure 3.6c) for pumping velocity = 0.009 m/s and 𝑎𝑎𝑝𝑝= 0.25

mm. Overall, flow directions and magnitudes show a very good agreement with the results of

Brand et al. (2013) with the exception of areas near the burrow where flow exchange is

considered differently in the integral approach.

Figure 3-6: Velocity magnitudes (colors) and flow directions (arrows) a) in the entire domain,

b) around pumping zone and c) between inlet and outlet branches

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

Figure 3.6 shows that there is a flow recirculation between inlet and outlet branch of the

chironomid burrow caused by the pressure distribution shown in Figure 3.4. This recirculation

is more intense close to the pumping zone (Figure 3.6b) and less intense further away from it.

Compared to Brand et al. (2013), lower pressures (e.g. ± 2.9 compared to ± 6 Pa in Brand et al.

(2013) at pumping zone) are observed in areas close to the burrow and velocities larger than

10-4 or 10-5 m/s, respectively, occur in smaller areas compared to results of Brand et al. (2013).

Figure 3-7: Dominant flow directions (white arrows) on z = -0.05 m plane, looking from a) top

b, c) lateral sides; blue plane is the cross section of a half cylinder

To make a comparison to the coupled approach of Brand et al. (2013), lateral velocities and

directions for z = -0.05 m are shown in Figure 3.7. We found out that dominant flow from outlet

to inlet branch (Figure 3.6b and 3.6c) is similar to the one shown in Brand et al. (2013).

However, a circular flow field occurs (Figure 3.7a and 3.7b). Water leaving the outlet burrow

branch is transported from the right side of the outlet branch to the left side of the inlet branch

and enters the burrow again. For 𝑎𝑎𝑝𝑝=0.25 mm and pumping velocity U = 0.01 m/s there is a

difference between the percentage of water leaving the burrow in different areas.

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

3.4.2 Comparison of flow patterns to experiment results

A further step towards better understanding the hydraulics of ventilated U-shaped burrows as

well as proof of plausibility of the integral approach is achieved by focusing on flow above

burrow inlet and burrow outlet. Using Case 2, the modelled results are compared to flow

patterns of experiments of Roskosch et al. (2010) and Morad et al. (2010). Velocity vector fields

shown by Morad et al. (2010) map flow magnitudes and directions. Higher magnitudes occur

close to inlet/outlet. Figure 3.8 shows a comparison of modeling results and results of Morad et

al. (2010) for flow directions and flow fields above inlet and outlet. Values obtained in the

modeling effort (≤ 38 mm/s over inlet/outlet) are among the range of Morad et al. (2010) (~35

mm/s over inlet/outlet). It can be seen that the model has successfully simulated the outward

(outlet) and inward (inlet) flows at thresholds.

Figure 3-8: Comparison of the modeled flow pattern (left) above inlet (top row) and outlet

(bottom row) with results of Morad et al. (2010) (right); dashed lines in modeled images mark

the sediment surface; note: values differ due to unidentical burrow lengths and pumping

velocities. Images based on the same measurements/raw data as images shown in Morad et al.

(2010)/Roskosch et al. (2010)

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

We tested the integral approach (Case 1, 𝑎𝑎p (effective grain size) = 0.25 mm) and related the

maximum inlet/outlet velocity (1.23 mm/s) to the triggering pumping velocity of 7.98 mm/s.

Using a wider range of pumping velocities (5 to 25 mm/s, Figure 3.9), we tried to find a

relationship where we could predict pumping velocity for given inlet/outlet velocities for this

geometry. Figure 3.9 displays the relation between flow at inlet/outlet to pumping velocity.

Throughout all scenarios, pumping velocity is ~6.5 times bigger than the inlet/outlet velocity.

This ratio depends on the burrow geometry, i.e. its depth. It can be generally said that having a

shallower burrow will also reduce the difference between the velocity at inlet/outlet and the

velocity at the pumping zone for identical pumping velocities.

Figure 3-9: Velocity at pumping zone to velocity at inlet/outlet (measured at outlet) for 5

different scenarios; 𝑎𝑎𝑝𝑝 (effective grain size) = 0.25 mm

An important point that Roskosch et al. (2010) have addressed is the accumulation of advective

flow as the result of pumping activity (mentioned as bioirrigation in their paper, also shown in

Morad et al. 2010). With the integral approach, this was observed in a setup which allows good

mass exchange between burrow and sediment. There is a slight difference between flow

velocity at inlet and outlet imminent overlying water reported by Roskosch et al. (2010) and

illustrated by Morad et al. (2010). Similar patterns where velocity at outlet is slightly (<1%)

higher than the one at inlet were seen. Velocity patterns over inlet and outlet branches agree

with observations of Morad et al. (2010).

0.5

1.5

2.5

3.5

4.5

0 5 10 15 20 25 30

Flow at inlet/outlet (mm/s)

Flow at pumping zone (mm/s)

3. Flow simulations in and around a ventilated U-shaped chironomid dwelled burrow using integral approach

3.5 Conclusion

This study aimed at gaining a better understanding of flow and exchange processes in, above

and around U-shaped chironomid-dwelled burrows and to present an innovative method in

dealing with such ecohydraulic topics. For that purpose the integral solver of Oxtoby et al.

(2013) was applied and compared to experiments and other modeling approaches. In the integral

approach mass exchange and feedback effects between burrow and sediment as well as

sediment and overlying water column are directly included in the model concept, while in

coupled approaches they are forced by ‘internal’ boundary conditions and therefore cannot

account for feedback effects.

Plausible results were obtained for advective flow in the burrow (Roskosch et al. 2010).

Sensivity analysis on pressure, effective grain size and pumping velocity resulted in a linear

relation between pressure at the pumping zone and the ratio of pumping velocity to effective

grain size. Using the identical setup as a prevalent coupled approach (Brand et al. 2013), flow

pattern and magnitudes qualitatively agreed. Further, the formation of the semi-circular flow

between the burrow branches was investigated.

Considering the overlying water, the missing link between the experiments of Roskosch et al.

(2010) and the modeling of Brand et al. (2013) was filled. Roskosch et al. (2010) focused on

areas above inlet and outlet and Brand et al. (2013) on flow in the sediment surrounding the

burrows, while we could do both investigations together with the integral approach. The

pumping activities of the Chironomus plumosus larvae control directly the flow in the burrows

and in the surround sediment as investigated in the present study. This alters not only the local

biogeochemistry in the sediment by oxygen supply (Lewandowski et al. 2007), particle

reworking (Kristensen et al. 2012) and long-term fixation of phosphate (Hupfer et al. 2019) but

also the water body above the sediment and thus, has severe impacts on the entire lake

ecosystem including the food web (Hölker et al. 2015). Our approach seems to be interesting

in quantifying the biogeochemical relevance of burrows for surface water quality. On the other

hand one has to take into account that the computational effort is huge. The drag term of the

integral solver has been derived for coarse and sandy soils. its application to very fine material

such as clay is hardly recommended. A more realistic comparison of the formation of advective

flow using this approach to bioturbation of animals which dwell in sandy sediment like Nereis

diversicolor or Arenicola marina (Kristensen et al. 2012) is expected.

4. Flow and transport modeling in heterogenous sediments using integral approach

4. Flow and transport modeling in heterogenous sediments using

integral approach

This study was published in Groundwater as:

___________________________________________________________________________

Sobhi Gollo, V., Broecker, T., Lewandowski, J., Nützmann, G. & Hinkelmann R. 2022. Flow

and transport modeling in heterogeneous sediments using an integral approach. Groundwater.

https://doi.org/10.1111/gwat.13275.

This is an open access article under the terms of the Creative Commons Attribution License,

which permits use, distribution and reproduction in any medium, provided the original work is

properly cited.

This is the peer reviewed version of the following article: Sobhi Gollo, V., Broecker, T.,

Lewandowski, J., Nützmann, G. & Hinkelmann R. 2022. Flow and transport modeling in

heterogeneous sediments using an integral approach. Groundwater., which has been published

in final form at https://doi.org/10.1111/gwat.13275. This article may be used for non-

commercial purposes in accordance with Wiley Terms and Conditions for Use of Self-Archived

Versions. This article may not be enhanced, enriched or otherwise transformed into a derivative

work, without express permission from Wiley or by statutory rights under applicable legislation.

Wiley’s version of record on Wiley Online Library and any embedding, framing or otherwise

making available the article or pages thereof by third parties from platforms, services and

websites other than Wiley Online Library must be prohibited.

___________________________________________________________________________

This is the postprint version of the article.

The test cases’ setups are listed in Appendix B (B9 one-dimensional tracer transport in

heterogenous groundwater, B10 transport simulations for heterogenous groundwater-surface

water interactions at rippled streambeds).

4. Flow and transport modeling in heterogenous sediments using integral approach

4.1 Abstract

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.

4.2 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, Pryschlak et al. 2015,

Lewandowski et al. 2019). This zone includes riverbeds, riverbanks, saturated sediments under

dry bars and riparian areas (Edwards 1998). A multi-directional, multi-scale 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, Pryschlak 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 to 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 and Cardenas 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. Elliot and Brooks 1997, Stonedahl et al. 2010, Broecker et al. 2019). Due to its

4. Flow and transport modeling in heterogenous sediments using integral approach

major impact on hyporheic exchange processes in streams, heterogeneity should be given high

priority in presence of morphology (Freeze and Cherry 1979, Pryschlak et al. 2015). In studying

hyporheic flow through morphological structures such as bed forms, heterogeneity adjests 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 modelled (Broecker et al. 2021).

Fig. 4-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 require 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.

Figure 4-1: Schematic view of hyporheic zone. qh – hyporheic flow from the main channel into

the streambed sediment and returning to the stream, qG – gaining flux from groundwater to

surface water, qL – losing flux from surface water to groundwater.

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 modelled the flume experiment of Fox et al. (2014) with a newly

developed integral approach which has advantages compared to 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

4. Flow and transport modeling in heterogenous sediments using integral approach

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. (2021 and 2022) have shown that compared to 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.

4.3 Materials and methods

4.3.1 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 to the original experimental setup

(6.4 m). The model width (in y-direction, see Fig. 4-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 (Fig. 4-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

4. Flow and transport modeling in heterogenous sediments using integral approach

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 72000 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 m2 compared to maximum of 8 × 10-3 m2 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.

Figure 4-2: Model setup and geometry for modeling flow and transport in heterogenous

sediment based on Fox et al. (2016); sediment heterogeneity is explained in Fig. 4-3.

4.3.2 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.

Fig. 4-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 Fig. 4-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

m3/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.

4. Flow and transport modeling in heterogenous sediments using integral approach

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

(qG) and losing (qL) conditions a 100 cm/d 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 m3/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.

Figure 4-3: (top) Grain diameters and (bottom) porosities.

In all cases pre-runs of 5 minutes with initial surface water flow of 0.15 m/s in x-direction are

carried out in order to reach quasi steady state. For losing and gaining cases, qL and qG are

applied as well. Results of these pre-runs 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 m2/s

enters the system through the inlet boundary.

4.3.3 Governing equations and computational tool

Figure 4-4: Model boundary conditions.

4. Flow and transport modeling in heterogenous sediments using integral approach

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. The equations for the conservation of mass and

momentum are written as ([ ]f is averaged parameter over void region):

φ𝛻𝛻[𝑣𝑣]𝑓𝑓= 0 (4-1)

φ�𝜕𝜕[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓

𝜕𝜕𝑡𝑡 +𝛻𝛻[𝜌𝜌]𝑓𝑓[𝑣𝑣]𝑓𝑓[𝑣𝑣]𝑓𝑓�=−φ𝛻𝛻[𝛻𝛻]𝑓𝑓+φ[𝜇𝜇]𝑓𝑓𝛻𝛻2[𝑣𝑣]𝑓𝑓+φ[𝜌𝜌]𝑓𝑓∙g + 𝛼𝛼 (4-2)

where 𝑣𝑣 [m/s] is the flow velocity, 𝑡𝑡 [s] is time, 𝛻𝛻 [Pa] is pressure, 𝜌𝜌 [kg/m3] is the density of

the fluids, 𝜇𝜇 [kg/(m s)] is the dynamic viscosity (𝜇𝜇= 𝜇𝜇𝑝𝑝ℎ + 𝜇𝜇𝑡𝑡; physical + turbulent viscosity),

g [m/s2] is the gravitational acceleration, φ [-] is the effective porosity and 𝛼𝛼 [kg/(m2s2)] 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

Eq. (3):

𝛼𝛼 = 𝐴𝐴 + 𝐵𝐵 (4-3)

𝐴𝐴 = −�150 (1−𝜑𝜑)

𝑑𝑑𝑝𝑝𝜑𝜑[𝜇𝜇]𝑓𝑓+ 1.75[𝜌𝜌]𝑓𝑓�[𝑣𝑣]𝑓𝑓��(1−𝜑𝜑)

𝑑𝑑𝑝𝑝[𝑣𝑣]𝑓𝑓 (4-4)

𝐵𝐵 = −0.34 (1−𝜑𝜑)

𝜑𝜑𝜕𝜕[𝑣𝑣]𝑓𝑓[𝜌𝜌]𝑓𝑓

𝜕𝜕𝑡𝑡 (4-5)

where 𝑎𝑎𝑝𝑝 [m] is the effective grain diameter.

Eq. (4-4) defines pressure loss due to fluid friction with sediment after Ergun (1952). In

sediment, compared to 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. Eq. (4-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 𝐴𝐴=𝐵𝐵=𝛼𝛼= 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 FVM (Finite Volume Method)

discretization scheme and for time it uses FDM (Finite Difference Method). The interface

between water and air phase is captured via the VoF (Vulume of Fluid) method. In OpenFOAM,

4. Flow and transport modeling in heterogenous sediments using integral approach

velocity is discretized using sparse linear solver (smoothSolver, smoother: symGaussSeidel).

Pressure is calculated by using GAMG (geometric agglomerated algebraic multigrid

preconditioner) solver. Coupling of velocity and pressure is done via the PIMPLE 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 [m2/s]:

𝜕𝜕[𝐶𝐶]𝑓𝑓

𝜕𝜕𝑡𝑡 +𝛻𝛻 ∙��𝑣𝑣𝑝𝑝

�



�𝑓𝑓[𝑣𝑣]𝑓𝑓�+𝛻𝛻∙([𝐷𝐷]𝑓𝑓𝛻𝛻[𝑣𝑣]𝑓𝑓)= 0 (4-6)

Where 𝑣𝑣𝑝𝑝

�



[m/s] is pore velocity, C is tracer concentration in [mg/l] and D [m2/s] is diffusion

coefficient. In 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 Sim

ulation; Smagorinsky

1963) which is more robust than the more commonly used RANS model for this study

Modeling the setup of Fig. 4-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-minute study of tracer propagation. Tracer expansion modelling with this

turbulence model led to comparable results as in laboratory experiment serving as a justification

for this approach.

4.4 Results and discussion

4.4.1 Validation of flow and transport in heterogenous sediment with a 1D analytical

solution

A one-dimensional 10 m long domain which is separated into two 5 m long sections with

porosities of 0.44 and 0.3 (compliant with Fox et al. 2016) is shown in Fig. 4-5. An inert tracer

with a diffusion coefficient of 10-9 m2/s is continuously injected from the left side of the domain.

Fig. 4-5 shows the tracer progression under constant velocity of 0.01 m/s. Simulation results

4. Flow and transport modeling in heterogenous sediments using integral approach

are compared to one-dimensional continuous injection analytical solutions of Kinzelbach

(1992). A very good agreement between analytical solution results and simulations can be seen

in Fig. 4-5 for different timesteps.

Figure 4-5: Comparison of integral model to analytical solutions for a continuous tracer

injection in a simple 1D model domain: Model domain and porosities (left top), tracer fronts

modelled with integral model after t=400 and 600 s (left bottom) and comparison of tracer

concentrations for different timesteps for the integral model and the analytical solution.

4.4.2 Comparison of numerical to experimental results for tracer transport

Figure 4-6: Tracer propagation for the investigated 6th ripple after 10, 40 and 60 min. Results

shown for both the lab experiment (a) and the integral model (b; purple). For 40 min, tracer

fronts simulated via MIN3P groundwater model (from Fox et al. 2016) are as well illustrated.

4. Flow and transport modeling in heterogenous sediments using integral approach

The ripples of the lab experiments deviated slightly from the shape at the start of the experiment

because of minor sediment movement during previous experiments.

In this section, for the transport of a tracer into heterogeneous sediment the lab observations of

Fox et al. (2016) are compared to the results of the integral model. Simulated tracer propagation

is compared to experiment in Fig. 4-6 under the 6th ripple. Numerical results of an alternative

model (MIN3P, Mayer et al. 2002; modelled by Fox et al. 2016) after 40 min tracer propagation

are also shown in Fig. 4-6 and compared to the integral model in Table 4-2.

Fig. 4-6 shows that the tracer first spreads around the green zone in the upper left with its small

hydraulic conductivity (smallest grain size diameter; Fig. 4-6, neutral 10 and 40 min, losing 10

min, gaining 10, 40 and 60 min). Compared to the integral modelling, 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.

Figure 4-7: Contribution of advection and diffusion to transport under different conditions at

60 min tracer propagation for the 6th ripple. Peclet number = Advective transport rate/diffusive

transport rate. In dark red areas (Peclet number > 1) advection is predominant. Other areas

(Peclet number < 1) are low-turbulent zones where diffusion is dominant.

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 Fig. 4-7 for tracer propagation

at 60 minutes. 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 (Fig. 4-7).

4. Flow and transport modeling in heterogenous sediments using integral approach

Figure 4-8: Differences of the real sediment setup in the experiment of Fox et al. (2016) (top

panel) and the originally planned sediment setup which was used in the integral model (middle

panel). The hyporheic zone beneath the 6th ripple, i.e. the discussed ripple in Fig. 4-9 is

magnified (lower panel) to display local sediment packing differences (𝑅𝑅𝑙𝑙). Note: compared to

Fig. 4-2, the first ripple from left is excluded from the model setup so it matches the image

provided by Fox et al. (2016).

Figure 4-8 shows how minor differences in dyed areas between experimental runs and modelled

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 modelled tracer fronts in heterogeneous sediments they

can be important. In Fig. 4-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 (𝑅𝑅𝑙𝑙) compared to 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 to the original plan the generally heterogeneous sediment pattern is

4. Flow and transport modeling in heterogenous sediments using integral approach

relatively intact (except 𝑅𝑅𝑙𝑙), a small bedform migration (𝑅𝑅𝑎𝑎1−5) and deformation (𝑅𝑅𝑡𝑡1−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 (𝑅𝑅𝑡𝑡) is also displayed in Fig. 4-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 Fig. 4-8.

Table 4-1: Tracer propagation comparison with time between integral model and experiment

under neutral, losing and gaining conditions.

To compare results of experimental runs and simulations quantitively, 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 4-1. Tracers are most quickly propagating under losing conditions and slowest

under gaining conditions. For gaining conditions the propagation rate drastically decreases with

time and results show very little further tracer propagation between 40 and 60 min (Table 4-1,

Fig. 4-9). RMSE (root mean square error) and N-RMSE (normalized RMSE; 0 < N-RMSE <

1; N-RMSE closer to 0 than 1 shows plausible fitting) values show that model and experiment

agree best for losing conditions 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

Dyed area [cm2]

Neutral

Losing

Gaining

Experiment

Simulation

Experiment

Simulation

Experiment

Simulation

10 min

45.726

43.114

58.175

59.013

26.920

25.737

40 min

92.620

84.713

135.711

136.145

55.077

57.323

60 min

101.042

96.360

200.417

197.353

55.580

57.431

RMSE

9.55

1.85

3.14

N-RMSE

0.17

0.01

0.1

4. Flow and transport modeling in heterogenous sediments using integral approach

undefined experimental setup modifications, they fairly match experimental tracer propagation

trends both qualitatively (Fig. 4-6) and quantitively (Table 4-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.

Table 4-2: Tracer propagation comparison after 40 minutes between the MIN3P and the integral

model (porousInter) under neutral, losing and gaining conditions.

Comparing the integral model results at 40 minutes to an alternative groundwater model

(MIN3P; Fig 4-6; obtained from Fox et al. 2016) in Table 4-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 to MIN3P), next, we will describe the modelled pressure distribution and velocity

fields in the heterogenous sediment.

4.4.3 Pressure and flow distribution in heterogenous sediment

Hydrodynamic pressure

Hydrodynamic pressure distributions under neutral, losing and gaining conditions on 4th, 5th

and 6th ripples are shown in Fig. 4-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 m3/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 m3/s through the bottom of the flume is applied equal to the flux used in Fox et

al. (2016).

Dyed area [cm2]

Neutral

Losing

Gaining

MIN3P

porousInter

MIN3P

porousInter

MIN3P

porousInter

84.873

84.713

136.030

136.145

56.007

57.323

4. Flow and transport modeling in heterogenous sediments using integral approach

Figure 4-9: Hydrodynamic pressure distribution under neutral, losing and gaining conditions

and the location of a horizontal cross section (z = 0.195 m, white line) as well as grain diameter

maps at 4th, 5th and 6th ripples.

Figure 4-10: Hydrodynamic pressures under neutral, losing and gaining conditions for a

horizontal cross section through the sediment bed in the depth marked in Fig. 4-9 as white line.

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 Fig. 4-9. Compared to 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 Fig. 4-9 are displayed in Fig. 4-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 Fig. 4-10,

4. Flow and transport modeling in heterogenous sediments using integral approach

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 Pa and 0.9 Pa, under losing conditions between 3 Pa and

8 Pa and under gaining conditions between 0.3 Pa 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. Fig. 4-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 to 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 to

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 Fig. 4-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 Fig. 4-3).

4. Flow and transport modeling in heterogenous sediments using integral approach

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 4th ripple is allowed to enter deeper into the sediment compared to the 5th

and the 6th 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 5th and the 6th ripple. At less intense (v > 1 × 10-

5 m/s) flow under neutral conditions, a blue block underneath the 5th ripple and partly the 6th

ripple has limited hyporheic flow to the upper 5-6 cm as opposed to the 4th ripple in which flow

can enter deeper into the sediment. Under the 4th 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 5th and 6th ripples, for hyporheic flow through the 4th ripple longer residence times

occur.

Figure 4-11: Flow vectors (dimensionless) mapped over grain size diameters of the sediment

underneath 4th, 5th and 6th ripples for v > 1 × 10-3 m/s, v > 1 × 10-4 m/s and v > 1 × 10-5 m/s for

neutral, losing and gaining conditions.

As was previously mentioned, hyporheic zone is thicker under losing conditions compared to

gaining and neutral conditions. Under losing conditions and at higher flow rates (v > 1 × 10-4

m/s), deflected hyporheic flow under the 4th 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

4. Flow and transport modeling in heterogenous sediments using integral approach

beneath the 5th and 6th ripple is being pulled down left (under 5th ripple) and right (under 6th

ripple) of the blue block (Fig. 4-11). The low grain size diameter/low porosity blue block under

the 5th 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 to 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 4th ripple due

to higher grain size diameter/higher porosity underneath this ripple, results in more impact of

gaining conditions on the hyporheic flow compared to the 5th 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, i.e. neutral, gaining and losing.

Compared to 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 to 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.

4. Flow and transport modeling in heterogenous sediments using integral approach

4.5 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, i.e. a laboratory flume. Even in the presence of groundwater flow and with

complex streambed heterogeneity, the modelled 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 (Broecker

et al. 2019) and transport through homogeneous sediment under ambient groundwater

conditions (Broecker et al. 2021). 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.

5. Supplementary contributions

5.1 Integral Flow simulations for groundwater-surface water interactions at

rippled streambeds

This study was published in Water as:

___________________________________________________________________________

Broecker, T., Teuber, K., Sobhi Gollo, V., Nützmann, G., Lewandowski, J. & Hinkelmann, R.

(2019). Integral Flow Modelling Approach for Surface Water-Groundwater Interactions along

a Rippled Streambed. Water, 11(7), 1517, https://doi.org/10.3390/w11071517.

article distributed under the terms and conditions of the Creative Commons Attribution

license (http://creativecommons.org/licenses/by/4.0/). The definitive peer-reviewed and edited

version of this article is published in doi:10.3390/w11071517 and available at

www.mdpi.com/journal/water

___________________________________________________________________________

This is the abstract of the article (postprint).

Abstract

Exchange processes of surface and groundwater are important for the management of water

quantity and quality as well as for the ecological functioning. In contrast to most numerical

simulations using coupled models to investigate these processes, we present a novel integral

formulation for the sediment-water-interface. The computational fluid dynamics (CFD) model

OpenFOAM was used to solve an extended version of the three-dimensional Navier–Stokes

equations which is also applicable in non-Darcy-flow layers. Simulations were conducted to

determine the influence of ripple morphologies and surface hydraulics on the flow processes

within the hyporheic zone for a sandy and for a gravel sediment. In- and outflowing exchange

fluxes along a ripple were determined for each case. The results indicate that larger grain size

diameters, as well as ripple distances, increased hyporheic exchange fluxes significantly. For

higher ripple dimensions, no clear relationship to hyporheic exchange was found. Larger ripple

5. Supplementary contributions

lengths decreased the hyporheic exchange fluxes due to less turbulence between the ripples. For

all cases with sand, non-Darcy-flow was observed at an upper layer of the ripple, whereas for

gravel non-Darcy-flow was recognized nearly down to the bottom boundary. Moreover, the

sediment grain sizes influenced also the surface water flow significantly.

5.2 Transport simulations for groundwater-surface water interactions under

neutral, losing and gaining flow conditions for homogenous streambed

This study was published in Groundwater as:

_______________________________________________________________

Broecker, T., Sobhi Gollo, V., Fox, A., Lewandowski, J., Nützmann, G., Arnon, S. &

Hinkelmann, R. (2021). High Resolution Integrated Transport Model for Studying Surface

Water–Groundwater Interaction, Groundwater, https://doi.org/10.1111/gwat.13071.

National Ground Water Association. This article is an open access article distributed under

the terms and conditions of the Creative Commons Attribution license

(https://creativecommons.org/licenses/by-nc-nd/4.0/). The definitive peer-reviewed and edited

version of this article is published in https://doi.org/10.1111/gwat.13071.

___________________________________________________________________________

This is the abstract of the article (postprint).

Abstract

Transport processes that lead to exchange of mass between surface water and groundwater play

a significant role for the ecological functioning of aquatic systems, for hydrological processes

and for biogeochemical transformations. In this study, we present a novel integral modeling

approach for flow and transport at the sediment-water interface. The model allows us to

simultaneously simulate turbulent surface and subsurface flow and transport with the same

conceptual approach. For this purpose, a conservative transport equation was implemented to

an existing approach that uses an extended version of the Navier-Stokes equations. Based on

previous flume studies which investigated the spreading of a dye tracer under neutral, losing

5. Supplementary contributions

and gaining flow conditions the new solver is validated. Tracer distributions of the experiments

are in close agreement with the simulations. The simulated flow paths are significantly affected

by in‐ and outflowing groundwater flow. The highest velocities within the sediment are found

for losing condition, which leads to shorter residence times compared to neutral and gaining

conditions. The largest extent of the hyporheic exchange flow is observed under neutral

condition. The new solver can be used for further examinations of cases that are not suitable for

the conventional coupled models, for example, if Reynolds numbers are larger than 10.

Moreover, results gained with the integral solver provide high resolution information on

pressure and velocity distributions at the rippled streambed, which can be used to improve flow

predictions. This includes the extent of hyporheic exchange under varying ambient groundwater

flow conditions.

6. Synthesis

6.1 Summary and conclusions

Due to a recent shift of perception towards considering groundwater and surface water as a

single resource, many field and laboratory techniques utilizing advanced measurement methods

are developed. Nevertheless, numerical models are not keeping pace with these developments

and look upon groundwater and surface water separately. The present thesis is an endeavor to

introduce, validate, apply and further develop an integral approach to capture continuous

transport and exchange between groundwater and surface water. For that matter, using a single

set of mathematical equations (modified three-dimensional Navier-Stokes equations),

groundwater and surface water are modeled in one domain. Simulations are based on CFD

software OpenFOAM. The benefit of utilizing this software is spatiotemporally high-resolution

modeling of groundwater-surface water exchange processes which however limits the spatial

scale due to its high computational effort.

6.1.1 General outcomes

Prior to detailed description of the outcomes of model comparison as well as validation and

further development of the integral approach, a general summary of this thesis is put forward.

In sections 2 and 3, integral approach is compared to coupled approach for two cases: flow and

exchange across rippled streambed and ventilation in U-shaped burrows:

• For both cases, steep velocity gradients and pressure variations occur at the immediate

interface between groundwater and surface water. The comparison of flow and pressure

results indicate that:

o In the case of flow across rippled streambed, by defining identical overlying surface

water flow velocity and water table for coupled and integral models, identical total

pressure fields in surface water are achieved. However, it is displayed that integral

approach shows smaller velocities and hydrodynamic pressures in the vicinity of the

ripples.

o In the case of ventilation through U-shaped burrow, by defining a fully saturated

sediment and burrow, identical hydrostatic pressure fields in the domain are

achieved. In the pumping zone however, identical hydrodynamic pressures were

only obtained when applying a higher flow velocity in the integral model compared

to coupled model. On the other hand, if identical flow velocities were exhibited,

6. Synthesis

integral model’s hydrodynamic pressure was smaller than the one of the coupled

model.

• In the immediate interface between groundwater and surface water, a continuous exchange

of flow occurs. Therefore, considering feedback between groundwater and surface water is

crucial:

o In the flow across rippled streambed case, surface water entering the rippled area

can continuously flow in from one side of the ripple and flow out at the other side

using the integral model. On the contrary, using coupled approach, the effect of flow

over the ripple is translated into pressure distributions to be transferred into the

groundwater model.

o In the case of ventilation through U-shaped burrow, unlike coupled model, a

continuous flow from inlet to outlet branch through the sediment is captured via

integral model.

• Very high spatiotemporal resolution of the integral approach authorizes its application for

small scale groundwater-surface water interface phenomena. Its high computational burden

should however be considered:

o In the flow across rippled streambed case, on account of a comparison with a simpler

shallow water equation model, the use of a CFD Navier-Stokes model for modeling

the surface water part of the coupled model is proved essential. Therefore,

computational resources required for the coupled and integral models are in fact

similar as they both require high performance computing clusters.

o In the ventilation through U-shaped burrow case, due to symmetric geometry of the

burrow, in the coupled model Navier-Stokes equations are applied for the flow

through the burrow without using a CFD model. Consequently, computational

burden of the integral model is relatively higher than the coupled model.

In section 4, the integral approach which was formerly validated for transport through

homogenous sediment by Broecker et al. (2021) is validated for modeling tracer transport

through heterogenous sediment. Through comparison with analytical solutions and

experimental results, it is demonstrated that:

• Complex hydraulics of flow and tracer propagation through a heterogenous sediment

can be plausibly captured using the integral approach.

This study was embedded in the collaborative research common topic group of “Interface

Hyporheic zones” of the 2nd cohort of the Urban Water Interfaces Research Training Group

(UWI). Next to the model development provided in this study, a series of field and laboratory

6. Synthesis

experiments which examined different biogeochemical processes of the hyporheic zone and

bank filtration processes are provided in other projects.

6.1.2 Outcomes of comparing integral to coupled approach for flow across rippled

streambed

For the case of flow across a rippled streambed, an integral model (using modified Navier-

Stokes equations and LES turbulence model) to model flow and exchange processes in

groundwater, surface water and groundwater-surface water interface is compared to a coupled

model of groundwater (using Richards equation model) and surface water (using Navier-Stokes

equations and LES turbulence model) flow. Furthermore, a shallow water equations model is

used to model flow in surface water. For accurate compatibility of the integral model to the

coupled model, next to using identical model domain details like model length, mesh resolution,

weir placement and ripple geometry in surface water part of both models, same flow conditions

are applied which result in very similar water tables and pressure distributions. However,

allowing continuous exchange near interface in integral model yields velocity distributions very

dissimilar to coupled model. Near the interface in surface water, velocity vectors meet the ripple

with different impact angles. In coupled model, these vectors are deflected as they cannot enter

the rigid bottom boundary of the surface water model causing turbulent eddies to form. In

integral model however, the interface is connected to the porewater and depending on the

sediment porosity, a portion of the flow enters the sediment causing considerably less deflection

and formation of less turbulence near the interface in surface water. In groundwater, coupled

model is bound to flow velocities within Darcy’s law. Contrarily, integral model can detect

areas where flow as the result of close contact to surface water has velocities outside of the

range defined by Darcy’s law. As continuous flow across interface in connected groundwater-

surface water streams constantly occurs, results of integral model for flow and exchange

processes in this zone are more plausible compared to coupled model and application of coupled

models on the interface results in overestimation of flow fields in interface region. The extent

of discrepancies between integral and coupled models can depend on factors like stream flow

velocity, streambed morphology, sediment characteristics and its heterogeneity.

Next to more plausible estimation of flow fields near interface, the use of the modified Navier-

Stokes equations for modeling groundwater-surface water exchange has further advantages.

Under proper hydrogeologic conditions (such as surface water flow, sediment porosity,

streambed morphology) surface water eddies are to some extent transferred into the sediment

close to the interface. Unlike coupled model, the equations governing the integral model can

6. Synthesis

detect these eddies. Nevertheless, the extent of turbulence penetration in relation to the

hydrogeological conditions needs to be further verified for the mathematical model of the

integral approach.

The residence time of the hyporheic water is a crucial topic in investigations on water quality

and purification potential of the hyporheic zone. The integral model continuously traces the

hyporheic water flow conditions and its pathways as it enters and leaves the hyporheic zone.

Such capability is far more superior to coupled models which only function by transferring the

resulting pressure distributions over the common boundaries from one model to the other.

By comparing surface water model outcomes of Navier-Stokes equations-based model to

shallow water equations-based model, application of the Navier-Stokes equations-based model

is proven necessary. Consequently, for modeling flow across rippled streambed, applying

computationally expensive CFD models is necessary for both approaches. Although run time

of the integral model is a few hours higher than of the surface water part of the coupled model,

both need access to supercomputers. Concurrently, modeling results are very similar a few

decimetres away from the interface. Although the difference in computational requirement is

not significant in this study, the outcomes suggest that if investigations for deeper groundwater

are intended, coupled model is computationally more efficient compared to integral model.

Detailed information of this study can be found in section 2.

6.1.3 Outcomes of comparing integral to coupled approach for ventilation through U-

shaped burrow

In the next step, by examining a distinctive case, application area of the integral approach is

extended and further advantages to coupled approach are discussed. Here, a case of flow

initiated by pumping activities inside a U-shaped burrow (by bioturbation of the tube-dwelling

invertebrates) that generates flow and exchange across surrounding sediment (porewater)-

burrow (surface water) interface is modelled with integral solver and compared to an existing

identical coupled model. For accurate compatibility of the integral model to the coupled model,

identical model domain details like burrow geometry, sediment characteristics and flow source

placement (at the pumping zone) are carefully implemented. It was however realized that to

generate the same pressure conditions at the pumping zone in the integral model, higher

velocities compared to coupled model should be deployed. Similar to rippled streambed case,

here allowing continuous exchange of flow between burrow and surrounding sediment by using

integral model results in unequal pressure-velocity fields in the pumping zone compared to

coupled model. In coupled model, the generated pressure (resulting from solving Navier-Stokes

6. Synthesis

equation for surface water flow) in the burrow which is used as coupling parameter is

independent from the sediment characteristics. By continuous burrow-porewater exchange in

the integral model, sediment characteristics participate in determining flow and pressure

conditions. Integral modeling of burrow ventilation processes with different flow velocities and

grain diameters show that higher velocities and smaller grain diameters result in higher

pressures in the pumping zone (which then spreads to the whole domain).

Next to plausible modeling of flow and pressure fields, integral model allows continuous

investigation of the transport of flow through the sediment, burrow and on the overlying water

and their interactions. As biochemical composition of water in the sediment, in the overlying

water and on the burrow branches are often different, monitoring their interactions gives insight

to better understanding the hydro-biogeochemical processes across tube-dwelling invertebrates

habitat.

Although required computational capacity of the integral model is higher than of the coupled

model in this case, negligence of important hydrogeological parameters in determining

velocity-pressure fields using the coupled model points towards the absolute necessity of using

integral model for plausible results acquirement. It should be noted that the velocity values used

in both coupled and integral numerical simulations are limited to the averages over pumping

and resting periods and do not mimic the actual generated velocities during burrow ventilation.

Detailed information of this study can be found in section 3.

6.1.4 Outcomes of transport modeling for heterogenous groundwater-surface water-

interactions

After gathering significant evidence of superiority of the integral approach to prevalent methods

for modeling groundwater-surface water interactions and pointing out its ability to track flow

across the interface, it is very beneficial to extend the solver for modeling transport processes

that it can step into the vast topic of contamination remediation, nutrient transformation and

water quality modification. The initial validation step for transport across homogenous

streambed (now backed by profound flow and exchange model validation and comparison

carried out in this thesis) has been taken as stated in section 5.2. Next to streambed morphology

(studied in 5.1), ambient groundwater conditions and sediment heterogeneity are decisive

factors in tracking transport in the sediment. Throughout this step, the integral solver, formerly

validated for transport through homogenous sediment by Broecker et al. (2021), was extended

for advective-diffusive transport in heterogenous sediment and under ambient groundwater

conditions. The mathematical model was initially validated next to analytical solutions. A flume

6. Synthesis

experiment conducted to track the tracer propagation through heterogenous sediment was

modeled via integral approach. Velocity and hydrodynamic pressure results indicated that under

losing condition more water is transported through the sediment compared to neutral and

gaining conditions resulting in an extended hyporheic zone. The transported water is deflected

in various areas as it reaches stretches with varying porosities and grain diameters resulting in

very different residence times. The collective impact of sediment heterogeneity and ambient

groundwater conditions is higher for losing condition as water meets more layers of

heterogeneity while entering the sediment. These findings were validated by comparing tracer

propagation in the model to the experiment. Plausible comparison of results authorizes the use

of the integral solver for transport modeling in heterogenous sediment and determining the

effect of sediment heterogeneity under ambient groundwater conditions on flow intensity,

direction and residence times. Detailed information of this study can be found in section 4.

6.1.5 Recommendations for modeling and data processing

Qualitative advantages of the integral approach compared to coupled approach should be

investigated next to comparing velocity-pressure and intrinsic computational effort differences

of the models. This study carries out pre-processing, processing and post-processing of two

methodically distinctive models. Here a few additional notes for model selection are presented:

• Groundwater models that resolve seepage flow in x-z axis and are capable of meshing

complex structures like ripples are usually not freely available. For this study, a

commercial groundwater model was licenced through scientific collaborations.

• Utilizing a coupled model requires a very good knowledge in modeling both

groundwater and surface water processes as well as coupling techniques.

• Groundwater and surface water models used in coupled approach oftentimes use

different pre- and post-processing tools. For instance, generating a numerical mesh with

identical mesh nodes on the boundaries is sometimes challenging. Applying the resulted

pressure distribution from the surface water model to the groundwater model is very

dependent of precise placement of nodes.

• In comparing necessary time to run coupled and integral models that essentially use the

same computational resource (inevitability to access supercomputers), the coupling step

is often neglected. Depending on the models used for coupling, this step can be time-

consuming.

All the discussions above indicate that integral approach is a very powerful method and should

be prioritized to coupled approach in addressing small-scale high-resolution groundwater-

6. Synthesis

surface water interaction processes. This approach is now capable of modeling flow and

transport of tracers from the surface water to heterogenous groundwater. As it was formerly

mentioned, the advantages of the integral approach vanish in modeling very deep groundwater

and its use is only recommended where the effects of groundwater-surface water interaction are

present. By comparing k-ε and LES turbulence models for the case of transport in the

heterogenous sediment, it was realized that the use of LES turbulence model for accurate

modeling of such complex case is necessary. For simpler cases like homogenous sediment (see

case in section 5.2), computational effort can decrease by choosing k-ε as the turbulence model.

With the use of 96 processors of the HLRN supercomputer, each 1 minute of model run for the

heterogenous streambed case took ~8 hours. With bigger domain and by including other factors

like overlying air, higher run times are expected. It is therefore advised to limit model

dimensions to be as small as possible, while still representative, and consider modeling options

such as symmetry planes (as used in U-shaped burrow case) that can noticeably reduce

computational effort.

Although many cases and aspects have been studied in this thesis, there are cases such as

breakwater structure flow and aspects such as sediment transport and reactive transport that can

draw the next steps of model application and development. These are explained in the next

section.

6.2 Outlook

To apply the integral approach to a wider range of scenarios, modeling of wave intensity

reduction through shoreline breakwater structures is recommended. This can be done by

investigating wave and turbulence penetration under different hydrogeological conditions.

Validation of wave and turbulence penetration can be done by comparing model results to

experiments of Roche et al. (2018) and Blois et al. (2014) who investigated the depth in which

turbulence can penetrate sediment with relatively large grain size. As for this case flow

fluctuations are a decisive factor, transient flow conditions (unlike what was assumed in the

previous studies) can be applied. The integral solver is currently embedded in OpenFOAM

2.4.0 version. A wider and more efficient range of libraries, boundary conditions and tools are

available in the newer versions of OpenFOAM which can be very useful in modeling wave

impact on breakwater structures. The integral solver has been recently validated which opens

up the potential to use the extended libraries and boundary conditions of the newest

OpenFOAM version.

6. Synthesis

Using integral solver in section 2 and under assumption of no bedform movement, a small area

of very high Reynolds number has been detected. Depending on sediment coherency, at a

certain Reynolds number in the sediment (under increasing surface water flow velocities),

bedform can presumably move. For sandy sediment, it is very beneficial to use the integral

solver to realize the maximum Reynolds numbers where sediment transport starts. This

can be done by sediment characteristics sensitivity analysis and modeling sediment transport

experiment setups. As long as the realized maximum Reynolds numbers are exceeded for a cell

in the model, sediment transport equations can be applied to them. As a basis “SedFoam”

(Chauchat et al. 2017) can be used to include sediment transport processes into the integral

solver.

To this date, advective-diffusive tracer propagation for heterogenous sediment and under

ambient groundwater condition has been verified and validated for the integral solver. The next

step in developing the solver is to include and validate it for advective-diffusive-dispersive

transport. In Appendix A, the solver is mathematically verified for 2D advective-diffusive-

dispersive transport next to analytical solutions of Kinzelbach (1992). For validation, a

comparison to experiments which include dispersion into their transport investigations is

recommended. Thereafter, reactive components of transport across groundwater-surface water

interface like phosphorus, nitrogen, dissolved oxygen and organic carbon can be included. As

a basis “BioChemFOAM” (Hernandez Murcia 2014) can be used to include reactive transport

processes into the integral solver.

Further development and validation of the integral solver through using experimental results

via collaborative research in the 3rd cohort of UWI is strongly advised. A detailed model of

flow and transport processes benefits field and laboratory experiment by providing indication

of flow rates and pathways which are crucial in determining biogeochemical processes.

Appendix A: Validation of porousInter for dispersive transport

This study was partly presented in Groundwater Quality 2019 conference as:

___________________________________________________________________________

Sobhi Gollo, V., Broecker, T., Lewandowski, J., Nützmann G. and Hinkelmann, R. (2019).

Integral modelling approach for surface water-groundwater interface domain

___________________________________________________________________________

This is the extended version of the presented study.

The study in section 4 as well as Broecker et al. (2021) who have modeled transport of tracer

through homogenous sediment, have both neglected dispersion. In addition, validation in

section 6 as well as in Broecker et al. (2021) has been done by using a 1D analytical solution.

Here, by implementing the general balance equation for transport of tracer with initial

concentration of “C” into the porousInter solver, we validate the solver for a 2D dispersive

tracer by a comparison to 2D analytical solutions of Kinzelbach (1992) for further use in future

cases. The general balance equation is written as follows:

𝜕𝜕(𝜑𝜑𝜌𝜌𝑤𝑤𝐶𝐶)

𝜕𝜕𝑡𝑡 +𝛻𝛻



[𝑣𝑣 𝑣𝑣 − 𝜑𝜑𝐷𝐷ℎ𝑦𝑦𝑑𝑑𝛻𝛻(𝜌𝜌𝑤𝑤𝑣𝑣)] = 𝑞𝑞𝐶𝐶 Eq. (A-1)

Where 𝜌𝜌𝑤𝑤 [kg/m3] is density of water, 𝑞𝑞𝐶𝐶 [kg/(m3s)] is sink or source of concentration and 𝐷𝐷ℎ𝑦𝑦𝑑𝑑

[m2/s] is the hydrodynamical dispersion tensor and is written as follows:

𝐷𝐷ℎ𝑦𝑦𝑑𝑑 = 𝐷𝐷𝑐𝑐𝑚𝑚𝑙𝑙 +𝐷𝐷𝑐𝑐𝑚𝑚𝑐𝑐ℎ=𝐷𝐷𝑐𝑐𝑚𝑚𝑙𝑙Ɨ+Đ𝑐𝑐𝑚𝑚𝑐𝑐ℎ=�𝐷𝐷𝑥𝑥𝑥𝑥 𝐷𝐷𝑥𝑥𝑦𝑦

𝐷𝐷𝑦𝑦𝑥𝑥 𝐷𝐷𝑦𝑦𝑦𝑦� Eq. (A-2)

Đ𝑐𝑐𝑚𝑚𝑐𝑐ℎ[m2/s] is mechanical dispersion tensor and Ɨ is the unit tensor. The components of the

dispersion tensor on the right side of Eq. (A-2) are written as follows:

𝐷𝐷𝑥𝑥𝑥𝑥 = 𝛼𝛼𝜌𝜌𝑣𝑣𝑥𝑥

|𝑣𝑣| + 𝛼𝛼𝑇𝑇𝑣𝑣𝑦𝑦

|𝑣𝑣| + 𝛼𝛼𝑇𝑇𝑣𝑣𝑧𝑧

|𝑣𝑣| + 𝐷𝐷𝑐𝑐𝑚𝑚𝑙𝑙

𝐷𝐷𝑦𝑦𝑦𝑦 = 𝛼𝛼𝑇𝑇𝑣𝑣𝑥𝑥

|𝑣𝑣| + 𝛼𝛼𝜌𝜌𝑣𝑣𝑦𝑦

|𝑣𝑣| + 𝛼𝛼𝑇𝑇𝑣𝑣𝑧𝑧

|𝑣𝑣| + 𝐷𝐷𝑐𝑐𝑚𝑚𝑙𝑙 Eq. (A-3)

𝐷𝐷𝑥𝑥𝑦𝑦 = 𝐷𝐷𝑦𝑦𝑥𝑥 = (𝛼𝛼𝜌𝜌 - 𝛼𝛼𝑇𝑇)𝑣𝑣𝑥𝑥 𝑣𝑣𝑦𝑦

|𝑣𝑣|

Appendix A: Validation of porousInter for dispersive transport

In Eq. (A-3) 𝛼𝛼𝜌𝜌 [m] is longitudinal and 𝛼𝛼𝑇𝑇 [m] is transversal dispersion length. These two

parameters are defined using 𝐷𝐷𝜌𝜌 [m2/s] (longitudinal dispersion coefficient) and 𝐷𝐷𝑇𝑇 [m2/s]

(transversal dispersion coefficient):

𝛼𝛼𝜌𝜌 = 𝐷𝐷𝜌𝜌/𝑣𝑣 , 𝛼𝛼𝑇𝑇 = 𝐷𝐷𝑇𝑇/𝑣𝑣 Eq. (A-4)

Figure A-1: Model domain of 2D dispersive transport validation. Small red square centered at

x = 1 m and y = 5 m shows the location of initial pulse injection with C = 1 mg/l. Flow is 1D

in x-direction.

Figure A-2: Model of dispersive pulse tracer transport progression using porousInter after 1 s,

50 s and 100 s.

Appendix A: Validation of porousInter for dispersive transport

A two-dimensional 10×10 m2 domain which is filled with a sediment with porosity of 0.3 is

shown in Figure A-1. A tracer with diffusion coefficient of 10-9 m2/s, longitudinal dispersion

length of 0.3 m and transversal dispersion length of 0.03 m is injected over the area of 1 m ≤ x

≤ 1.1 m and 4.95 m ≤ y ≤ 5.05 m. We display the simulated tracer propagation after 1, 50 and

100 s in Fig. A-2. As the tracer progresses forward with the flow, dispersion results in the

extension of the affected area. Higher longitudinal than transversal dispersion results in wider

tracer propagation in the 𝑥𝑥 direction. For each timestep, in Fig. A-3, we compare the simulation

results (left) to the analytical solutions (right).

Figure A-3: Simulated (left) and analytical (right) tracer concentrations after 1, 50 and 100 s

over the model domain. Tracer concentration is displayed in the vertical axis.

Appendix A: Validation of porousInter for dispersive transport

Furthermore, we calculate the correlation between the simulated and analytical tracer

propagations by looking at the tracer concentrations across several cross sections (Fig. A-3).

These cross sections are highlighted with the dotted red lines in Fig. A-3. For all three timesteps,

cross-sections are located at 𝑦𝑦 = 5 m. Perpendicular to these cross sections, for t = 1 s, 50 s and

100 s respectively, cross-sections are located at 𝑥𝑥 = 1 m, 6 m and 10 m. Both cross sections

intercept at the maximum tracer concentration in each timestep.

Figure A-4: Comparison of the simulated and analytical tracer propagation after 1, 50 and 100

s for cross sections in 𝑥𝑥 and 𝑦𝑦 directions. These cross-sections are highlighted as red dotted

lines in Fig. A-3.

In Fig. A-4, we used root mean squared error (RMSE) and normalized-root mean squared error

(N-RMSE) as correlation variables to verify the likeness of our simulations to the Kinzelbach

(1992) analytical solutions. Both variable values are very low (RMSE << tracer concentration

ranges; N-RMSE << 1) for all the cross-sections confirming a very good agreement between

the simulations and analytical solutions. Therefore, the extended mathematical code of the

integral approach is affirmed to be capable of plausibly simulating the dispersion of a tracer in

the porous medium.

Appendix B: Test case overview

B1 Seepage through a dike modeled with porousInter and PCSiWaPro®

Table B 1: Model setup for test case B1 (seepage through a dike modeled with porousInter and

PCSiWaPro®)

General

Referred to in section

References

Published in

2.7

Casagrande (1937), Lattermann (2010)

Broecker et al. (2019), Sobhi Gollo et al. (2022a)

Domain discretization

Dimensions

Mesh generator

Number of cells

1*) two-dimensional (length: 120 m/1.2 m, height: 25 m/2.5 m),

2*) two-dimensional (length: 8.7 m, height: 2.2 m)

1) blockMesh, 2) PCSiWaPro mesh generator

1) 75000, rectangular, 2) 27728, triangular

Turbulence model

1) laminar, 2) -

Hydrodynamic simulations

Solver

Time step

Simulation time

1) porousInter, 2) PCSiWaPro® (Richards equations)

1) variable, 2) min: 0.0001 s, max: 1 s

1) 9000 s, 2) 3 h

Transport simulations

Not performed

*1 = porousInter, 2 = PCSiWaPro®

Appendix B: Test case overview

Table B 2: Boundary conditions for test case B1 - porousInter

alpha.water (-)

p_rgh (kg/(ms²))

U (m/s)

effPackingRadius (m)

porosity (-)

inlet_water

inletOutlet

inletValue: uniform 1

value: uniform 1

fixedValue

value uniform 186390/18639

zeroGradient

inlet_air

inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

outlet

zeroGradient

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

atmosphere

inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

pressureInletOutletVelocity

value uniform (0 0 0)

zeroGradient

ground

zeroGradient

fixedFluxPressure

value uniform 0

fixedValue

value: uniform (0 0 0)

zeroGradient

sidewalls

empty

initial

conditions

uniform 0

define initial water level of

19 m/1.9 m the first 100 m

/10 m using boxToCell in

setFields

(volScalarFieldValue

alpha.water 1)

uniform 0

define initial water level of

19 m/1.9 m the first 100

m/10 m using boxToCell in

setFields

(volScalarFieldValue p_rgh

186390/18639)

uniform (0 0 0)

uniform 0.0159/0.002

uniform 1

in the

sediment 0.25

defined by

setFields

Appendix B: Test case overview

100

Table B 3: Boundary conditions for test case B1 – PCSiWaPro®.

left

Linear pressure head distribution to account for constant water level of 1.9 m

right

Free drainage

top

Free drainage

bottom

No flow

initial

conditions

Dam fully saturated (water content = 0.25), material: DIN 4220 pure sand, porosity is

manually set to 0.25

Appendix B: Test case overview

101

B2 Seepage through a rectangular dam modeled with porousInter and

PCSiWaPro®

Table B 4: Model setup for test case B2 (Seepage through a rectangular dam modeled with

porousInter and PCSiWaPro®).

General

Referred to in section

References

Published in

2.7

Westbrook (1985), Aitchison and Coulson (1972), Kobus

and Keim (2001), Di Nucci (2015)

Broecker et al. (2019), Sobhi Gollo et al. (2022a)

Domain discretization

Dimensions

Mesh generator

Number of cells

1*) two-dimensional (length: 250 m, height: 48 m),

2*) two-dimensional (length: 16 m, height: 24 m)

1) blockMesh, 2) PCSiWaPro mesh generator

1) 191744, 2) 14266

Turbulence model

1) laminar, 2) -

Hydrodynamic simulations

Solver

Time step

Simulation time

1) porousInter, 2) PCSiWaPro® (Richards equations)

1) variable, 2) min: 0.0001 s, max: 1 s

1) 9000 s, 2) 3 h

Transport simulations

Not performed

1 = porousInter, 2 = PCSiWaPro®

Appendix B: Test case overview

102

Table B 5: Boundary conditions for test case B2 - porousInter

alpha.water [-]

p_rgh [kg/(ms²)]

U [m/s]

effPackingRadius [m]

porosity [-]

inlet_water

inletOutlet

inletValue: uniform

value: uniform

fixedValue

value uniform

235440

zeroGradient

inlet_air

inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

outlet

zeroGradient

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

atmosphere

inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

pressureInletOutletVelocity

value uniform (0 0 0)

zeroGradient

ground

zeroGradient

fixedFlusPressure

value uniform 0

fixedValue

value: uniform (0 0 0)

zeroGradient

sidewalls

empty

initial

conditions

uniform 0

define initial water

level of 24 m at the left side

of the dam, to 4 m at the right side and stepwise

from 24 to 4 m inside the dam

using boxToCell

in setFields

(volScalarFieldValue

alpha.water 1)

uniform 0

define

pressure according

initial water

level using

boxToCell in

setFields

uniform (0 0 0)

uniform 0.0159/0.002

uniform 1

in the sediment

0.25 defined by

setFields

Appendix B: Test case overview

103

Table B 6: Boundary conditions for test case B2 – PCSiWaPro®

left

Linear pressure head distribution to account for constant water level of 24 m

right

Linear pressure head distribution to account for constant water level of 4 m

top

Free drainage

bottom

No flow

initial

conditions

Dam fully saturated (water content = 0.25), material: DIN 4220 pure sand, porosity is

manually set to 0.25

Appendix B: Test case overview

104

B3 One-dimensional one-phase modeling of surface water flow over rippled

streambed with hms (shallow water equations)

Table B 7: Model setup for test case B3 (surface water flow over ripples using hms)

General

Referred to in section

References

Published in

2.3.5

Broecker et al. (2019)

Sobhi Gollo et al. (2022a)

Domain discretization

Dimensions

Mesh generator

Number of cells

one-dimensional (length: 15 m)

hms

1500

Turbulence model

laminar

Hydrodynamic simulations

Solver

Time step

Simulation time

hms (shallow water equations)

forward_Euler, 10 s (max)

1200 s

Transport simulations

Not performed

Table B 8: Boundary conditions for test case B3 – PCSiWaPro®.

Inlet (left)

HSurfaceFlowOrthogonalUnitDischargeValues: Volumetric inflow rate: 0.5 m3/s

Outlet (right)

HSurfaceFlowFreeOutflowValues: free outflow

bottom

HSurfaceFlowNoFlowValues: No flow

initial conditions

0.5 m water table for the entire domain, Manning coefficient: 0.03 s/m1/3 to account for clean and

straight natural streams

Appendix B: Test case overview

105

B4 Three-dimensional two-phase modeling of surface water flow over rippled

streambed with interFoam

Table B 9: Model setup for test case B4 (surface water flow over ripples using interFoam)

General

Referred to in section

References

Published in

2.3.5

Broecker et al. (2019)

Sobhi Gollo et al. (2022a)

Domain discretization

Dimensions

Mesh generator

Number of cells

three-dimensional (length: 15 m, height: 1 m, depth: 1 m)

gmsh

1160900

Turbulence model

LES

Hydrodynamic simulations

Solver

Time step

Simulation time

interFoam

variable

60 s

Transport simulations

Not performed

Appendix B: Test case overview

106

Table B 10: Boundary conditions for test case B4

alpha.water [-]

p_rgh [kg/(ms²)]

U [m/s]

inlet_water

inletOutlet

inletValue: uniform

value: uniform

fixedFluxPressure

value: uniform

flowRateInletVelocity

volumetricFlowRate 0.

5 [m3/s]

value

uniform (0 0 0)

inlet_air

inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

outlet

zeroGradient

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

atmosphere

inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

pressureInletOutletVelocity

value

uniform (0 0 0)

walls

zeroGradient

fixedFluxPressure

value

: uniform 0

fixedValue

value: uniform (0 0 0)

initial conditions

uniform 0

define initial water

level of 0.5 m

using

boxToCell in

setFields (volScalarFieldValue alpha.water 1)

uniform 0

uniform (0 0 0)

Appendix B: Test case overview

107

B5 Two-dimensional modeling of groundwater flow in sediment under rippled

streambed with PCSiWaPro®

Table B 11: Model setup for test case B5 (surface water flow over ripples using PCSiWaPro®)

General

Referred to in section

References

Published in

2.3.5

Broecker et al. (2019) (for ripple geometry)

Sobhi Gollo et al. (2022a)

Domain discretization

Dimensions

Mesh generator

Number of cells

two-dimensional (length: 7 m, height: 0.556 m)

PCSiWaPro mesh generator

15524

Turbulence model

Hydrodynamic simulations

Solver

Time step

Simulation time

PCSiWaPro® (Richards equations)

min: 0.0001 s, max: 1 s

steady state

Transport simulations

Not performed

Table B 12: Boundary conditions for test case B5

left

Linear pressure head distribution to account for constant water level of 0.495 m over the top

node

right

Linear pressure head distribution to account for constant water level of 0.487 m over the top

node

top

Pressure head distribution from bottom wall modeling results of case B4 – applied node to

node

bottom

No flow

initial

conditions

Domain fully saturated (water content = 0.25), material: DIN 4220 pure sand, porosity is

manually set to 0.25, hydraulic conductivity is set to 4.64 × 10-3 m/s

Appendix B: Test case overview

108

B6 Three-dimensional two-phase integral modeling of flow across rippled

streambed with porousInter

Table B 13: Model setup for test case B6 (Integral modeling of groundwater-surface water

interactions across rippled streambeds using porousInter)

General

Referred to in section

References

Published in

2.3.5

Broecker et al. (2019)

Broecker et al. (2019), Sobhi Gollo et al. (2022a)

Domain discretization

Dimensions

Mesh generator

Number of cells

three-dimensional (length: 15 m, height: 1.5 m (max),

depth: 1 m)

gmsh

1806450

Turbulence model

LES

Hydrodynamic simulations

Solver

Time step

Simulation time

porousInter

variable

300 s

Transport simulations

Not performed

Appendix B: Test case overview

109

Table B 14: Boundary conditions for test case B6

alpha.water [-]

p_rgh [kg/(ms²)]

U [m/s]

effPackingRadius

[m]

Porosity [-]

inlet_water

inletOutlet

inletValue: uniform

value: uniform

fixedFluxPressure

value: uniform

flowRateInletVelocity

volumetricFlowRate 0.

5 m3/s

value

: uniform (1 0 0)/(0.5 0 0)

zeroGradient

inlet_air

inletOutlet

inletValue:

uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

outlet

zeroGradient

totalPressure

p0: uniform 0

value: uniform 0

zeroGradient

atmosphere

inletOutlet

inletValue:

uniform 0

value: uniform 0

totalPressure

p0: uniform 0

value: uniform 0

pressureInletOutletVelocity

value

uniform (0 0 0)

zeroGradient

walls

zeroGradient

fixedFluxPressure

value: uniform

fixedValue

value uniform (0 0 0)

zeroGradient

side

slip

initial

conditions

uniform 0

define initial water

level

from x = 0 m to x = 12 m and

y =

0.5 m to 0.5 m and z = 0 m to z = 1 m

using

boxToCell in

setFields (volScalarFieldValue

alpha.water 1)

uniform 0

uniform (0 0 0)

uniform 0.002

uniform 1

in the sediment

0.25 defined by

setFields

Appendix B: Test case overview

110

B7 Three-dimensional flow simulation in and around a ventilated U-shaped

burrow using porousInter – case 1

Table B 15: Model setup for test case B7 (Integral modeling of U-shaped burrow ventilation

using porousInter - case 1)

General

Referred to in section

References

Published in

3.3.1

Brand et al. (2013)

Sobhi Gollo et al. (2021)

Domain discretization

Dimensions

Mesh generator

Number of cells

three-dimensional (half-cylinder: height: 0.25 m, radius: 0.25 m)

Salome 9.3.0

450399

Turbulence model

laminar

Hydrodynamic simulations

Solver

Time step

Simulation time

porousInter

variable

100 s

Transport simulations

Not performed

Appendix B: Test case overview

111

Table B 16: Boundary conditions for test case B7

alpha.water [-]

p_rgh [kg/(ms²)]

U [m/s]

effPackingRadius

[m]

porosity [-]

top

zeroGradient

totalPressure

p0: uniform 0

value: uniform 0

pressureInletOutletVelocity

value

uniform (0 0 0)

zeroGradient

symmetry

symmetryPlane

Walls

(including

bottom)

zeroGradient

fixedFluxPressure

gradient

: uniform 0

value: uniform 0

fixedValue

value: uniform (0 0 0)

zeroGradient

baffle

zeroGradient

fixedFluxPressure

gradient

: uniform 0

value: uniform 0

fixedValue

value: uniform (0.005 0 0) up

to (0.025 0 0)

zeroGradient

initial

conditions

uniform 1

uniform 0

uniform (0 0 0)

uniform 0.000125 to

0.00025

Defined 1 for U-shaped

burrow and

in the sediment 0.39 using

SurfaceToCell in setFields

Appendix B: Test case overview

112

B8 Three-dimensional flow simulation in and around a ventilated U-shaped

burrow using porousInter – case 2

Table B 17: Model setup for test case B8 (Integral modeling of U-shaped burrow ventilation

using porousInter - case 2)

General

Referred to in section

References

Published in

3.4.2

Morad et al. (2010), Roskosch et al. (2010)

Sobhi Gollo et al. (2021)

Domain discretization

Dimensions

Mesh generator

Number of cells

three-dimensional (half-cylinder: height: 1.1 m, radius: 1 m)

Salome 9.3.0

399531

Turbulence model

laminar

Hydrodynamic simulations

Solver

Time step

Simulation time

porousInter

variable

40 s

Transport simulations

Not performed

Appendix B: Test case overview

113

Table B 18: Boundary conditions for test case B8

alpha.water

[-]

p_rgh

[kg/(ms²)] U [m/s]

effPackingRadius

[m]

porosity

[-]

top

slip

totalPressure

p0: uniform 0

value: uniform 0

pressureInletOutletVelocity

value

uniform (0 0 0)

zeroGradient

symmetry

symmetryPlane

Walls

(including

bottom)

zeroGradient

fixedFluxPressure

gradient

: uniform 0

value: uniform 0

fixedValue

value: uniform (0 0 0)

zeroGradient

baffle

fixedValue

value: uniform 1

fixedFluxPressure

gradient

: uniform 0

value: uniform 0

fixedValue

value: uniform (0.04 0 0)

zeroGradient

initial

conditions

uniform

(0 0 0)

uniform

0.00025

Defined 1 for U

-shaped

burrow and overlying water,

in the sediment 0.39 using

SurfaceToCell in setFields

Appendix B: Test case overview

114

B9 One-dimensional tracer transport in heterogenous groundwater with

porousInterTracer

Table B 19: Model setup for test case B9 (One-dimensional tracer transport in groundwater

with two porosities)

General

Referred to in section

References

Published in

4.4.1

Kinzelbach (1992)

Sobhi et al. (2022b)

Domain discretization

Dimensions

Mesh generator

Number of cells

one-dimensional (length: 10 m, width: 1 m)

blockMesh

10000

Turbulence model

laminar

Hydrodynamic simulations

Not performed

Transport simulations

Solver

Time step

Simulation time

Tracer diffusivity

porousInterTracer

variable

700 s

-9

Appendix B: Test case overview

115

Table B 20: Boundary conditions for test case B9.

alpha.water (-)

p_rgh (kg/(ms²))

U (m/s)

C (kg/m³)

effPackingRadius (m)

porosity (-)

inlet

inletOutlet

inletValue:

uniform

value: uniform

zeroGradient

fixedValue

value

uniform

(0.01 0 0)

fixedValue

value: uniform

zeroGradient

outlet

zeroGradient

fixedValue

value: uniform 0

zeroGradient

Upper,

lower and

sidewalls

empty

zeroGradient

initial

conditions

uniform 1

uniform 0

uniform

(0.01 0 0)

uniform 0

uniform 0.01

Define sediment porosity 0.44 from x = 0

m to x= 5 m and 0.33 for x = 5 m to x =

10 m through y = 0 m to 1 m and z = 0 m

to z = 1 m

using boxToCell in setFields

(volScalarFieldValue)

Appendix B: Test case overview

116

B10 Transport simulations for groundwater-surface water interactions at

heterogenous rippled streambeds

Table B 21: Model setup for test case B10 (Transport simulations for groundwater-surface

water interactions at heterogenous rippled streambeds).

General

Referred to in section

References

Published in

4.3.1, 4.3.2

Fox et al. (2016)

Sobhi Gollo et al. (2022b)

Domain discretization

Dimensions

Mesh generator

Number of cells

three-dimensional (length: 2.375 m,

height: 0.29 m, width: 0.29 m)

gmsh

hexahedra: 17110, prisms: 55014

Turbulence model

LES

Hydrodynamic simulations

Solver

Time step

Simulation time

porousInterTracer

variable

300 + 3600 s

Transport simulations

Solver

Time step

Simulation time

Tracer diffusivity

porousInterTracer

variable

3600 s

-9

Appendix B: Test case overview

117

Table B 22: Boundary conditions for test case B10.

alpha.water

(-)

p_rgh (kg/(ms²))

U (m/s)

C (kg/m³)

effPackingRadius

(m)

porosity (-)

inlet

fixedValue

value:

uniform

zeroGradient

flowRateInletVelocity

volumetricFlowRate

003045

value

: uniform

(

0.0362069 0 0)

inletOutlet

inletValue:

uniform

value:

uniform

zeroGradient

outlet

inletOutlet

inletValue:

uniform

value:

uniform

fixedValue

value: uniform 0

zeroGradient

atmosphere

slip

zeroGradient

lower walls

zeroGradient

fixedValue

value:

uniform (0 0 0)

zeroGradient

sidewalls

slip

Appendix B: Test case overview

118

groundwater

inflow (for

gaining and

losing

conditions)

fixedValue

value:

uniform

fixedFluxPressure

flowRateInletVelocity

volumetricFlowRate

± 4.85

× 10-6 m3/s

value

: uniform (0

± 1.15738

× 10-5

inletOutlet

inletValue:

uniform

value:

uniform

zeroGradient

initial

conditions

uniform 1

uniform 0

for surface water:

uniform (0.15 0 0), for

the sediment:

uniform

(0 0 0) for neutral

conditions, uniform

(

0 ± 1.15738

× 10-5

0) for gaining and

losing conditions

defined by setFields

uniform 0

nonuniform, in

the sediment, a

combination of

0.00038, 0.0013

and 0.0023

defined by

setFields

(surfaceTo Cell)

uniform 1,

in the

sediment, a

combination

of 0.33 and

0.44 uniform

0.33 defined

by setFields

(surfaceTo

Cell)

Bibliography

119

Bibliography

Aller, R. C. 1980. Quantifying Solute Distributions in the Bioturbated Zone of Marine

Sediments by Defining an Average Microenvironment, Geochimica et Cosmochimica Acta,

44(12), 1955–1965.

Arnon, S., K. Yanuka, and A. Nejidat. 2013. Impact of overlying water velocity on ammonium

uptake by benthic biofilms. Hydrological Processes 27, no. 4: 570-578.

Arnon, S., L. P. Marx, K. E. Searcy, and A. I. Packman. 2010. Effects of overlying velocity,

particle size, and biofilm growth on stream–subsurface exchange of particles. Hydrological

processes. 24, 108.

Ashby, S.F. and R.D. Falgout. 1996. A Parallel Multigrid Preconditioned Conjugate Gradient

Algorithm for Groundwater Flow Simulations. Nuclear Science and Engineering 124, 145–

159. https://doi.org/10.13182/NSE96-A24230.

Bardini, L., F. Boano, M. B. Cardenas, R. Revelli, and L. Ridolfi. 2012. Nutrient cycling in

bedform induced hyporheic zones. Geochimica et Cosmochimica Acta, 84, 47-61,

doi:10.1016/j.gca.2012.01.025.

Bartram, J., R. Balance. 1996. United Nations, and World Health Organization, eds. Water

Quality Monitoring: A Practical Guide to the Design and Implementation of Freshwater Quality

Studies and Monitoring Programmes. 1st ed. London ; New York: E & FN Spon chapter 12.

ISBN:0419217304.

Bayon, A., D. Valero, R. García-Bartual, F. J. Vallés-Morán, and P. A. López-Jiménez. 2016.

Performance assessment of OpenFOAM and FLOW-3D in the numerical modeling of a low

Reynolds number hydraulic jump. Environmental Modelling & Software 80: 322–35,

doi:10.1016/j.envsoft.2016.02.018.

Bencala, K. E. 2000. Hyporheic zone hydrological processes. Hydrological Processes, 14(15),

2797-2798, doi:10.1002/1099-1085(20001030)14:15<2797::Aid-hyp402>3.0.Co;2-6.

Bijeljic, B., and M. J. Blunt. 2006. Pore-scale modeling and continuous time random walk

analysis of dispersion in porous media. Water Resources Research 42 (1):5,

doi:10.1029/2005wr004578.

Biksey, T. M., and E. D. Gross. 2001. The hyporheic zone: linking groundwater and surface

water - Understanding the paradigm. Remediat. J. 12, 55.

Bibliography

120

Blois, G., J. L. Best, G. H. Sambrook Smith, and R. J. Hardy. 2014. Effect of bed permeability

and hyporheic flow on turbulent flow over bed forms. Geophysical Research Letters, 41(18),

6435-6442, doi:10.1002/2014gl060906.

Boano, F., C. Camporeale, R. Revelli, L. Ridolfi. 2006. Sinuosity‐driven hyporheic exchange

in meandering rivers. Geophysical Research Letters, 33. doi:10.1029/2006GL027630.

Boano, F., C. Camporeale, R. Revelli. 2010. A linear model for the coupled surface‐subsurface

flow in a meandering stream. Water Resources Research, 46, W07535,

https://doi.org/10.1029/2009WR008317.

Boano, F., J. W. Harvey, A. Marion, A. I. Packman, R. Revelli, L. Ridolfi, and A. Wörman.

2014. Hyporheic flow and transport processes: Mechanisms, models, and biogeochemical

implications, Rev. Geophysics, 52, 603–679, https://doi.org/10.1002/2012RG000417.

Bobba, A. Ghosh. 2012. Ground Water-Surface Water Interface (GWSWI) Modeling: Recent

Advances and Future Challenges. Water Resources Management 26, Nr. 14 4105–31, doi:

10.1007/s11269-012-0134-x.

Bottacin‐Busolin, A., and A. Marion. 2010. Combined role of advective pumping and

mechanical idispersion on tme scales of bed form— Induced hyporheic exchange. Water

Resources Research, 46, W08518, https://doi.org/10.1029/2009WR008892.

Boulton, A.J., S. Findlay, P. Marmonier, E.H. Stanley, and H.M. Valett. 1998. The functional

significance of the hyporheic zone in streams and rivers. Ecological systems. 29, 59.

Boulton, AJ., T. Datry, T. Kashara, M. Mutz, and J. A. Stanford. 2010. Ecology and

management of the hyporheic zone: stream groundwater interactions of running waters and

their floodplains. Journal of the North American Benthological Society, vol. 29, no. 1, p. 26-

40, http://dx.doi.org/10.1899/08-017.1.

Boutron, O., C. Margoum, J.-M., Chovelon, , C. Guillemain, and V. Gouy. 2011. Effect of the

submergence, the bed form geometry, and the speed of the surface water flow on the mitigation

of pesticides in agricultural ditches. Water Resources Research, 47(8),

doi:10.1029/2011wr010378.

Brand, A., J. Lewandowski, E. Hamann, and G. Nützman. 2013. Advection around ventilated

U-Shaped burrows: A model study: Advection around ventilated U-shaped burrows. Water

Resources Research 49, no. 5 (May 2013): 2907–17.

Bibliography

121

Broecker, T., K. Teuber, V. Sobhi Gollo, G. Nützmann, J. Lewandowski, and R. Hinkelmann.

2019. Integral flow modelling approach for surface water-groundwater interactions along a

rippled streambed. Water 11 no. 7: 1517, doi:10.3390/w11071517.

Broecker, T., V. Sobhi Gollo, A. Fox, J. Lewandowski, G. Nützmann, S. Arnon, and R.

Hinkelmann. 2021. High‐resolution integrated transport model for studying surface water–

groundwater interaction. Groundwater 59, No:4 488-502, doi:10.1111/gwat.13071.

Brunke, M. 2001. Wechselwirkungen zwischen Fließgewässer und Grundwasser: Bedeutung

für aquatische Biodiversität, Stoffhaushalt und Lebensraumstrukturen. Wasserwirtschaft, 90,

32-37.

Brunke, M., and T. Gonser. 1997. The ecological significance of exchange processes between

rivers and groundwater. Freshwater Biology, 37(1), 1-33, doi:10.1046/j.1365-

2427.1997.00143.x.

Brunner, P., R. Therrien, P. Renard, C. T. Simmons and H.‐J. H. Franssen. 2017. Advances in

understanding river–groundwater interactions. Reviews of Geophysics, 55, 818–854,

https://doi.org/10.1002/2017RG000556.

Buffington, J.M., and D. Tonina. 2009. Hyporheic exchange in mountain rivers II: Effects of

channel morphology on mechanics, scales, and rates of exchange. Geography Compass 3, no.

3: 1038-1062.

Cardenas, M. B. 2009. Stream-aquifer interactions and hyporheic exchange in gaining and

losing sinuous streams. Water Resources Research, 45(6), doi:10.1029/2008wr007651.

Cardenas, M. B., and J. L. Wilson. 2007a. Dunes, turbulent eddies, and interfacial exchange

with permeable sediments. Water Resources Research, 43(8), W08412,

doi:10.1029/2006wr005787.

Cardenas, M. B., and J. L. Wilson. 2007b. Effects of current‐bed form induced fluid flow on

thermal regime of sediments. Water Resources Research, 43, W08431,

https://doi.org/10.1029/2006WR005343.

Cardenas, M.B. 2009. Stream-aquifer interactions and hyporheic exchange in gaining and

losing sinuous streams. Water Resources Research 45, no. 6: W06429.

Cardenas, M.B., J.L. Wilson, and V.A. Zlotnik. 2004. Impact of heterogeneity, bed forms, and

stream curvature on subchannel hyporheic exchange. Water Resources Research 40, no. 8:

W08307.

Bibliography

122

Caretto, L. S., A. D. Gosman, S. V. Patankar, D. B. Spalding. 1973. Two calculation procedures

for steady, three-dimensional flows with recirculation. Paper presented at the Proceedings of

the Third International Conference on Numerical Methods in Fluid Mechanics, Berlin,

Heidelberg.

Chapman, S. W., B. L. Parker, J. A. Cherry, R. Aravena and D. Hunkeler. 2007. Groundwater‐

surface water interaction and its role on TCE groundwater plume attenuation. Journal of

Contaminant Hydrology, 91(3‐4), 203–232. https://doi.org/10.1016/j.jconhyd.2006.10.006.

Chen, X., M. B. Cardenas, and L. Chen. 2015. Three-dimensional versus two-dimensional bed

form-induced hyporheic exchange. Water Resources Research, 51(4), 2923-2936,

doi:10.1002/2014wr016848.

Chen, X., M. B. Cardenas, L. Chen. 2018. Hyporheic Exchange Driven by Three-Dimensional

Sandy Bed Forms: Sensitivity to and Prediction from Bed Form Geometry. Water Resources

Research, 54(6), 4131-4149, doi:10.1029/2018wr022663.

Coleman, J. C., M. C. Miller, and F. L. Mink. 2011. Hydrologic disturbance reduces biological

integrity in urban streams. Environmental monitoring. 172, 663.

Coluccio, K. and L. K. Morgan. 2019. A review of methods for measuring groundwater–surface

water exchange in braided rivers. Hydrology and Earth System Sciences 23(10):4397-4417,

doi:10.5194/hess-23-4397-2019.

Conant, B., J. Cherry, and R. Gillham. 2004. A PCE groundwater plume discharging to a river:

Influence of the streambed and near‐river zone on contaminant distributions. Journal of

Contaminant Hydrology, 73(1‐4), 249–279, doi:10.1016/j.jconhyd.2004.04.001.

Cuthbert, M., R. Mackay, V. Durand, M-F. Aller, R.B. Greswell, M. O. Rivett. 2010. Impacts

of river-bed gas on the hydraulic and thermal dynamics of the hyporheic zone. Advances in

Water resources management 33.

Darwin, C. 1881. The formation of vegetable mould through the action of worms with

observation of their habits. John Murray.

Derx, J., A.P. Blaschke, G. Blöschl. 2010. Three dimensional flow patterns at the river-aquifer

interface—a case study at the Danube. Advances in Water resources management 33.

DIN 4220:2008-11. 2008. Pedologic Site Assessment - Designation, Classification and

Deduction of Soil Parameters (Normative and Nominal Scaling), doi:10.31030/1436635.

Bibliography

123

Doyle, M. W. 2005. Incorporating hydrologic variability into nutrient spiraling. Journal of

Geophysical Research: Biogeosciences, 110(G1), doi:10.1029/2005jg000015.

Edwards, R. T. 1998. The Hyporheic zone. In R. J. Naiman, & R. E. Bilby (Eds.), River ecology

and management: Lessons from the Pacific coastal ecoregion, (pp. 399–429). New York:

Springer, https://doi.org/10.1007/978‐1‐4612‐1652‐0_16.

Edwards, R.T. 1998. The hyporheic zone, in River Ecology and Management: Lessons from

the Pacific Coastal Ecoregion. edited by R. J. Naiman and R. E., Bilby, pp. 399–429, Springer,

New York.

Elliott, A. H., and N. H. Brooks. 1997a. Transfer of nonsorbing solutes to a streambed with

bed forms: Laboratory experiments. Water Resources Research, 33(1), 137-151,

doi:10.1029/96wr02783.

Elliott, A. H., and N. H. Brooks. 1997b. Transfer of nonsorbing solutes to a streambed with

bed forms: Theory. Water Resources Research, 33(1), 123-136, doi:10.1029/96wr02784.

Ellis, P. A., R. Mackay, and M. O. Rivett. 2007. Quantifying urban river aquifer fluid exchange

processes: A multi‐scale problem. Journal of Contaminant Hydrology, 91(1‐2), 58–80,

https://doi.org/10.1016/j.jconhyd.2006.08.014.

Ergun, S. 1952. Fluid Flow Through Packed Column. Chemical Engineering Progress 48: 89-

94.

Feris, K., P. Ramsey, C. Frazar, J. N. Moore, J. E. Gannon, and W. E. Holben. 2003a.

Differences in Hyporheic-Zone Microbial Community Structure along a Heavy-Metal

Contamination Gradient. Applied and Environmental Microbiology, 69(9), 5563-5573,

doi:10.1128/aem.69.9.5563-5573.2003.

Feris, K., P. Ramsey, C. Frazar, M. C. Rillig, J. E. Gannon, W. E. Holben. 2003b. Structure and

seasonal dynamics of hyporheic zone microbial communities in free-stone rivers of the estern

United States. Microbial Ecology, 46(2), 200-215, doi:10.1007/BF03036883.

Fischer, H., F. Kloep, S. Wilzcek, and M.T. Pusch. 2005. A river’s liver—Microbial processes

within the hyporheic zone of a large lowland river. Biogeochemistry 76, 349.

Fleckenstein, J. H., S. Krause, , D. M. Hannah, and F. Boano. 2010. Groundwater–surface

water interactions: New methods and models to improve understanding of processes and

dynamics, Advances Water Resources, 33, 1291–1295, doi:10.1016/j.advwatres.2010.09.011.

Bibliography

124

Fox, A., F. Boano, S. Arnon. 2014. Impact of losing and gaining streamflow conditions on

hyporheic exchange fluxes induced by dune-shaped bed forms. Water Resources Research 50,

1895–1907. doi:10.1002/2013WR014668.

Fox, A., G. Laube, C. Schmidt, J.H. Fleckenstein, and S. Arnon. 2016. The effect of losing and

gaining flow conditions on hyporheic exchange in heterogeneous streambeds: Hyporheic

exchange in heterogenous streambeds. Water Resources Research 52, No. 9: 7460–77, doi:

10.1002/2016WR018677

Fraser, B.G., and D.D. Williams. 1998. Seasonal boundary dynamics of a groundwater/surface

water ecotone. Ecology 79(6):2019-2031, doi:10.2307/176706.

Freeze, R. A., and J. A. Cherry. 1979. Groundwater 604 pp., Prentice Hall, Englewood Cliffs,

N. J.

Furukawa, Y., S. J. Bentley, and D. L. Lavoie. 2001. Bioirrigation Modeling in Experimental

Benthic Mesocosms, Journal of Maritime Research., 59(3), 417–452.

Gent, M. 1995. Wave Interaction with Permeable Coastal Structures, Elsevier Science:

Amsterdam, The Netherlands, Volume 95. ISBN:90-407-1182-8.

Geuzaine, C., and J.-F. Remacle. 2009. Gmsh: A 3-D finite element mesh generator with built-

in pre- and postprocessing facilities. International Journal for Numerical Methods in

Engineering 79 no. 11: 1309-1331.

Gibert, J., and L. Deharveng. 2002. Subterranean ecosystems: A truncated functional

biodiversity. Bioscience 52, 473.

Gomez, J. D., J. L. Wilson, M. B. Cardenas. 2012. Residence time distributions in sinuosity‐

driven hyporheic zones and their biogeochemical effects. Water Resources Research, 48,

W09533, https://doi.org/10.1029/2012WR012180.

Gomez-Velez, J.D., S. Krause, and J.L. Wilson. 2014. Effect of low-permeability layers on

spatial patterns of hyporheic exchange and groundwater upwelling. Water Resources Research

50, no. 6: 5196–5215.

González-Pinzón, R., A. S. Ward, C. E. Hatch, , A. N. Wlostowski, K. Singha, M.N. Gooseff,

R. Haggerty, J. W. Harvey, O. A. Cirpka, and J. T. Brock. 2015. A field comparison of multiple

techniques to quantify groundwater–surface-water interactions, Freshwater Science., 34, 139–

160, https://doi.org/10.1086/679738.

Bibliography

125

Gorelick, S. M., editor. 1986. Conjunctive water use: understanding and managing

groundwater-surface water interactions. International Association of Hydrologic Sciences,

Publication no. 156, Wallingford, 547pp. ISBN: 9780947571658.

Goyeau, B., D. Lhuillier, D. Gobin, and M. Velarde. 2003. Momentum transport at a fluid-

porous interface, International journal of Heat Mass Transfer, 46, 4071–4081,

doi:10.1016/S0017-9310(03)00241-2.

Greenshields, C. 2010. OpenFOAM 1.7.0 Released. Retrieved from

http://openfoam.org/release/1-7-0/

Guo, J., R. Blankenburg, X. Geng, and P-W. Graeber. 2017. Hydrological process analysis in

earth dams using the PCSiWaPro® as a basis for stability analysis. Journal of Geological

Resource and Engineering 5, no. 3, doi: 10.17265/2328-2193/2017.03.004.

Harvey, J. W., and B. J. Wagner. 2000. Quantifying hydrologic interactions. In: Jones J.B.,

Mulholland P.J. (eds) Streams and groundwaters. Academic press.

Harvey, J. W., and K. E. Bencala. 1993. The effect of streambed topography on surface‐

subsurface water exchange in mountain catchments. Water Resources Research, 29(1), 89–98,

https://doi.org/10.1029/92WR01960.

Harvey, J. W., J. K. Böhlke, M. A. Voytek, D. Scott, and C. R. Tobias. 2013. Hyporheic zone

denitrification: Controls on effective reaction depth and contribution to whole stream mass

balance. Water Resources Research, 49, 6298–6316, https://doi.org/10.1002/wrcr.20492.

Hassanizadeh, S. M., M. A. Celia, H. K. Dahle. 2002. Dynamic Effect in the Capillary Pressure-

Saturation Relationship and its Impacts on Unsaturated Flow. Vadose Zone Journal 1 (1):38-

57, doi:10.2113/1.1.38.

Hayashi, M., and D.O. Rosenberry. 2002. Effects of ground water exchange on the hydrology

and ecology of surface water. Ground Water 40(3):309-316, doi: 10.1111/j.1745-

6584.2002.tb02659.x.

Higuera, P., J. L. Lara, and I. J. Losada. 2014. Three-dimensional interaction of waves and

porous coastal structures using OpenFOAM®. Part II: Application. Coastal Engineering 83:

259–70, doi:10.1016/j.coastaleng.2013.08.010.

Bibliography

126

Hill, A. R., C. F. LaBadia, and K. Sanmugadas. 1998. Hyporheic zone hydrology and nitrogen

dynamics in relation to the streambed topography of a N-rich stream. Biogeochemistry 42:285-

310, doi:10.1023/A:1005932528748.

Hill, A., and M. Cardaci. 2004. Denitrification and organic carbon availability in riparian

wetland soils and subsurface sediments. Soil Science Society of America Journal 68, 320-325.

Hinkelmann R. 2005. Efficient Numerical Methods and Information-Processing Techniques for

Modeling Hydro- and Environmental Systems. Lecture Notes in Applied and Computational

Mechanics 21.

Hölker, F., M. J. Vanni, J. J. Kuiper, C. Meile, H-P. Grossart, P. Stief, R. Adrian, Lorke, A.

Dellwig, A. Brand, M. Hupfer, W. M. Moolj, G. Nützmann, and J. Lewandowski. 2015. Tube-

Dwelling Invertebrates: Tiny Ecosystem Engineers Have Large Effects in Lake Ecosystems,

Ecological Monographs, 85(3), Pp. 333–351.

Hupfer, M., S. Jordan, C. Herzog, C. Ebeling, R. Ladwig, M. Rothe, and J. Lewandowski. 2019.

Chironomid Larvae Enhance Phosphorus Burial in Lake Sediments: Insights from Long-Term

and Short-Term Experiments. Science of The Total Environment 663: 254–64,

https://doi.org/10.1016/j.scitotenv.2019.01.274.

Ibisch, R., S. I, and D. Borchardt. 2009. Influence of periphyton biomass dynamics on

biological colmation processes in the hyporheic zone of a gravel bed river (River Lahn,

Germany). Advances in Limnology, 61, 87-104.

Issa, R. I. 1986. Solution of the implicitly discretised fluid flow equations by operator-splitting.

Journal of Computational Physics, 62(1), 40-65, doi:10.1016/0021-9991(86)90099-9.

Janssen, F., M. B. Cardenas, A. H. Sawyer, T. Dammrich, J. Krietsch de Beer DA. 2012.

Comparative Experimental and Multiphysics Computational Fluid Dynamics Study of Coupled

Surface–Subsurface Flow in Bed Forms. Water resources reearch, 48,

doi:10.1029/2012WR011982.

Jeppesen, E., J. P. Jensen, M. Sondergaard, T. Lauridsen, , P. H. Moller, K. Sandby. 1998.

Changes in nitrogen retention in shallow eutrophic lakes following a decline in density of

cyprinids. Archiv für Hydrobiologie 142, 129‐151.

Jin, G., H. Tang, B. Gibbes, L. Li and D. Barry. 2010. Transport of nonsorbing solutes in a

streambed with periodic bedforms, Advances in Water Resources 33: 1402-1416, doi:

10.1016/j.advwatres.2010.09.003.

Bibliography

127

Jones, J. P., E. A. Sudicky, R. G. McLaren. 2008. Application of a fully-integrated surface-

subsurface flow model at the watershed-scale: a case study. Water resources research 44,

doi: 10.1029/2006WR005603.

Jones, J.B., and P.J. Mulholland. 2000. In Streams and Ground Waters. San Diego: Academic

Press.

Jones, J.B., and R.M. Holmes. 1996. Surface-subsurface interactions in stream ecosystems.

Trends in Ecology and Evolution 11:239-242, doi: 10.1016/0169-5347(96)10013-6.

Kalbus, E., C. Schmidt, M. Bayer‐Raich, S. Leschik, F. Reinstorf, G. Balcke, and M. Schirmer.

2007. New methodology to investigate potential contaminant mass fluxes at the stream‐aquifer

interface by combining integral pumping tests and streambed temperatures. Environmental

Pollution, 148(3), 808–816, https://doi.org/10.1016/j.envpol.2007.01.042.

Kalbus, E., F. Reinstorf, and M. Schirmer. 2006. Measuring methods for groundwater and

surface water interactions: Hydrol. Earth Syst. Sci., 10(6), 873-887, doi:10.5194/hess-10-873-

2006.

Kennedy, C. D., D. P. Genereux, D. R. Corbett, and H. Mitasova. 2009. Spatial and temporal

dynamics of coupled groundwater and nitrogen fluxes through a streambed in an agricultural

watershed. Water Resources Research, 45(9), doi:10.1029/2008wr007397.

Khan, H. H., and A. Khan. 2019. Chapter 14 - Groundwater and Surface Water Interaction. In

S. Venkatramanan, M. V. Prasanna, & S. Y. Chung (Eds.), GIS and Geostatistical Techniques

for Groundwater Science (pp. 197-207): Elsevier.

Kinzelbach, W. 1992. Numerische Methoden zur Modellierung des Transports von

Schadstoffen im Grundwasser: Oldenbourg Wissenschaftsverlag.

Klein, S., E. Worch, and T. P Knepper. 2015. Occurrence and Spatial Distribution of

Microplastics in River Shore Sediments of the Rhine-Main Area in Germany. Environmental

Science & Technology, 49(10), 6070-6076, doi:10.1021/acs.est.5b00492.

Kolditz, O., U.-J. Görke, H. Shao, W. Wang. (Eds.). 2012. Thermo-Hydro-Mechanical-

Chemical Processes in Porous Media: Benchmarks and Examples, Lecture Notes in

Computational Science and Engineering. Springer Berlin Heidelberg, Berlin, Heidelberg.

https://doi.org/10.1007/978-3-642-27177-9.

Kollet, S. J., R. M. Maxwell. 2006. Integrated surface-groundwater flow modeling: a

freesurface overland flow boundary condition in a parallel groundwater flow model. Advances

in Water resources management 29:945–58, doi:10.1016/j.advwatres.2005.08.006.

Bibliography

128

Koretsky, C. M., C. Meile, and P. Van Cappellen. 2002. Quantifying Bioirrigation Using

Ecological Parameters: A Stochastic Approach’, Geochemical Transactions, Vol. 3, No. 3, Pp

17-30.

Krause, S., A. Bronstert, and E. Zehe. 2007. Groundwater‐surface water interactions in a North

German lowland floodplain—Implications for the river discharge dynamics and riparian water

balance. Journal of Hydrology, 347(3‐4), 404–417, https://doi.org/10.1016/j.

jhydrol.2007.09.028.

Krause, S., D. M. Hannah, J. H. Fleckenstein, C. M. Heppell, D. Kaeser, R. Pickup, G. Pinay,

A. L. Robertson, and P. J. Wood. 2011. Inter-disciplinary perspectives on processes in the

hyporheic zone. Ecohydrology, 4(4), 481-499, doi:10.1002/eco.176.

Kristensen, E., G. Penha-Lopes, M. Delefosse, T. Valdemarsen, C. O. Quintana, and G. T.

Banta. 2012. What is bioturbation? The need for a precise definition for fauna in aquatic

sciences. Marine Ecology Progress series. 446:285-302.

Kunkel, U., and M. Radke. 2008. Biodegradation of Acidic Pharmaceuticals in Bed Sediments:

Insight from a Laboratory Experiment. Environmental Science & Technology 42, no. 19 7273–

79.

Larkin, R.G., and J.M. Sharp. 1992. On the relationship between river basin geomorphology,

aquifer hydraulics, and groundwater flow direction in alluvial aquifers. Geol. Soc. Am. Bull.

104, 1608–1620, doi:10.1130/0016-7606(1992)104<1608:OTRBRB>2.3.CO;2.

Laskov, C, C. Herzog, J. Lewandowski, and M. Hupfer. 2007. Miniaturised Photometrical

Methods for the Rapid Analysis of Phosphate, Ammonium, Ferrous Iron and Sulfate in Pore

Water of Aquatic Sediments. Limnology and Oceanography: Methods.

Lawrence, J., M. Skold, F. Hussain, D. Silverman, V. Resh, D. Sedlak, R. Luthy, and J. McCray.

2013. Hyporheic Zone in Urban Streams: A Review and Opportunities for Enhancing Water

Quality and Improving Aquatic Habitat by Active Management. Environmental Engineering

Science, 30, 480-501, doi:10.1089/ees.2012.0235.

Lewandowski, J., C. Laskov, and M. Hupfer. 2007. The Relationship between Chironomus

Plumosus Burrows and the Spatial Distribution of Pore-Water Phosphate, Iron and Ammonium

in Lake Sediments. Freshwater Biology 52, 331–43.

Bibliography

129

Lewandowski, J., K. Meinikmann, and S. Krause. 2020. Groundwater–Surface Water

Interactions: Recent Advances and Interdisciplinary Challenges. Water, 12(1), 296.

doi:10.3390/w12010296.

Lewandowski, J., K. Rüter, and M. Hupfer. 2002. Two-Dimensional Small-Scale Variability of

Pore Water Phosphate in Freshwater Lakes: Results from a Novel Dialysis Sampler.

Environmental Science Technology, 36, 2039-2047.

Lewandowski, J., M. Hupfer. 2005a. Effect of Macrozoobenthos on Two-Dimensional Small-

Scale Heterogeneity of Pore Water Phosphorus Concentrations in Lake Sediments: A

Laboratory Study. Limnology and Oceanography 50:1106–1118.

Lewandowski, J., M. Hupfer. 2005b. Effect of Chironomus Plumosus on Spatial Distribution

of Porewater Phosphate Concentration in Lake Sediments: A Laboratory Experiment. Verh.

Internat. Verein. Theoret. Angew. Limnol. 29: 937-940.

Lewandowski, J., M. Schadach, and M. Hupfer. 2005c. Impact of Macrozoobenthos on Two-

Dimensional Small-Scale Heterogeneity of Pore Water Phosphorus Concentrations: In-Situ

Study in Lake Arendsee (Germany). Hydrobiologia 549, no. 1, 43–55.

Lewandowski, J., S. Arnon, E. Banks, O. Batelaan, A. Betterle, T. Broecker, C. Coll, J. D.

Drummond, J. Gaona Garcia, J. Galloway, J. Gomez-Velez, R. C. Grabowski, S. P. Herzog, R.

Hinkelmann, A. Höhne, J. Hollender, M. A. Horn, A. Jaeger, S. Krause, A. Löchner Prats, C.

Magliozzi, K. Meinikmann, B. B. Mojarrad, B. M. Mueller, I. Peralta-Maraver, A. L. Popp, M.

Posselt, A. Putschew, M. Radke, M. Raza, J. Riml, A. Robertson, C. Rutere, , J. L. Schaper, M.

Schirmer, H. Schulz, M. Shanafield, T. Singh, A. S. Ward, P. Wolke, A. Wörman, and L. Wu.

2019. Is the Hyporheic Zone Relevant beyond the Scientific Community? Water 11, no. 11:

2230.

Li, T., P. Troch, J. D. and Rouck. 2004. Wave overtopping over a sea dike, Journal of

Computational Physics 198 : 686-726, doi: 10.1016/j.jcp.2004.01.022.

Lopez-Rodriguez, M.J., J.M. Tierno de Figueroa, S. Fenoglio, T. Bo, and J. Alba-Tercedor.

2009. Life strategies of 3 Perlodidae species (Plecoptera) in a Mediterranean seasonal stream

in southern Europe. Journal of American Benthology. 28, 611.

Luckner, L., M. Th. Van Genuchten, and D. R. Nielsen. 1989. A Consistent set of parametric

models for the two-phase flow of immiscible fluids in the subsurface. Water Resources

Research 25, No. 10 2187–93, doi:10.1029/WR025i010p02187.

Bibliography

130

Magliozzi, C., R. C. Grabowski, A. I. Packman, and S. Krause. 2018. Toward a conceptual

framework of hyporheic exchange across spatial scales, Hydrol. Earth Syst. Sci., 22, 6163–

6185, https://doi.org/10.5194/hess-22-6163-2018.

Marzadri, A., D. Tonina, A. Bellin, G. Vignoli, and M.Tubino. 2010. Effects of bar topography

on hyporheic flow in gravel‐bed rivers. Water Resources Research, 46, W07531,

https://doi.org/10.1029/2009WR008285.

Marzadri, A., D. Tonina, A. Bellin. 2011. A semianalytical three‐dimensional process‐based

model for hyporheic nitrogen dynamics in gravel bed rivers. Water Resources Research, 47,

W11518, https://doi.org/10.1029/2011WR010583.

Mayer, K. U., E. O. Frind, and D. W. Blowes. 2002. Multicomponent reactive transport

modeling in variably saturated porous media using a generalized formulation for kinetically

controlled reactions, Water Resources Research 38, no. 9: 1174.

McKnight, D. M., E. W. Boyer, P. K. Westerhoff, P. T. Doran, T. Kulbe, and D. T. Andersen.

2001. Spectrofluorometric characterization of dissolved organic matter for indication of

precursor organic material and aromaticity. Limnology and Oceanography 46, no. 1: 38–48.

Medina, M. A., R. L. Doneker, N. Grosso, D. M. Johns, W. Lung, M. F. N. Mohsen, A. I.

Packman, and P. J. Roberts. 2002. Surface water–ground water interactions and modeling

applications, in Environmental Modeling and Management: Theory, Practice and Future

Directions, edited by, C. C. Chien et al., pp. 1–62, DuPont Company, Wilmington, Del.

Meile, C., P. Berg, P. Van Cappellen, and K. Tuncay. 2005. Solute-Specific Pore Water

Irrigation: Implications for Chemical Cycling in Early Diagenesis, Journal of Maritime

Research, 63(3), 601–621.

Merill, L., and D.J. Tonjes. 2014. A Review of the Hyporheic Zone, Stream Restoration, and

Means to Enhance Denitrification. Technology & Society Faculty Publications. 5.

Meysman, F. J. R., O. S. Galaktinov, B. Gribsholt, and J. J. Middelburg. 2006. Bioirrigation in

Permeable Sediments; Advective Pore-Water Transport Induced by Burrow Ventilation,

Limnology and Oceanography, Vol. 51, No. 1, Pp 142-156.

Mojarrad, B. B., J. Riml, A. Wörman, and H. Laudon. 2019. Fragmentation of the hyporheic

zone due to regional groundwater circulation. Water Resources Research, 55,

https://doi.org/10.1029/2018WR024609.

Bibliography

131

Morad, M. R., A. Khalili, A. Roskosch, and J. Lewandowski. 2010. Quantification of Pumping

Rate of Chironomus Plumosus Larvae in Natural Burrows. Aquatic Ecology 44, no. 1. 143–53.

Mosthaf, K., K. Baber, B. Flemisch, R. Helmig, A. Leijnse, I. Rybak, and B. Wohlmuth. (2011).

A coupling concept for two-phase compositional porous-medium and single-phase

compositional free flow. Water Resources Research, 47(10), doi:10.1029/2011wr010685.

Moulton, J. D., M. Berndt, R. Garimella, L. Prichett-Sheats, G. Hammond, M. Day, and J. C. Meza.

2012. High-level design of Amanzi, themulti-process high performance computing simulator, ASCEM-

HPC-2011-03-1, Off. of Environ. Manage., U. S. Dep. of Energy, Washing-ton, D. C.

Mulholland, P. J., E. R. Marzolf, J. R. Webster, and D. R. Hart. 1997. Evidence that hyporheic

zones increase heterotrophic metabolismand phosphorus uptake in forest streams. Limnology

and Oceanography, 44(1), 230–231, https://doi.org/10.4319/lo.1999.44.1.0230.

Müller, E. N. 2014. Skalen, Schwerpunkte, Rückkopplungen and Herausforderungen der

ökohydrologischen Forchung in Deutschland. Hydrologie und Wasserwirtschaftung 58.

Jahrgang, Heft 4.

Nogaro, G., T. Datry, F. Mermillod-Blondin, S. Descloux, and B. Montuelle. 2010. Influence

of streambed sediment clogging on microbial processes in the hyporheic zone. Freshwater

Biology. 55, 1288.

Nützmann, G., and S. Mey. 2007. Model-based estimation of runoff changes in a small lowland

watershed of north-eastern Germany. Journal of Hydrology, 334(3), 467-476,

doi:10.1016/j.jhydrol.2006.10.026.

Oxtoby, O., J. Heyns, and R. Suliman. 2013. A finite-volume solver for two-fluid flow in

heterogeneous porous media based on OpenFOAM. In Open Source CFD International

Conference. Hamburg, Germany, doi: 10.13140/2.1.3075.8400.

Packman, A. I., and N. H. Brooks. 2001. Hyporheic exchange of solutes and colloids with

moving bed forms. Water resources research 37:2591–605, doi:10.1029/2001WR000477.

Packman, A. I., M. Salehin, and M. Zaramella. 2004. Hyporheic Exchange with Gravel Beds:

Basic Hydrodynamic Interactions and Bedform-Induced Advective Flows. Journal of

Hydraulic Engineering, 130(7), 647-656, doi:10.1061/(ASCE)0733-9429(2004)130:7(647).

Paul, M.J., and J.L. Meyer. 2001. Streams in the urban landscape. Ecological systems. 32, 333.

Bibliography

132

Peter, K. T., S. Herzog, Z. Tian, C. Wu, J. E. McCray, K. Lynch, and E. P. Kolodziej. 2019.

Evaluating emerging organic contaminant removal in an engineered hyporheic zone using high

resolution mass spectrometry. Water Research, 150, 140-152,

doi:10.1016/j.watres.2018.11.050.

Powell, D. M. 1998. Patterns and processes of sediment sorting in gravel-bed rivers, Progress

in Physical Geography 22, no. 1: 1–32.

Prudic, D. 1989. Documentation of a computer program to simulate stream-aquifer relations

using a modular, finite-difference, ground-water flow model. In: UGS (USGS), editor. USGS

Open-File Report, Carson City, doi:10.3133/ofr88729.

Pryshlak, T.T., A.H. Sawyer, S.H. Stonedahl, and M.R. Soltanian. 2015. Multiscale hyporheic

exchange through strongly heterogeneous sediments. Water Resources Research 51, no. 11:

9127–9140.

Revelli, R., F. Boano, C. Camporeale, L. Ridolfi. 2008. Intra‐meander hyporheic flow in

alluvial rivers. Water Resources Research, 44, W1242,.

https://doi.org/10.1029/2008WR007081.

Revelli, R., F. Boano, C. Camporeale, L. Ridolfi. 2008. Intra‐meander hyporheic flow in

alluvial rivers. Water Resources Research, 44, W12428, doi:10.1029/2008WR007081.

Reynolds O. 1883. An Experimental Investigation of the Circumstances Which Determine

Whether the Motion of Water Shall Be Direct or Sinuous, and of the Law of Resistance in

Parallel Channels.” Philosophical Transactions of the Royal Society of London 174.

Ribes, A., and C. Caremoli. 2007. Salomé platform component model for numerical

simulation, COMPSAC 07: Proceeding of the 31st Annual International Computer Software

and Applications Conference, pages 553-564, Washington, DC, USA, 2007, IEEE Computer

Society.

Richards, L. A. 1931. Capillary conduction of liquids through porous mediums. Physics 1, Nr.

5 318–33, doi:10.1063/1.1745010.

Riml, J., A. Wörman, U. Kunkel, and M. Radke. 2013. Evaluating the fate of six common

pharmaceuticals using a reactive transport model: Insights from a stream tracer test. Science of

The Total Environment, 458-460, 344-354, doi:10.1016/j.scitotenv.2013.03.077.

Bibliography

133

Roche, K. R., G. Blois, J. L. Best, , K. T. Christensen, A. F. Aubeneau, and A. I. Packman.

2018. Turbulence Links Momentum and Solute Exchange in Coarse-Grained Streambeds.

Water Resources Research, 54(5), 3225-3242, doi:10.1029/2017wr021992.

Rodrigues, M. A., L. Pardela, V. Geraldes, J. Santos, H. A. Matos, E. J. Azevedo. 2011.

Theophylline Polymorphs by Atomization of Supercritical Antisolvent Induced Suspensions.

Journal of Supercritical Fluids 58 (2), 303–312.

Rosenberry, D.O., and J.W. LaBaugh. 2008. Field techniques for estimating water fluxes

between surface water and ground water: U.S. Geological Survey Techniques and Methods 4–

D2. p. 128.

Rosenzweig, R., and U. Shavit. 2007. The laminar flow field at the interface of a Sierpinski

carpet configuration, Water Resources Research 43,W10402, doi:10.1029/2006WR005801.

Roskosch, A., M. Hupfer, G. Nützmann, and J. Lewandowski. 2011. Measurement Techniques

for Quantification of Pumping Activity of Invertebrates in Small Burrows. Fundamental and

Applied Limnology / Archiv Für Hydrobiologie 178, no. 2, 89–110.

Roskosch, A., M.R. Morad, A. Khalili, and J. Lewandowski. 2010. Bioirrigation by

Chironomus Plumosus : Advective Flow Investigated by Particle Image Velocimetry. Journal

of the North American Benthological Society 29, no. 3, 789–802.

Roskosch, A., N. Hette, M. Hupfer, and J. Lewandowski. 2012. Alteration of Chironomus

plumosus ventilation activity and bioirrigation-mediated benthic fluxes by changes in

temperature, oxygen concentration, and seasonal variations. Freshwater Science, 31(2), 269–

281.

Saenger, N., P. K. Kitanidis, and R. Street. 2005. A numerical study of surface-subsurface

exchange processes at a riffle-pool pair in the Lahn River, Germany. Water Resour. Res. 41

(12):12424, doi:10.1029/2004WR003875.

Salehin, M., A.I. Packman, and M. Paradis. 2004. Hyporheic exchange with heterogeneous

streambeds: Laboratory experiments and modeling. Water Resources Research 40, no. 11:

W11504.

Sawyer, A.H., M. Bayani Cardenas, and J. Buttles. 2011. Hyporheic exchange due to channel-

spanning logs. Water Resources Research 47, no. 8: W08502.

Schaper, J. L., M. Posselt, C. Bouchez, A. Jaeger, G. Nuetzmann, A. Putschew, G. Singer, and

J. Lewandowski. 2019. Fate of Trace Organic Compounds in the Hyporheic Zone: Influence of

Bibliography

134

Retardation, the Benthic Biolayer, and Organic Carbon. Environmental Science & Technology,

53(8), 4224-4234, doi:10.1021/acs.est.8b06231.

Schaper, J. L., M. Posselt, J. L. McCallum, E. W. Banks, A. Hoehne, K. Meinikmann, M. A.

Shanafield, O. Batelaan, and J. Lewandowski. 2018. Hyporheic Exchange Controls Fate of

Trace Organic Compounds in an Urban Stream. Environmental Science & Technology, 52(21),

12285-12294, doi:10.1021/acs.est.8b03117.

Schmitt, P., and B. Elsaesser. 2015. On the use of OpenFOAM to model oscillating wave surge

converters. Ocean Engineering 108: 98–104, doi: 10.1016/j.oceaneng.2015.07.055.

Schulze, L., and C. Thorenz. 2014. The multiphase capabilities of the CFD toolbox Openfoam

for hydraulic engineering applications. ICHE 2014, Hamburg. Bundesanstalt für Wasserbau.

ISBN 978-3-939230-32-8.

Shrivastava, S., M. J. Stewardson, M. Arora. 2021. Sediment Reworking in Streambeds With

Fine Sediment Deposits and Its Influence on Hyporheic Flow Regime, Water Resources

Research, doi:10.1029/2021WR030360, 57, 12.

Silliman, S. E., and D. F. Booth 1993. Analysis of time-series measurements of sediment

temperature for identification of gaining vs. losing portions of Juday Creek, Indiana. Journal of

Hydrology, 146, 131-148, doi:10.1016/0022-1694(93)90273-C.

Simons, F., T. Busse, J. Hou, I. Özgen, and Reinhard Hinkelmann. 2014. A model for overland

flow and associated processes within the hydroinformatics modelling system. Journal of

Hydroinformatics 16, no. 2 (March 1, 2014): 375–91.

Simons, F., T. Busse, J. Hou, I. Özgen, and Reinhard Hinkelmann. 2014. A model for overland

flow and associated processes within the hydroinformatics modelling system. Journal of

Hydroinformatics 16, no. 375–91, doi:10.2166/hydro.2013.173.

Singha, K., A. Pidlisecky, F.D. Day-Lewis, and M.N. Gooseff. 2008. Electrical characterization

of non-Fickian transport in groundwater and hyporheic systems. Water Resources Research 44,

no. 4: W00D07.

Smagorinsky, J. 1963. General circulation experiments with the primitive equations. Monthly

weather review 99-164.

Sobhi Gollo, V., Broecker, T., Lewandowski, J., Nützmann, G. & Hinkelmann R. 2021. An

integral approach to simulate three-dimensional flow in and around a ventilated U-shaped

chironomid dwelled burrow, Journal of Ecohydraulics, doi:10.1080/24705357.2021.1938258.

Bibliography

135

Sobhi Gollo, V., Broecker, T., Lewandowski, J., Nützmann, G. & Hinkelmann R. 2022. Flow

and transport modeling in heterogeneous sediments using an integral approach. Groundwater,

https://doi.org/10.1111/gwat.13275.

Sobhi Gollo, V., Broecker, T., Marx, C., Lewandowski, J., Nützmann, G. & Hinkelmann R.

2022. A comparative study of integral and coupled approaches for modeling hydraulic

exchange processes across a rippled streambed, International Journal on Geomathematics 13,

no 16: 1-27, doi: 10.1007/s13137-022-00206-5.

Sondergaard, M., E. Jeppesen, T. Lauridsen, C. Skov, E. H. Van Nes, R. Roijackers, E.

Lammens, and R. Portielje. 2007. Lake restoration: successes, failures and long‐term effects.

Journal of Applied Ecology 44,1095‐1105.

Song, J., X. Chen, and C. Cheng. 2010. Observation of bioturbation and hyporheic flux in

streambeds. Environmental scientific engineering. China 4, 340.

Sophocleous, M. 2002. Interactions between groundwater and surface water: the state of the

science. Hydrogeology Journal, 10(1), 52-67, doi:10.1007/s10040-001-0170-8.

Stanford, J. A., and J. V. Ward. 1988. The hyporheic habitat of river ecosystems. Nature,

335(6185), 64–66, https://doi.org/10.1038/335064a0.

Stief, P., and A. Schramm. 2010b. Regulation of nitrous oxide emission associated with benthic

invertebrates. Freshwater Biology 55:1647–1657.

Stief, P., L. Polerecky, M. Poulsen, and A. Schramm. 2010a. Control of Nitrous Oxide Emission

from Chironomus Plumosus Larvae by Nitrate and Temperature, Limnology and

Oceanography., 55(2), 872–884.

Stonedahl, S. H., J. W. Harvey, A. Wörman, M. Salehin, and A. I. Packman. 2010. A multiscale

model for integrating hyporheic exchange from ripples to meanders, Water Resources Research

46, no. 12: W12539.

Stonedahl, S. H., J. W. Harvey, J. Detty, A. Aubeneau, and A. I. Packman. 2012. Physical

controls and predictability of stream hyporheic flow evaluated with a multiscale model. Water

Resources Research, 48, W10513, https://doi.org/10.1029/2011WR011582.

Stonedahl, S.H., A.H. Sawyer, F. Stonedahl, C. Reiter, and C. Gibson. 2018. Effect of

Heterogeneous Sediment Distributions on Hyporheic Flow in Physical and Numerical Models.

Groundwater 56, no. 6: 934–46.

Bibliography

136

Storey, R.G., R.R. Fulthorpe, and D.D. Williams. 1999. Perspectives and predictions on the

microbial ecology of the hyporheic zone. Freshwater Biology 41:119-130, doi:10.1046/j.1365-

2427.1999.00377.x.

Stubbington, R., P.J. Wood, and I. Reid. 2011. Spatial variability in the hyporheic zone

refugium of temporary streams. Aquatic science. 73, 499.

Swain, E. D., and E. J. Wexler. 1993. A coupled surface-water and ground-water flow model

for simulation of stream-aquifer interaction. USGS Open-file report 92–138. Tallahassee, FL:

US Geological Survey, doi:10.3133/ofr92138.

Tonina, D., and J. M. Buffington. 2007. Hyporheic exchange in gravel bed rivers with pool‐

riffle morphology: Laboratory experiments and three‐dimensional modeling. Water Resources

Research, 43, W01421, doi:10.1029/2005WR004328.

Tonina, D., and J. M. Buffington. 2009. A three‐dimensional model for analyzing the effects of

salmon redds on hyporheic exchange and egg pocket habitat. Canadian Journal of Fisheries

and Aquatic Sciences, 66(12), 2157–2173, https://doi.org/10.1139/F09‐146.

Trauth, N., C. Schmidt, M. Vieweg, S.E. Oswald, and J.H. Fleckenstein. 2015. Hydraulic

controls of in-stream gravel bar hyporheic exchange and reactions. Water Resour. Res. 51:

2243-2263, doi:10.1002/2014WR015857.

Trauth, N., C. Schmidt, M. Vieweg, U. Maier and J. H. Fleckenstein. 2014. Hyporheic transport

and biogeochemical reactions in pool-riffle systems under varying ambient groundwater flow

conditions. Journal of Geophysical Research: Biogeosciences, 119(5),

doi:10.1002/2013jg002586.

Trauth, N., C. Schmidt, U. Maier, M. Vieweg, and J. Fleckenstein. 2013. Coupled 3-D stream

flow and hyporheic flow model under varying stream and ambient groundwater flow conditions

in a pool-riffle system. Water Resources Research, 49(9), 5834-5850, doi:10.1002/wrcr.20442.

Valett, H. M., J. A. Morrice, C. N. Dahm, and M. E. Campana. 1996. Parent lithology, surface-

groundwater exchange, and nitrate retention in headwater streams. Limnology and

Oceanography 41, no. 2: 333–345.

Van Gent, M. 1995. Wave Interaction with Permeable Coastal Structures, Elsevier Science:

Amsterdam, The Netherlands, Volume 95.

Bibliography

137

Van Loo, D., L. Bouckaert, O. Leroux, E. Pauwels, M. Dierick, L. Van Hoorebeke, V. Cnudde,

S. De Neve, S. Sleutel. 2014. Contrast agents for soil investigation with X-ray computed

tomography. Geoderma 213:485-491, doi:10.1016/j.geoderma.2013.08.036.

Vaux, W. G. 1968. Intragravel flow and interchange of water in a streambed, Fishery Bulletin

66, no. 3: 479–489.

Weller, H. G., G. Tabor, H. Jasak, and C. Fureby. 1998. A tensorial approach to computational

continuum mechanics using object-oriented techniques, Computers in Physics 12 : 620, doi:

10.1063/1.168744.

Weishaupt, K., V. Joekar-Niasar, R. Helmig. 2019. An efficient coupling of free flow and

porous media flow using the pore-network modeling approach. Journal of Computational

Physics: X 1, 100011. https://doi.org/10.1016/j.jcpx.2019.100011.

Westhoff, M. C., T. A. Bogaard Savenije HHG. 2010. Quantifying the effect of in-stream rock

clasts on the retardation of heat along a stream. Advances in Water resources management 33.

Williams, D.D., C.M. Febria, and J.C.Y. Wong. 2010. Ecotonal and other properties of the

hyporheic zone. Fundumental applications in Limnology. 176, 349.

Winter, T. C., J. W. Harvey, O. L. Franke, and W. M. Alley. 1998. Ground water and surface

water; a single resource (1139). Retrieved from, doi:10.3133/cir1139.

Woessner, W. W. 2000. Stream and Fluvial Plain Ground Water Interactions: Rescaling

Hydrogeologic Thought. Groundwater, 38(3), 423-429, doi:10.1111/j.1745-

6584.2000.tb00228.x.

Wood, P.J., A.J. Boulton, and R. Stubbington. 2010). Is the hyporheic zone a refugium for

aquatic macroinvertebrates during severe low flow conditions? Fundumental applications in

Limnology. 176, 377.

Young, I. M., K. Ritz. 2000. Tillage, habitat space and function of soil microbes. Soil Tillage

Research 53 (3-4):201-213, doi:10.1016/s0167-1987(99)00106-3.

Zarnetske, J. P., R. Haggerty, S. M. Wondzell, M. A. Baker. 2011. Dynamics of nitrate

production and removal as a function of residence time in the hyporheic zone. Journal of

Geophysical Research, 116, G01025, https://doi.org/10.1029/2010JG001356.