NUMERICAL SIMULATIONS OF PERMEABLE-WALL TURBULENCE WITH
APPLICATIONS IN HYPORHEIC EXCHANGE
By
Guangchen Shen
A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering โ Doctor of Philosophy
2022
ABSTRACT
NUMERICAL SIMULATIONS OF PERMEABLE-WALL TURBULENCE WITH
APPLICATIONS IN HYPORHEIC EXCHANGE
By
Guangchen Shen
In aquatic environments such as rivers, the exchange of solutes across the interface between the sedi-
ment and the overlying water plays a significant role in controlling biogeochemical processes, which
are important for an array of topics from nutrient transport and cycling to release of greenhouse
gases such as nitrous oxide. Most previous studies on characterizing this exchange are focused on
flows with sediment bedforms much larger than individual sediment grains. The physics at the
pore or grain scale were typically not resolved. The effects of grain roughness on the sediment
bed surface on the transport across the sediment-water interface (SWI), isolated from those of bed
permeability and bedforms, are not well understood. In this work, direct numerical simulations
(DNS) of the connected system of turbulent open-channel flow and pore-resolved sediment flow
are carried out, with different arrangements of grains at the sediment surface.
First, the statistics and structure of the mean flow and turbulence are characterized in flows
with a friction Reynolds number of 395 and a permeability Reynolds number of 2.6 over sediments
with either regular or random grain packing on a macroscopically flat bed. It is shown that, even
in the absence of any bedform, the subtle details of grain roughness alone can significantly affect
the dynamics of turbulence and the time-mean flow. Such effects translate to large differences in
penetration depths, apparent permeabilities, vertical mass fluxes and subsurface flow paths. The
less organized distribution of mean recirculation regions near the interface with a random packing
leads to a more isotropic form-induced stress tensor, which plays a significant role in increasing
mixing and wall-normal exchange of mass and momentum.
Next, the mass exchange is characterized in detail for macroscopically flat river beds, focusing
on the transit timeโthe time spent by a fluid particle in the sedimentโwhich determines the role
of hyporheic zones in transforming the chemical signature of stream water. Results show that bed
roughness leads to interfacial pressure variations, which induces deep subsurface flow paths that
yield a transit time distribution with a heavy tail. Furthermore, the addition of molecular diffusion
is accounted for and is shown to increase transit times regardless of roughness texture. The results
demonstrate that particle roughness on a macroscopically flat sediment bed can induce significant
hyporheic exchange that is fundamentally similar to that induced by bedforms.
Lastly, to identify possible interaction between the effect of grain roughness and that of a
bedform, DNSs of open channels with a friction Reynolds number of 1580 on a porous dune
with two different roughnesses are conducted. Results show that the roughness modifies the wall
friction, shear penetration depth and pressure distribution along the interface. Unlike the case
on a macroscopically flat bed where the random roughness induces more intense roughness-scale
pressure variation than the regular roughness, over a bedform the random roughness reduces the
macroscopic pressure distribution at the interface instead due to its higher hydrodynamic drag.
The weaker pressure variation in turn weakens the pumping and shortens transit times. The
results highlight the nonlinear interaction between the effects of bed morphological features of
different scales. Pore-resolved simulations such as the ones herein can be used in the future in
direct characterization of pore-scale dynamics to provide insights for pore-unresolved modeling of
biogeochemical processes.
Copyright by
GUANGCHEN SHEN
2022
ACKNOWLEDGEMENTS
I am fortunate to have had support and guidance from many people during the pursuit of my doctorate
degree. First and foremost, I would like to thank my advisor, Professor Junlin Yuan, for her support,
encouragement and guidance throughout my education at the Michigan State University. Thank you
for always believing in me and teaching me how to be an excellent researcher. I am also grateful to
Professor Mantha S. Phanikumar. Your expertise and knowledge in hyporheic exchange has helped
me shape my current understanding in this field. Thanks also go to my dissertation committee.
Professor Farhad Jaberi, whose courses on Turbulent Flows made a significant impact on my study
of turbulence. I would also like to thank my committee member Professor Ricardo Mejia-Alvarez
for your valuable comments and suggestions on my thesis. Working with Professor Mejia-Alvarez
as two semesters teaching assistant has been very memorable and enjoyable.
I would like to thank my colleagues at our TSM lab Saurabh Pargal, Sai Chaitanya and Mostafa
Aghaei Jouybari for their suggestions and help throughout the years. I will never forget the days
we spent in Atalanta, Maryland, Seattle and Phoenix. Special thanks to my friends who help and
share happiness with me throughout these years: Amirreza Gandomkar, Sheng Chen, Yuan Yao,
and Joy Mojumder.
Finally, I would like to thank my parents for their endless love and support and also my sister
who always believes in her brother.
v
TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . xiv
CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Hyporheic exchange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Turbulence across the sediment-water interface . . . . . . . . . . . . . . . 2
1.2.2 Conservative-solute transport in the sediment . . . . . . . . . . . . . . . . 5
1.3 Research objectives and outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
CHAPTER 2 DIRECT NUMERICAL SIMULATIONS OF TURBULENCE AND HY-
PORHEIC MIXING NEAR SEDIMENTโWATER INTERFACES . . . . . . 11
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.3 Synthesis of porous beds . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.4 Geometrical comparison of the two interface types . . . . . . . . . . . . . 17
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.1 Validation of sediment synthesis . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.2 Effects of interface irregularity on flow statistics . . . . . . . . . . . . . . . 23
2.3.2.1 DA velocity and friction . . . . . . . . . . . . . . . . . . . . . . 23
2.3.2.2 Turbulent fluctuations . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.2.3 Form-induced fluctuations . . . . . . . . . . . . . . . . . . . . . 28
2.3.2.4 Reynolds-stress budgets . . . . . . . . . . . . . . . . . . . . . . 30
2.3.3 Implications for the mechanism of flat-bed hyporheic exchange . . . . . . . 34
2.3.3.1 Variation of apparent permeability . . . . . . . . . . . . . . . . . 34
2.3.3.2 Momentum transport mechanisms and dependence on rough-
ness geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.3.3 Mechanisms producing vertical mass flux . . . . . . . . . . . . . 39
2.4 Conclusions and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
CHAPTER 3 QUANTIFYING THE EFFECTS OF BED ROUGHNESS ON TRANSIT
TIME DISTRIBUTIONS VIA DIRECT NUMERICAL SIMULATIONS
OF TURBULENT HYPORHEIC EXCHANGE . . . . . . . . . . . . . . . . 45
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2.1 Case configuration and parameters . . . . . . . . . . . . . . . . . . . . . . 46
3.2.2 Transit time calculation based on particle tracking . . . . . . . . . . . . . . 49
vi
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.3.1 Effect of bed roughness geometry on transit times . . . . . . . . . . . . . . 53
3.3.2 Interfacial ๐ฃ distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3.3 Effect of molecular diffusion on TTD . . . . . . . . . . . . . . . . . . . . 62
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
CHAPTER 4 PORE-RESOLVED DIRECT NUMERICAL SIMULATIONS OF HY-
PORHEIC EXCHANGE INDUCED BY BEDFORMS AND BED ROUGH-
NESS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2.1 Configurations and parameters . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2.2 Transit time calculation based on particle tracking . . . . . . . . . . . . . . 71
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.3.1 Effects of roughness geometry on velocity . . . . . . . . . . . . . . . . . . 73
4.3.2 Effects of roughness geometry on pressure field . . . . . . . . . . . . . . . 74
4.3.3 Effects of roughness geometry on transit times . . . . . . . . . . . . . . . . 77
4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
CHAPTER 5 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
APPENDIX A DETAILS OF POROUS-BED SYNTHESES . . . . . . . . . . . . . 84
APPENDIX B SUPPORT INFORMATION FOR CHAPTER 3 . . . . . . . . . . . 86
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
vii
LIST OF TABLES
Table 2.1: Summary of simulations. ๐ ๐๐ฃ๐ is the average porosity in the bulk of the
sediment; ๐ท is the sphere diameter; ๐ป๐ is the sediment depth; ๐ฟ ๐ฅ๐ is the
domain size in ๐ฅ๐ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Table 2.2: Flow parameters for Cases 2 and 3. ๐ฟ ๐ and ๐ฟ ๐ are Brinkman-layer thickness
and Reynolds-stress penetration depth, both measured from ๐ฆ = โ๐; ๐ฟโ๐ and ๐ฟโ๐
are the same lengths measured from crest; ๐พ โ is apparent permeability. . . . . . 22
Table 3.1: Summary of simulation parameters. ๐ท is the grain diameter; ๐ฟ ๐ฅ๐ is the simu-
lation domain size in ๐ฅ๐ . ฮ๐ฅ + , ฮ๐ฆ + , and ฮ๐ง+ are DNS grid sizes in ๐ฅ, ๐ฆ and ๐ง,
respectively, normalized using the viscous length scale ๐/๐ข ๐ . . . . . . . . . . . . 47
Table 3.2: Various characteristics of transit times calculated based on mean advection
only and both mean advection and molecular diffusion, at two elevations of
particle entry (๐ฆ 1 ) and exit (๐ฆ 2 ). Transit times are normalized by ๐ฟ/๐ข ๐ . Note.
โAโ denotes calculation based on time-mean advection alone; โMDโ denotes
calculation accounting for both time-mean advection and molecular diffusion
(discussed in Section 3.3.3). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
Table 4.1: Summary of parameters. ๐ท is the grain diameter; ๐ฟ ๐ฅ๐ is the simulation domain
size in ๐ฅ๐ ; ๐
๐ ๐พ = 2.5, ๐
๐ ๐ = 1580, ๐ท + = 79. ฮ๐ฅ + , ฮ๐ฆ + , and ฮ๐ง+ are DNS grid
sizes in ๐ฅ, ๐ฆ and ๐ง, respectively, normalized using the viscous length scale ๐/๐ข ๐ . . 69
Table 4.2: Various characteristics of transit times calculated based on mean advection
only. Transit times are normalized by ๐ฟ/๐ข ๐ . Note. ๐50 and ๐95 represent 50th
(median) and 95th percentiles of the TTD distributions. . . . . . . . . . . . . . 77
viii
LIST OF FIGURES
Figure 1.1: Perspectives in (a) reach scale (Boano et al., 2014), (b) local scale (Stonedahl
et al., 2010) and (c) pore scale (Voermans et al., 2018) of surface-subsurface
water interactions and hyporheic exchange. . . . . . . . . . . . . . . . . . . . 1
Figure 2.1: (a) Simulation domain and (b) positions of sediment crest, ๐ ๐ , and zero-plane
displacement, ๐, defined in Section 2.3.2.1. . . . . . . . . . . . . . . . . . . . 14
Figure 2.2: Synthesized sediment beds colored by ๐ฆ/๐ฟ with (a) the regular and (b) the
random interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Figure 2.3: Wall-normal profiles of plane-averaged porosity profile for regular ( ) and
random ( ) interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Figure 2.4: Bed surface height fluctuations for the (a) regular and (b) random interfaces.
(c) Autocorrelation of height fluctuations and Taylor micro scales. . . . . . . . 18
Figure 2.5: Histogram of bed height fluctuations for the (a) regular and (b) random interfaces. 19
Figure 2.6: Comparison of (a) mean velocity and (b-d) streamwise, wall-normal and shear
components of the Reynolds stress tensor. โฆ Experiment, DNS. . . . . . 20
Figure 2.7: Comparison of (a) streamwise, (b) wall-normal and (c) shear component
of the form-induced stress tensor. โฆ Experiment, DNS. Black lines
represent the values spatially averaged over the entire horizontal plane; grey
lines emulate the sampling used in the experiment. . . . . . . . . . . . . . . . 21
Figure 2.8: (a) Diagnostic function for fitting the logarithmic profile and (b) DA velocity
profiles for regular ( ) and random ( ) interfaces; (red) fitted
logarithmic profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Figure 2.9: DA velocity along ๐ฆ in linear scaling for regular ( ) and random ( )
interfaces. Inset magnifies the distribution at ๐ฆ โ โ๐. . . . . . . . . . . . . . . 23
Figure 2.10: Components of (a) the Reynolds-stress tensor and (b) the anisotropy tensor
for regular ( ) and random ( ) interfaces. . . . . . . . . . . . . . . . . 25
Figure 2.11: Contours of instantaneous ๐ขโฒ+ in the (๐ฅ, ๐ง) plane at the respective โจ๐ขโฒ2 โฉ-peak
elevations for the (a) regular and (b) random interfaces. . . . . . . . . . . . . . 26
ix
Figure 2.12: Integral length scales of ๐ขโฒ motions along the streamwise directions for regular
( ) and random ( ) interfaces. . . . . . . . . . . . . . . . . . . . . . . 26
Figure 2.13: Components of the form-induced stress tensor: (a) streamwise, (b) wall-
normal and (c) spanwise components, for regular ( ) and random ( )
interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
Figure 2.14: (a) Time-averaged streamlines for the regular interface. (b) Streamlines on
Slice A (๐ฅโ๐ฆ plane across the crests of the grains). (c) Streamlines distribution
on Slice B (๐ฅ โ ๐ฆ plane across the valley between two rows of grains). Dash
lines indicate ๐ฆ = โ๐. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Figure 2.15: (a) Time-averaged streamlines for the random interface. (b) and (c) show
streamlines on two ๐ฅ โ ๐ฆ planes, Slice C and Slice D. Dash lines indicate ๐ฆ = โ๐. 29
Figure 2.16: Shear stresses for cases with (a) regular and (b) random interfaces:
Reynolds shear stress, viscous shear stress and (red) form-induced
shear stress. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
Figure 2.17: Terms in the TKE budget equation. (a) Production and viscous dissipation
for regular (symbols and thin lines) and random-interface (thick lines) cases;
๐๐ ( , โณ) , ๐๐ + ๐๐ค ( , โ, red) and ๐ ( , โก). Transport terms for
(b) regular and (c) random interfaces; ๐๐ค ( , red), ๐๐ก ( ), ๐๐ ( ) and
๐๐ ( 3
). All terms are normalized by ๐ข ๐ /๐ฟ. . . . . . . . . . . . . . . . . . . . 31
Figure 2.18: (a) Form-induced production and (b) pressure-strain-rate term in regular ( )
and random interface ( ) cases for the budgets of โจ๐ขโฒ2 โฉ๐ (โ, red), โจ๐ฃ โฒ2 โฉ๐
(โก) and โจ๐ค โฒ2 โฉ๐ (โณ, blue); normalization is done with ๐ข ๐ and ๐ฟ. . . . . . . . . . 32
Figure 2.19: Time-averaged spanwise vorticity ๐ ๐ง normalized by ๐ข ๐ /๐ฟ at (a) Slice B of the
regular interface and (b) Slice D of the random interface. . . . . . . . . . . . . 33
Figure 2.20: Horizontal variation of the apparent permeability ๐พ โ in the Brinkman layer
in the random-interface case. . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Figure 2.21: Terms in the DA ๐ข-momentum equation (equation (2.16)), normalized by
๐ข 2๐ /๐ฟ for (a) regular and (b) random interfaces: Reynolds-stress term( ,
red), viscous stress term ( , red), total drag ( ), form-induced stress
term ( ) and pressure-gradient term( ). . . . . . . . . . . . . . . . . . . . . 36
Figure 2.22: Horizontal variation of local volume-averaged D (a) and |I| (b), both nor-
malized by ๐ข 2๐ /๐ฟ. (c) sediment surface height fluctuations โโฒ normalized by
๐ฟ for the random-interface case. (d) Correlations of D โฒ (red) and |I|โฒ (black)
with โโฒ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
x
Figure 2.23: (a) Wall-normal flux at ๐ฆ = โ๐ normalized by ๐ข ๐ ๐ฟ2 with different averaging
window sizes ๐๐ข ๐ /๐ฟ for regular ( ) and random ( ) interfaces. (b)
Pressure-work term of โจe 2 3
๐ฃ โฉ๐ budget normalized by ๐ข ๐ /๐ฟ for regular( , โ)
and random( ) interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Figure 2.24: Time-averaged pressure variation and 3D mean streamlines for (a) regular
and (b) random interfaces. Blue arrows indicate examples of hyporheic flow
paths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Figure 3.1: Simulation domain and synthesized sediment beds colored by ๐ฆ in (a) the reg-
ular and (b) the random cases. (c) One-dimensional power spectral densities
๐ธ ๐ of the interface height fluctuations for Case 1 (regular case, ) and
Case 2 (random case, , red). (red) Slope of power-law decay for
Case 2; ๐ ๐๐๐ is the root-mean-square (rms) of interface height fluctuations. . . 46
Figure 3.2: (a) Sketch of typical tracked flow paths for the random interface (Case 2). The
cumulative backward transit time distribution (b) and its pdf (c) are shown for
two sets of entry and exit locations: ๐ฆ 1 = ๐ฆ 2 = 0 ( ) and ๐ฆ 1 = ๐ฆ 2 = โ3๐ท
( ) for the regular interface (black lines) and the random interface (red
lines). (d) Modified pdf showing the probability density in evenly spaced
logarithmic increments. In (b), the horizontal dashed lines ( , blue) mark
the 50th (median) and 95th percentiles. In (c), the dashed blue lines ( ,
blue) show the fitted power law with slopes indicated. . . . . . . . . . . . . . . 52
Figure 3.3: Time-mean subsurface flow paths corresponding to the median (50th-percentile)
transit time (a,b,e,f) and the 95th-percentile transit time (c,d,g,h) for entry/exit
elevations at ๐ฆ 1 = ๐ฆ 2 = 0 (a,b,c,d) and ๐ฆ 1 = ๐ฆ 2 = โ3๐ท (e,f,g,h), for the regular
(a,c,e,g) and random (b,d,f,h) interfaces. . . . . . . . . . . . . . . . . . . . . . 55
Figure 3.4: Wall-normal variation of ๐ฃ magnitude conditionally averaged along down-
ward (red) or upward (black) flow paths tracked at ๐ฆ 1 = ๐ฆ 2 = 0, for regular
( ) and random ( ) cases. . . . . . . . . . . . . . . . . . . . . . . . . . . 57
Figure 3.5: Median transit time normalized by ๐ฟ/๐ข ๐ with independent variations of ๐ฆ 1
(entry) and ๐ฆ 2 (exit) elevations, for regular (a) and random (b) cases. Similar
variations of ๐95 are shown in the Appendix B Figure B.5. . . . . . . . . . . . 58
Figure 3.6: Sketch showing how local high-pressure (stagnation point) and low-pressure
(recirculation) regions lead to entry and exit of fluid parcels (red spheres and
red lines), specifically, in Case 2. Roughness geometry. ๐(๐ฅ, ๐ฆ, ๐ง) and
๐ฃ(๐ฅ, ๐ฆ, ๐ง) are the time-averaged values of pressure and wall-normal velocity;
e ๐ฆ, ๐ง) = ๐(๐ฅ, ๐ฆ, ๐ง) โ โจ๐โฉ(๐ฆ) is the form-induced pressure, i.e., deviations
๐(๐ฅ,
from the plane-averaged value. . . . . . . . . . . . . . . . . . . . . . . . . . . 59
xi
Figure 3.7: Bed-parallel contours for the random interface (Case 2) at ๐ฆ = 0: (a) bed-
normal mean velocity, (b) form-induced pressure, (c) bed-normal gradient of
mean pressure, and (d) transit time of flow paths entering at each location.
Only 1/4 of the domain is shown. All quantities in (a-d) are normalized by ๐ฟ
and ๐ข ๐ . Joint probability density functions of ๐๐/๐๐ฆ and ๐ฃ in the interface
region, for Case 1 (e) and Case 2 (f). . . . . . . . . . . . . . . . . . . . . . . . 60
Figure 3.8: Cumulative backward transit time distribution (a) and its pdf (b) based on
mean advection only ( ) and based on both mean advection and molecular
diffusion ( ), for Case 1 (black) and Case 2 (red) at ๐ฆ 1 = ๐ฆ 2 = 0. In (b),
(blue) shows the fitted power law with slope indicated. . . . . . . . . . . . . 63
Figure 4.1: Simulation domain and synthesized bedforms with (a) regular and (b) random
roughnesses, colored by ๐ฆ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Figure 4.2: Roughness surfaces colored by ๐ฆ (a,b) and local height fluctuations measured
from the bedform surface (c,d) in Cases 1 (a,c) and Case 2 (b,d). (e) One-
dimensional power spectral densities ๐ธ ๐ of the local height fluctuations for
Case 1 (regular roughness, ) and Case 2 (random roughness, , red).
๐ ๐๐๐ is the root-mean-square (rms) of interface height fluctuations. . . . . . . . 70
Figure 4.3: Conceptual sketch to show the enter and exit locations of tracked particles. . . 71
Figure 4.4: (a) Time- and spanwise-averaged streamwise velocity โจ๐ขโฉ +๐ง for cases with
regular ( ) and random ( , red) roughnesses. โจ๐ขโฉ +๐ง contours for (b)
regular and (c) random roughness. . . . . . . . . . . . . . . . . . . . . . . . . 72
Figure 4.5: Profiles of (a) viscous shear stresses โ๐โจ๐ขโฉ ยฏ ๐ง+ /๐๐ฆ + and (b) Reynolds shear
stresses โโจ๐ขโฒ๐ฃ โฒโฉ๐ง+ for cases with regular ( ) and random ( , red)
roughnesses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Figure 4.6: Time- and spanwise-averaged pressure field for cases with (a) regular and (b)
random roughnesses. Pressure distributions along the SWI normalized using
bulk velocity (c) or friction velocity (d) for both cases. . . . . . . . . . . . . . . 75
Figure 4.7: Distribution of the Clauser parameter ๐ฝ in the two cases. The black dash line
shows the bedform interface geometry. . . . . . . . . . . . . . . . . . . . . . . 76
Figure 4.8: Cumulative backward transit time distributions (a) and their pdf distributions
(b) for regular ( ) and random ( , red) cases. In (a), the horizontal
dashed lines ( , blue) mark the 50th (median) and 95th percentiles. In
(b), the dashed blue lines ( , blue) show the fitted power law with slopes
indicated. Subsurface flow paths in the (c) regular and (d) random cases. . . . . 78
xii
Figure 4.9: Histogram of the number of flow paths reaching a depth of a given ๐ฆ value
for the two cases. Values are normalized such that the area integral is one. . . . 79
Figure A.1: Bed synthesis with random interface packing. . . . . . . . . . . . . . . . . . . 84
Figure A.2: Bed synthesis with regular interface packing. . . . . . . . . . . . . . . . . . . . 84
Figure B.1: Residence time distributions (RTD) normalized by ๐/๐ for ๐ฆ 1 = ๐ฆ 2 = 0
( ) and โ3๐ท ( ). For ๐๐ข ๐ /๐ฟ > ๐(106 ) when ๐ฆ 1 = ๐ฆ 2 = 0 or ๐(107 )
when ๐ฆ 1 = ๐ฆ 2 = โ3๐ท, the RTD rapidly decreases; this is not physical but an
artifact associated with the finite domain size in ๐ฅ (limited to 2๐ฟ ๐ฅ ) imposed
for particle tracking as the tracking procedure cannot be indefinitely long. . . . 86
Figure B.2: Validation of particle tracking against EB97, with various spatial resolution
and a sediment domain of ๐ฆ โ [โ๐/2, 0] (a) or ๐ฆ โ [โ๐, 0] (b). . . . . . . . . . 88
Figure B.3: Sketch to demonstrate that a particle entry elevation (๐ฆ 1 ) lower than the exit
elevation (๐ฆ 2 ) tracks a larger fraction of subsurface flow paths than in the case
with ๐ฆ 1 > ๐ฆ 2 , considering the same |๐ฆ 1 โ ๐ฆ 2 |. Here, (๐ฆโฒ1 , ๐ฆโฒโฒ2 ) tracks 2 out of
the 4 paths (50%) entering from ๐ฆโฒ1 , while (๐ฆโฒโฒ1 , ๐ฆโฒ2 ) tracks 2 out of the 2 paths
(100%) entering from ๐ฆโฒโฒ1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Figure B.4: Streamwise dispersive velocity intensities in the sediment bed: comparison
between results obtained from DNS data taken for simulation times of 10
( , red) and 20 ( ) large-eddy turn-over times, for the regular- (a) and
random-interface (b) cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Figure B.5: Transit time of the 95th percentile (๐95 ) normalized by ๐ฟ/๐ข ๐ with independent
variations of ๐ฆ 1 (entry) and ๐ฆ 2 (exit) elevations, for the regular- (a) and
random-interface (b) cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
xiii
KEY TO SYMBOLS AND ABBREVIATIONS
๐ด ๐ area occupied by fluid
๐ด๐ total area of fluid and solid
๐ฝ Clauser paramete
๐๐ ๐ Reynolds stress anisotropy tensor
๐ต mean-velocity profile intercept in the logarithmic region
๐ถ ๐ friction coefficient
๐ zero-plane displacement
๐ท grain diameter
๐ท ๐ molecular diffusion coefficient
๐ท + particle Reynolds number
ฮ๐ก time differential
ฮ๐ฅ, ฮ๐ฆ, ฮ๐ง grid spacing
ฮ๐ + roughness function
๐ฟ boundary layer thickness
๐ฟ ๐ Brinkman-layer thickness
๐ฟ ๐ Reynolds-stress penetration depth
๐ฟโ displacement thickness
๐ฟ ๐ viscous length scale
๐ฟ๐ ๐ Kronecker Delta
๐ viscous dissipation
๐ธ ๐ power spectral density
๐น๐ immersed-boundary-method body force in ๐ direction
๐น๐ท total drag
โ sediment height
๐ป๐ sediment depth
xiv
๐ Taylor micro-scale of local roughness-height fluctuations
๐ + roughness Reynolds number
๐ ๐ sediment crest elevation
๐
von Kรกrmรกn constant
๐พ permeability
๐พ โ aparent permeability
๐ wavelength
๐ฟ ๐ฅ , ๐ฟ ๐ฆ , ๐ฟ ๐ง domain size
๐ฟ ๐ ๐ integral length scale
๐ kinematic viscosity
๐๐ , ๐ ๐ , ๐ ๐ numbers of grid points in ๐ฅ, ๐ฆ, and ๐ง
๐
autocorrelation
๐ ๐ transit time probability density function
๐๐ transit time cumulative distribution function
๐ ๐ residence time probability density function
๐ or ๐ modified pressure (pressure over density)
๐๐ shear production
๐๐ค wake production
ฮ pressure work
๐ instantaneous flow variable
๐ vertical volumetric flux
๐ density
R pressure-strain-rate term
๐
๐ Reynolds number
๐
๐ ๐ friction Reynolds number
๐
๐ ๐พ permeability Reynolds number
๐ instantaneous total water storage
xv
๐โ root-mean-square of height fluctuation
๐ก time
๐ total simulation time
๐๐ pressure transport
๐๐ก turbulent transport
๐๐ฃ viscous transport
๐ transit time
๐ porosity
๐ ๐๐ฃ๐ average porosity in the bulk of the sediment
๐ข๐ instantaneous velocity components
๐ข ๐ friction velocity
๐ mean velocity
๐๐ฟ free-surface (or channel-half-height) flow velocity
๐ ๐ constant velocity in the sediment far below the SWI
๐ turbulent vorticity
๐ฅ๐ spatial direction
๐ด๐๐บ adverse pressure gradient
๐ถ ๐ด๐ท computer-aided design
๐ท ๐ด double-averaing
๐ท๐ ๐ direct numerical simulation
๐น๐๐บ favourable pressure gradient
๐ผ ๐ต๐ immersed boundary method
๐พ โ ๐ป Kelvin-Helmholtz
๐ฟ๐ธ ๐ large-eddy simulation
๐ฟ๐ธ๐๐๐ large-eddy turn-over time
๐ ๐ท molecular dynamics
๐ ๐๐ผ message passing interface
xvi
๐๐๐ probability distribution function
๐
๐ด๐ ๐ Reynolds-averaged Navier-Stokes
๐
๐๐ root-mean-square
๐
๐ ๐ท residence time distribution
๐๐ ๐ผ sediment-water interface
๐ ๐พ ๐ธ turbulent kinetic energy
๐๐ ๐ท transit time distribution
xvii
CHAPTER 1
INTRODUCTION
In this Chapter, first the hyporheic exchange and its importance are introduced. Then a detailed
literature review is conducted and limitations of the current understanding of the exchange processes
are discussed. The research objectives are proposed based on present limitations and the studies in
this dissertation to address these goals are outlined.
1.1 Hyporheic exchange
Aquatic environments such as rivers form an important part of the water circulation on earth. A
typical river flow can be considered as a turbulent water flow bounded by permeable sediment
beds at the sides and the bottom, as shown in Fig. 1.1 (a). Vertical and lateral exchanges of water
Figure 1.1: Perspectives in (a) reach scale (Boano et al., 2014), (b) local scale (Stonedahl et al.,
2010) and (c) pore scale (Voermans et al., 2018) of surface-subsurface water interactions and
hyporheic exchange.
1
and solutes with the more slowly moving flow in the sediments are important. Specifically, the
sediments slow down the overall movement of solutes and increase opportunities for biogeochemical
processing (Boano et al., 2014). The hyporheic zone is the portion of sediments surrounding the
stream that is permeated with stream water. The exchange of water, solutes and momentum between
the subsurface and the surface flows is called the hyporheic exchange (Bencala, 2000; Wondzell,
2006). Hyporheic exchange plays a significant role in controlling biogeochemical processes in
aquatic environments. As shown in Fig. 1.1, the interfacial exchange phenomena occur over a
wide range of spatial and temporal scales. The hyporheic exchange is affected by a number of
key processes and parameters including sediment roughness, grain size (Malard et al., 2002),
bedform, bed permeability, bed heterogeneity, water flow velocity and unsteadiness (Boano et al.,
2014; Cardenas, 2015) and directly influences the fluxes of water and solutes (nutrients, dissolved
organic carbon and oxygen, etc.) by transferring momentum and energy across the sediment-water
interface.
1.2 Literature Review
1.2.1 Turbulence across the sediment-water interface
At a fundamental level, a river flow can be modeled as a turbulent open channel flow of water
overlaying a porous wall. A review of existing studies on turbulence bounded by permeable walls
is provided next.
Studies of turbulent flow over permeable beds have addressed the effect of wall permeability
on canonical wall turbulence (such as a fully developed turbulent channel flow). These include, for
examples, experimental studies of Zagni and Smith (1976), Kong and Schetz (1982), Zippe and Graf
(1983), Manes et al. (2011), Kim et al. (2018) and numerical studies of Kuwata and Suga (2016a),
Kuwata and Suga (2016b), Fang et al. (2018). Wall permeability is shown to significantly affect the
dynamics and statistics of turbulence. Specifically, the effects include a higher friction coefficient,
a permeability-dependent von Kรกrmรกn constant (๐
), reduced streamwise velocity fluctuations, and
augmented wall-normal and spanwise velocity fluctuations due to a relaxation of wall blocking.
2
Another focus of study has been on the effects of heterogeneous permeability inside the sediment
(for example, see Laube et al. (2018)). A few studies emphasized the importance of the uppermost
sediment layer on the transport across the interface (for example, see Kalbus et al. (2009)).
The interfacial flows can be grouped into three regimes (Suga et al., 2011; Voermans et al., 2017)
according to the type of dominant flow processes. Such categorization is based on the permeability
โ
Reynolds number, ๐
๐ ๐พ = ๐ข ๐ ๐พ/๐ (where ๐ข ๐ is the friction velocity, ๐พ is the permeability, and
๐ is the kinematic viscosity): (1) an effectively impermeable regime (๐
๐ ๐พ โช 1) where attached
eddies dominate near-wall dynamics, displaying a broad range of turbulent scales (Cossu and
Hwang, 2017); (2) a highly permeable regime (๐
๐ ๐พ โซ 1) where the interface is characterized by a
distinct inflection point of the mean velocity and a predictable frequency of coherent motions due
to Kelvin-Helmholtz instability; and (3) a transitional regime (๐
๐ ๐พ โผ ๐(1)), where characteristics
of both limiting regimes exist. The current work is focused on a scenario similar to river flows over
sand beds, i.e., turbulent water flows over sediments with ๐
๐ ๐พ โผ ๐(1).
A challenge in understanding the momentum and mass transport across the sediment-water
interface is associated with the lack of detailed information on the fluid dynamics at the scale of
sediment grains (or particles), especially near the bed surface. This is due to both the difficulty of
experimental measurements at the scale of sediment grains as well as the high cost of numerical
simulations needed to resolve the full spectrum of scales from individual grains to bedformsโ
ripples, dunes, pool-riffle structures, etc. The wall-normal protuberances formed by individual
grains are characterized as bed roughness (or particle roughness). This is to be differentiated
from the protuberances of bed forms, which are formed by clusters of grains and are typically
orders-of-magnitude larger compared to the length scale associated with grains.
The bed roughness influences properties of time-average flow and turbulence over permeable
walls (Nikora et al., 1998). Many rough, permeable-bed flow studies focused on the role of
permeability using similarly structured sediments that vary in particle size; in this case, it is
difficult to separate the effect of bed roughness from that of permeability, except for cases with very
small particle sizes (Breugem et al., 2006) or cases where the grain packing structure is modified
3
inside the bed (Fang et al., 2018) only. Manes et al. (2009) compared permeable and impermeable
walls with the same roughness and showed that the penetration of turbulent flows into the porous
bed leads to higher flow resistance than an impermeable rough bed. Even in the fully rough regime,
the friction factor for the permeable bed increases with Reynolds number; this is different from an
impermeable-bed turbulent flow which, in the region near the wall, is Reynolds number independent
in the fully rough regime (Jimรฉnez, 2004).
The understanding on the effect of roughness isolated from that of permeability is limited.
Several existing studies are summarized below. Padhi et al. (2018a,b) experimentally compared
two gravel beds with different bed surface roughnesses. They observed that the bed with streamwise
oriented gravels which had a higher bed surface roughness generates significantly higher Reynolds
stresses, form-induced stresses and turbulent kinetic energy (TKE) flux than the sediment with
randomly oriented gravels . Han et al. (2018) conducted large-eddy simulations (LES) of solute
transport on synthesized sediments consisting of regularly packed spheres with a layer of small
spheres above to simulate roughness in different regimes of rough-walled turbulence. They found
that a high roughness Reynolds number (๐๐ข ๐ /๐, where ๐ is the roughness height) leads to disruption
of the diffusive layer and induces turbulent motions that more efficiently transport the solute than
does a lower roughness Reynolds number. Kim et al. (2020) experimentally investigated surface and
subsurface flows with idealized permeable beds made of cubically packed spheres. Two cases are
considered with one exhibiting a smooth interface while the other that embodied a hemispherical
surface topography at the SWI. The results show that the presence of bed roughness intensifies
the strength and penetration of flow into the permeable bed modulated by large-scale structures
in the surface flow, and linked to possible roughness-formed channelling effects and shedding of
smaller-scale flow structures from the roughness elements.
While most existing studies focusing on the effects of sediment bed roughness were conducted
using regularly packed particles (Manes et al., 2009; Kim et al., 2017; Fang et al., 2018; Kim et al.,
2020), natural permeable bed roughnesses are characterized by random shape, orientation, spacing
and arrangement of particles. These irregularities potentially affect the pattern and intensity of
4
turbulence, as well as the transport characteristics across the interface.
1.2.2 Conservative-solute transport in the sediment
In addition to the roughness effects on the water transport mentioned above, the transport of
solutesโsuch as dissolved oxygen and nitrateโacross the sediment-water interface may be also
important. The role of hyporheic zones in transforming the chemical signature of stream water de-
pends, among other factors, on transport time scales often described using transit time distributions
(TTDs) and residence time distributions (RTDs). For a solute mass pulse entering the river bed
at time ๐ก, the residence time ๐ is the age of a water parcel in the bed (Elliott and Brooks, 1997b).
Transit time, on the other hand, is the residence time (or water age) at the time a parcel of water
exits the system, and therefore represents the maximum residence time (or age) that can be attained
by a water parcel for a given subsurface flow path. RTDs and TTDs are related and they have both
been extensively used in hyporheic zone research. For example, residence times have been used to
evaluate the ability of the hyporheic zone to process nutrients (Boano et al., 2014) and to classify
hyporheic zones (Harvey et al., 2013). Transit times have been used to understand the kinematics
of age mixing processes in advection-dispersion models (Benettin et al., 2013), as well as relations
between flow, storage and chloride transport in a small watershed (Benettin et al., 2015; Harman,
2015), to infer changes in water cycle dynamics in managed landscapes (Danesh-Yazdi et al., 2016)
and to understand biogeochemical reactions in a dam-regulated river corridor (Song et al., 2020).
Median transit times have been used to assess the reaction efficiency of hyporheic flows by defining
a Damkรถhler number representing the ratio of transport and biogeochemical time scales (Song
et al., 2020).
Bed geometry and morphology exert a strong influence on exchange fluxes and transit times.
Bed geometries are often multiscale, with both large- and small-scale features contributing to the
exchange (Persson et al., 2005; Wรถrman et al., 2006; Aubeneau et al., 2015; Stubbs et al., 2018;
Roche et al., 2018; Lee et al., 2020). For example, large-scale features such as bars that scale
with bankfull depth and channel width, relatively small-scale features such as bedforms, as well
5
as even smaller-scale roughness details at the scale of individual grains all influence exchange
fluxes. Existing studies on flow and residence/transit times are typically focused on the effects of
morphological features such as ripples and dunes that are orders of magnitude larger than individual
sediment grains. Elliott and Brooks (1997a,b) conducted analytical studies and lab experiments on
the exchange with permeable bedforms. They summarized that bedforms โturnoverโ. Pumping is
the movement of fluid through the bed induced by static pressure variations over bedforms. Turnover
occurs as a moving bedform trap and releases interstitial fluid in the bed. The main mechanism
driving mass exchange when bedforms are stationary are based on the process of pumping. Elliott
and Brooks (1997b) (henceforth referred to as EB97) analytically solved for the velocity field and
residence time for an ideal subsurface flow model with a sinusoidal bedform. The solutions were
based on a two-dimensional (2D) Darcy equation with prescribed pressure distribution at the bed
surface; these solutions described the pumping mechanism of the exchange. The time a fluid
particle spends along its subsurface flow path was calculated using a numerical particle tracking
approach. For stationary bedforms, it was observed that the sinusoidal head pumping model (Elliott
and Brooks, 1997a) can only be used in the initial stage and underpredicts the mass exchange in
the later stage. They inferred that additional exchange may be related to bed inhomogeneity or
irregular variations in pressure at the bed surface which are not related to bedforms. Such effects are
likely to be more pronounced in actual river flows. The model of EB97 considered only advection
(by the spatially and temporally averaged Darcy velocity) for flow path tracking, while in Elliott
(1991) additional mixing mechanisms (turbulent mixing or molecular diffusion of solutes) were
also accounted for by using a random walk method. In all the above methods, a constant solute
concentration is assumed in the fluid column; in other words, the flow is assumed well mixed. It
is probably different in an actual river flow, as the bed morphology is likely to induced spatial
heterogeneity of local solute concentration.
Examples of the application of the Darcy model for subsurface flow calculation and particle
tracking to estimate residence time based on advective exchange (EB97) include the following,
among many others. Packman et al. (2000) applied a model of pumping exchange to colloid
6
transport and filtration. They calculated the residence time by superimposing advective transport
and particle settling in the bed and included the effect of physicochemical filtration by bed sediment
also. Wรถrman et al. (2002) developed a framework to integrate the residence time estimation
based on the method of Packman et al. (2000) into a model of longitudinal solute transport in
streams. Many studies, such as Salehin et al. (2004),Sawyer and Cardenas (2009) and Laube
et al. (2018), characterized hyporheic exchange within heterogeneous streambeds. It was shown
that, in a heterogeneous bed compared to a homogeneous one, more shorter and faster preferential
hyporheic flow paths are induced, leading to more rapid net hyporheic exchange and a shallower
hyporheic zone. Sawyer and Cardenas (2009) numerically studied the heterogeneity of cross-bedded
sediment in which sedimentary structures are roughly horizontal units composed of inclined layers.
The results show that the permeability structures do impact the distribution of residence times.
Specifically, permeability heterogeneity can increase the proportion of short or long residence
times. The shape of the residence time distribution appears to be sensitive to permeability patterns
near the sediment-water interface. However, the tails of the residence time distributions follow
a power law, regardless of specific permeability structure. Laube et al. (2018) further developed
a modeling strategy for an equivalent conductivity tensor for heterogenous sediments. Recently,
several studies analyzed sediment flows using pore-scale simulation data (Sun et al., 2015; Kim
and Kang, 2020). For example, Kim and Kang (2020) studied 2D pore-scale transport phenomena
across a surface-subsurface interface, based on a Reynolds-averaged Navier-Stokes (RANS) closure.
Most of the studies summarized above quantified residence time using ideal conceptual models
of bedform features (such as ripples or dunes). Detailed aspects of solute transport near the bed,
including turbulent mixing and pore-scale dynamics, are rarely studied. At the โmicroโ scale of
pore flow, the effects of bed roughness formed by the uppermost-layer of sediment grains and their
arrangement on biogeochemical processes are likely to be important, as the benthic biolayer where
important nutrient or pollutant processing occurs is limited to a thin layer (of 2-5 cm) beneath the
streambed surface (Battin et al., 2008; Harvey et al., 2013; Knapp et al., 2017). In addition, in a
warmer climate streambed roughness is expected to increase as a result of enhanced bioturbationโ
7
the effect of small animals, such as chironomid larvae and tubificid worms, within the sediment
bed as they rework the sediment (Baranov et al., 2016; Dairain et al., 2020).
The fundamental effects of wall roughness on turbulence were predominantly studied for flows
over impermeable walls (see reviews by Jimรฉnez (2004) and Flack and Schultz (2014)). Typically,
roughness is categorized as bed height variations in the range between 5๐ฟ ๐ and 0.1๐ฟ (where
๐ฟ ๐ = ๐/๐ข ๐ is the viscous length scale and ๐ฟ is the turbulent boundary layer thickness). Roughness
affects turbulence by increasing wall friction and enhancing turbulence production near the wall.
It is known to enhance the momentum transfer toward the wall (Jimรฉnez, 2004) and dynamically
modify the near-wall flow (Raupach and Shaw, 1982). The effects become more important in
non-equilibrium flows (Ghodke and Apte, 2016; Yuan and Piomelli, 2014a, 2015; Mangavelli
et al., 2021) where streamwise mean pressure gradientsโwhich can be induced by the presence of
bedformsโare significant. In the context of turbulent flow bounded by a permeable wall (such as
the region near a sediment-water interface), bed roughness likely affects the flow near the bed and
consequently the mass exchange also.
Some insights into how small-scale features of bed topography (e.g., grain-scale irregularities)
may affect the exchange were provided by a few studies on multiscale or fractal bed topographies.
Aubeneau et al. (2015) performed laboratory measurements of solute transport with sand beds
under different flow velocities. They found that fractal properties of the bed topography affect
solute residence time distributions. Larger bedforms induced more hyporheic exchange but shorter
residence times, as shown by the steeper slope of the washout breakthrough curve. He et al.
(2019) applied LES to study the effects of roughness Reynolds number on scalar transport at
SWI. They showed that from hydraulically smooth regime (0 < ๐
๐ ๐ โค 5, where ๐
๐ ๐ is the
roughness Reynolds number) to the transitionally rough regime (5 < ๐
๐ โค 70), the dominant
transport mechanism at the interface changes from molecular diffusion to turbulent diffusion. In
addition, Lee et al. (2020) analyzed the effect of fractal properties of the bed topography on
interfacial flux and residence time, based on sequentially coupled simulations using RANS for the
channel flow and a Darcy model for the subsurface flow. Residence times were calculated using the
8
random walk particle tracking method. They observed that the interfacial flux increases with the
fractal dimension of the bed (i.e., prominence of smaller topographical features). In such coupled
surface-subsurface model framework (used also by Cardenas and Wilson (2007), Cardenas et al.
(2008), and Chen et al. (2015), for examples), the two sets of governing equations are coupled
through the pressure distribution along the sediment-water interface. In this case, the accuracy
of RANS pressure solution along the interface dictates the accuracy of the coupled simulation.
Given the difficulties of RANS simulations in predicting accurately the near-wall flow and flow
separation in strongly non-equilibrium wall-bounded turbulence (Yuan et al., 2014; Dutta et al.,
2016), pore-resolved simulation or measurement data are needed to characterize near-wall dynamics
and provide validation for pore-unresolved modeling.
1.3 Research objectives and outline
The goal of this dissertation is to analyze the effect of morphological details of sediment beds on
interfacial turbulence and solute transport in turbulent flows that are representative of river flows
on sand beds. To this end, three-dimensional DNS simulations are carried out with synthesized
sediments. Both the flow in the water column and the that inside the sediment pores are resolved.
Specific questions are:
1. How different are the effects of roughness on a flow over porous wall versus its effects on
that over a non-permeable wall? How does sediment bed roughness on a macroscopically
flat bed affect the turbulence statistics and structure?
2. How does bed roughness on a macroscopically flat bed influence fluid parcel transit times in
the bed? Does small-scale bed roughness lead to significant fluid volumetric flux normal to
the interface?
3. To what extent will the grain-scale bed roughness affect the hyporheic exchange with the
presence of bedform?
The outline of this dissertation is as follows.
9
โข Chapter 2 describes the methodologies of the sediment synthesis with two different rough-
nesses, regular and random, and the numerical simulations. Simulation parameters are
introduced. The simulations are then validated against existing experimental measurements.
Results on how roughness on a flat bed, isolated from the effect of bed permeability, affects
the interfacial turbulent flow are discussed.
โข Chapter 3 introduces the methodology to calculate the transit time distribution and compares
the transit time distribution for the two different roughnesses. The effect of molecular
diffusion in the exchange is also quantified.
โข Chapter 4 focuses on turbulent flows over a bed with both bedform and roughness. The
bedform synthesis and simulation parameters are introduced. The effects of bed roughness
geometry on velocities, pressure and transit time are compared. Comparisons are made with
flat-bed results discussed in Chapters 3 and 4, to show a different role played by the roughness
when a bedform is present.
โข Chapter 5 summarizes the major findings and contributions, and describes the plan for future
work.
10
CHAPTER 2
DIRECT NUMERICAL SIMULATIONS OF TURBULENCE AND HYPORHEIC
MIXING NEAR SEDIMENTโWATER INTERFACES
2.1 Introduction
The goal of this Chapter is to characterize the detailed effects of different bed roughnesses on
aspects of turbulence that are responsible for scalar and momentum transfer, including averaged
statistics and flow structural information. Specific questions include the following:
(1) How different are the roughness characteristics corresponding to regular and random ar-
rangements of sediment grains in the top layer?
(2) How does the bed roughness texture influence the dynamics at the individual grain level
and the overall interfacial exchange of mass and momentum? Specifically, how important are the
additional, form-induced terms in the momentum and stress balances?
To address the above questions, we perform direct numerical simulations (DNS) of fully-
developed turbulent open-channel flows over grain-resolved permeable beds with an identical bulk
permeability (below the interface region) but different roughnesses produced by either random or
regular grain arrangements at the bed surface. The Reynolds numbers based on the open-channel
height (๐
๐ ๐ = ๐ข ๐ ๐ฟ/๐, where ๐ฟ is the open-channel height) and on the permeability (๐
๐ ๐พ ) are kept
constant. This chapter is organized as follows. Section 2.2 introduces the methodology of the simu-
lations and porous-bed synthesis and compares the two different interfaces generated. Section 2.3.1
validates the random sediment synthesis by comparing with the experimental measurements of Vo-
ermans et al. (2017). Section 2.3.2 compares the turbulence statistics and structure as affected
by the differences in the interface. Lastly, Section 2.3.3 discusses implications for the hyporheic
exchange mechanisms.
11
2.2 Methodology
2.2.1 Governing equations
The incompressible flow of a Newtonian fluid is governed by the equations of conservation of mass
and momentum:
๐๐ข๐
= 0, (2.1)
๐๐ฅ๐
๐๐ข ๐ ๐๐ข๐ ๐ข ๐ ๐๐
+ = โ + ๐โ2 ๐ข ๐ + ๐น ๐ . (2.2)
๐๐ก ๐๐ฅ๐ ๐๐ฅ ๐
Here, ๐ฅ 1 , ๐ฅ2 and ๐ฅ 3 (or ๐ฅ, ๐ฆ and ๐ง) are, respectively, the streamwise, wall-normal and spanwise
directions, and ๐ข ๐ (or ๐ข, ๐ฃ and ๐ค) are the velocity components in those directions; ๐ = ๐/๐
is the modified pressure, where ๐ the density. The term ๐น ๐ in Equation (2.2) is a body force
imposed by an immersed boundary method (IBM) to impose no-slip boundary conditions on the
fluid-solid interface. The grain geometry is well-resolved by the grid. The IBM is based on the
volume-of-fluid approach (Scotti, 2006); its detailed implementation and validation are provided
in Yuan and Piomelli (2014b,c). The ๐น๐ values are negligible except in the cells that are cut by the
immersed solid boundaries. The simulations are performed using a well-validated code that solves
the governing equations (2.1) and (2.2) on a staggered grid using second-order, central differences
for all terms, semi-implicit time advancement (with second-order Adams-Bashforth scheme for the
wall-normal diffusion term) and MPI parallelization (Keating et al., 2004).
Fully developed open-channel flows are simulated with symmetric boundary conditions applied
at both top and bottom boundaries of the simulation domain and periodic conditions are applied in
๐ฅ and ๐ง. A constant pressure gradient is used to drive the flow.
In the sediment, the presence of grains leads to spatial heterogeneity of the time-averaged
variables; these time-averaged fluctuations are separated from turbulent fluctuations using the
double-averaging (DA) decomposition introduced by Raupach and Shaw (1982),
๐(๐ฅ, ๐ก) = โจ๐โฉ(๐ฆ) + ๐e(๐ฅ) + ๐โฒ (๐ฅ, ๐ก), (2.3)
12
Interface type ๐ ๐๐ฃ๐ ๐
๐ ๐พ ๐
๐ ๐ ๐ท/๐ฟ ๐ท+ ๐ป๐ /๐ฟ (๐ฟ ๐ฅ , ๐ฟ ๐ง )/๐ฟ (ฮ๐ฅ + , ฮ๐ฆ +min , ฮ๐ง+ )
Case 1 Random 0.41 2.56 178 0.43 76 2.6 (12, 6) (2.0, 1.5, 2.0)
Case 2 Regular 0.41 2.62 395 0.20 79 0.8 (6, 3) (1.5, 0.19, 1.5)
Case 3 Random 0.41 2.62 395 0.20 79 0.8 (6, 3) (1.5, 0.19, 1.5)
Table 2.1: Summary of simulations. ๐ ๐๐ฃ๐ is the average porosity in the bulk of the sediment; ๐ท is
the sphere diameter; ๐ป๐ is the sediment depth; ๐ฟ ๐ฅ๐ is the domain size in ๐ฅ๐ .
where ๐ is an instantaneous flow variable, โจ๐โฉ is the intrinsic spatial average in the (๐ฅ, ๐ง)-plane,
โซ
โจ๐โฉ = 1/๐ด ๐ ๐ด ๐๐๐ด (where ๐ด ๐ is the area occupied by fluid), ๐ is the temporal average, ๐โฒ = ๐ โ ๐
๐
is the instantaneous turbulent fluctuation, and ๐e = ๐ โ โจ๐โฉ is the form-induced fluctuation. The area
averaging carried out in the total area of fluid and solid, ๐ด๐ , is termed superficial area averaging,
โซ
denoted by โจ๐โฉ๐ = 1/๐ด๐ ๐ด ๐๐๐ด; the two averaging approaches satisfy the relation โจโฉ๐ = ๐ (๐ฆ)โจโฉ,
๐
where ๐ (๐ฆ) is the plane-averaged porosity at elevation ๐ฆ,
๐ด ๐ (๐ฆ)
๐ (๐ฆ) = . (2.4)
๐ด๐
The plane-averaged total drag ๐น๐ท (๐ฆ), exerted by the grains on the fluid (including both viscous
and pressure drag contributions), is calculated using the IBM body force, ๐น1 :
โซ
๐
๐น๐ท (๐ฆ) = ๐น1 (๐ฅ, ๐ฆ, ๐ง)๐๐ฅ๐๐ง, (2.5)
๐ฟ๐ฅ ๐ฟ ๐ง ๐ด๐
where ๐ฟ ๐ฅ๐ is the domain size in ๐ฅ๐ . For a detailed explanation of this method, see Yuan and Piomelli
(2014b). The friction velocity, ๐ข ๐ , is calculated based on the maximum value of the total stress,
which are the sum of the the viscous shear stress, turbulent shear stress and the form-induced shear
stress.
2.2.2 Parameters
The detailed simulation parameters for all cases are listed in Table 2.1. Two types of SWI โ regular
and random interfaces โ are simulated; their definitions are given in Section 2.2.3. Case 1 is used
in Section 2.3.1 to validate the synthesis of the random interface generated with experimental data.
13
Cases 2 and 3 share a higher friction Reynolds number than that in Case 1; these two cases are used
to compare the effects of different interface geometries on the turbulent flow in Section 2.3.2.
In Table 2.1, ๐ ๐๐ฃ๐ is the volume-averaged porosity inside the sediment; its value is kept the same
for all cases. The bulk permeability, ๐พ, is estimated using the Kozeny-Carman model (Kozeny,
1927; C. Carman, 1937),
3
๐ ๐๐ฃ๐
๐พ= ๐ท 2, (2.6)
180(1 โ ๐ ๐๐ฃ๐ )2
where ๐ท is the grain diameter. The Kozeny-Carman model has been widely used for multiple types
of artificially-generated packed beds, for example, by Voermans et al. (2017) for a random packing
of spheres and by Fang et al. (2018) for periodic arrays of spheres. The elevation of ๐ฆ = 0 is chosen
at the sediment crest, for ease of simulation domain setup. However, ๐ฆ = 0 is not used as the virtual
origin of the ๐ฆ axis for comparison of flow statistics. Instead, the virtual origin is chosen at the
zero-plane displacement, โ๐, which is obtained by fitting the DA velocity to the logarithmic law
(shown in Section 2.3.2). ๐ป๐ is the sediment depth measured from the sediment crest to the bottom
boundary of the simulation domain. Figure 2.1(a,b) show the simulation domain and the relation
between various lengths and the ๐ฆ axis.
The domain sizes for Case 1 in the ๐ฅ and ๐ง directions are 12๐ฟ and 6๐ฟ, respectively. For
Cases 2 and 3 with higher ๐
๐ ๐ values, smaller domain sizes are used. Such domain sizes are
considered sufficiently large for the following reasons. (1) The streamwise and spanwise two-point
autocorrelations of the streamwise turbulent velocity fluctuations, calculated at ๐ฆ โค 0.5๐ฟ, fall below
Figure 2.1: (a) Simulation domain and (b) positions of sediment crest, ๐ ๐ , and zero-plane
displacement, ๐, defined in Section 2.3.2.1.
14
0.1 at half the domain length or width for all cases. (2) The domain is sufficient for the sediment
domains to contain 2000 โ 2500 randomly distributed grains to produce statistically converged
single-point and two-point turbulent flow statistics. A test is performed for Case 3 with a domain
size twice as large in both ๐ฅ and ๐ง, which yields similar results. (3) The sediment depth, ๐ป๐ , is
much larger than the penetration depth of turbulence into the sediment, as shown in Section 2.3.2.
The total simulation time for each case is at least ๐ = 25๐ฟ/๐ข ๐ after the transient.
The number of grid points are 1024(๐ฅ) ร 380(๐ฆ) ร 512(๐ง) for Case 1 and 1536(๐ฅ) ร 478(๐ฆ) ร
768(๐ง) for Cases 2 and 3. For Cases 2 and 3, we have verified that coarsening ฮ๐ฅ๐+ two-fold does
not noticeably affect the results shown in this paper. Here, superscript + indicates normalized by
viscous length scale, ๐ฟ ๐ = ๐/๐ข ๐ . Uniform spacing of grids are used in ๐ฅ and ๐ง; ฮ๐ฅ = ฮ๐ง. The grain
geometry is resolved by 36 grid points along each direction in Case 1 and 50 grid points in Cases 2
and 3. In the ๐ฆ direction, 120 grids are used to refine the grid within one diameter distance below
the crest; deeper inside the sediment, uniform ฮ๐ฆ same as the grid size in ๐ฅ and ๐ง is used. Above
the crest, the ๐ฆ grid is stretched with finer resolution close to the interface. The minimum vertical
grid spaces ฮ๐ฆ +min in wall units is 1.5 for Case 1 and 0.19 for Cases 2 and 3. The maximum ฮ๐ฆ + is
5.3 for Case 1 and 4.2 for Cases 2 and 3.
To provide further evidence that the resolution of 50 grid points per diameter captures well
the velocity distribution and forces produced by the spheres, we conducted another simulation of a
uniform flow past a single sphere with a set-up similar to Mittal et al. (2008) and compared results
to the experimental measurements of Taneda (1956). Two values of the Reynolds number based
on grain diameter, ๐ท๐/๐ (where ๐ is the uniform undisturbed velocity), of 75 and 100 are used.
These values are similar to or higher than the Reynolds number experienced by the spheres near
the interface in this work. Results (not shown) quantify the difference between our DNS and the
experimental results as 1-4% for the flow characteristics (separation angle, recirculation bubble size
and recirculation center location) and 1 โ 1.5% for the drag force. The comparison gives confidence
that the results herein, especially the local mean shear layers (important for turbulence production)
and the drag force (important for form-induced-stress production (Raupach and Shaw, 1982)) are
15
well captured by the spatial resolution.
2.2.3 Synthesis of porous beds
Numerical studies of permeable walls usually employ idealized models characterized by regular
packing of grains and a constant bed height due to their easy implementation. Examples include a
Cartesian grid of cubes used by Breugem and Boersma (2005), a regular packed bed of Breugem
et al. (2006), or a simple cubic packing of spheres used in Stoesser et al. (2007) and Fang et al.
(2018). Numerical approximations of natural sediments were constructed by, for example, Stubbs
et al. (2018), who used a CAD model to design and manufacture an artificial gravel riverbed
to approximate a natural bed for both experimental and numerical studies. Here, we implement
a different approach based on molecular dynamics (MD) simulations, which are widely applied
in flows through porous media with applications in materials science (Khirevich, 2011; Amadio,
2014). The simulations are carried out using the open-source code LAMMPS (Plimpton, 1995)
to generate sediments composed of randomly packed, monodisperse hard spheres. The process of
pouring hard spheres into a tank is simulated; this process is designed to reproduce the sediment
used in the experimental studies of Voermans et al. (2017) (henceforth referred to as VGI17).
Figure 2.2: Synthesized sediment beds colored by ๐ฆ/๐ฟ with (a) the regular and (b) the random
interfaces.
16
Details of the MD simulations are included in Appendix A.
Two types of distribution of the sphere-center locations of the uppermost-layer spheres are
generated: (1) regular distribution in (๐ฅ, ๐ง) with constant height, ๐ฆ, called a regular interface, and
(2) random values of both (๐ฅ, ๐ง) and ๐ฆ directions, called a random interface. Below the interface,
sediments with both types of interface are synthesized using a random sphere distribution that is
statistically similar. The total number of spheres are 2633, 2351 and 2116, for Cases 1, 2, and 3,
separately. Figure 2.2 compares the sphere distributions in Cases 2 and 3; the difference at the
top of the sediments is apparent. These characteristics of the interface can also be expressed in
terms of the mean and variance of the permeability in the uppermost sediment layer as well as the
correlation length scales in the ๐ฅ and ๐ง directions (Section 2.2.4).
2.2.4 Geometrical comparison of the two interface types
The characteristics of the grain distributions in Cases 2 and 3 are compared. First, the plane-
averaged porosity profiles along the vertical direction are shown in Figure 2.3. In the interface
region (๐ฆ โ 0), the random interface yields monotonically decreasing porosity with decreasing
๐ฆ, with an inflection point as observed in VGI17. In comparison, the regular sediment yields a
Figure 2.3: Wall-normal profiles of plane-averaged porosity profile for regular ( ) and random
( ) interfaces.
17
different pattern, with large-amplitude fluctuations. This is due to the constant sphere height in the
uppermost layer. The vertically-averaged porosity inside the sediment is the same between the two
cases.
To characterize the positioning of the sediment grains located in the vicinity of the interface,
we define the local bed surface height, โ(๐ฅ, ๐ง), as the highest sphere-surface elevation at each (๐ฅ, ๐ง)
location. Note that this is different from using the sediment-grain center height (which would
yield constant height distribution for the regular interface). Since the penetration depths are within
one sphere diameter as shown later in Table 2.2, only sphere surfaces within one sphere diameter
distance below the crest are considered in calculating โ(๐ฅ, ๐ง). The number of top-layer spheres
identified in this way are 900 for the regular interface (Case 2) and 632 for the random interface
(Case 3). The sediment-height fluctuations, โโฒ (๐ฅ, ๐ง), are defined as the deviations of the local height
values โ(๐ฅ, ๐ง) from its plane-averaged value. Figure 2.4(a,b) compare โโฒ (๐ฅ, ๐ง) distributions in Cases
2 and 3, normalized by ๐ฟ. Interpolation is carried out to show a smooth variation for demonstration.
The magnitude of โ(๐ฅ, ๐ง) fluctuations can be characterized by the root-mean-square value of โโฒ, ๐โ .
The ๐โ+ values are 17 and 23 for Cases 2 and 3, respectively. The horizontal distribution of โโฒ (๐ฅ, ๐ง)
15 28
100
10-1
10-2
101 102 103
Figure 2.4: Bed surface height fluctuations for the (a) regular and (b) random interfaces. (c)
Autocorrelation of height fluctuations and Taylor micro scales.
18
is characterized here by a Taylor micro scale, ๐, determined from a polynomial fit of the two-point
autocorrelation ๐
๐ฅ of โโฒ, similar to the method used by Padhi et al. (2018a). ๐
๐ฅ is calculated as
โซ
1
๐
๐ฅ = โโฒ (๐ฅ, ๐ง)โโฒ (๐ฅ + ๐ ๐ฅ , ๐ง)๐๐ด, (2.7)
๐ด๐ ๐ฅ,๐ง
where ๐ ๐ฅ is the separation length in ๐ฅ. Figure 2.4(c) shows that ๐+ , estimated using the two-point
correlation with streamwise separation, is 15 and 28 for Cases 2 and 3, respectively. The significant
difference in roughness length scales affects the apparent permeability in the uppermost sediment
layer, as shown in Section 2.3.3.1. The random interface is more heterogeneous with a larger mean
and variance of the apparent permeability, although the permeability in the bulk of the sediment is
the same.
The probability density functions of โโฒ (๐ฅ, ๐ง) are shown in Figure 2.5. Case 2 (regular) gives a
highly skewed probability distribution with a skewnwss of -2.5 and a kurtosis of 8.7, while Case
3 (random) displays a less skewed but still asymmetric distribution with a skewness of -0.075 and
kurtosis of 1.9. The skewness and a kurtosis of Case 3 are well within the ranges for gravel-beds
summarized in Nikora et al. (1998), confirming that the interface is random. These observations
show that the major differences between the sediment surfaces in Cases 2 and 3 are the larger
wall-normal and horizontal length scale for the random case. These lengths are small compared
to the grain size (๐ท + โ 70); thus, such height fluctuations are considered as roughness, instead
Figure 2.5: Histogram of bed height fluctuations for the (a) regular and (b) random interfaces.
19
of bedform, whose length scales are much larger than the grain scale. We have also verified
that different realizations of the random interface (obtained from separate MD simulations) share
similar geometrical characteristics of the sediment surface. For example, another realization of the
sediment in Case 1 gives matching values of ๐โ /๐ท and ๐/๐ท within 1% difference.
2.3 Results
2.3.1 Validation of sediment synthesis
To test how the random interface mimic sediment beds in natural and laboratory settings, we
compare Case 1 with the measurements of Case L12 from VGI17. The present values of ๐ = 0.41
and ๐
๐ ๐พ = 2.56 match those in VGI17 and fall in the range in real aquatic systems. To be consistent
with VGI17, in this section we define the location of ๐ฆ = 0 to be the location where ๐ 2 ๐/๐๐ฆ 2 = 0.
Note that for the rest of this paper, ๐ฆ = 0 is defined at the sediment crest (while the virtual origin is
chosen as ๐ฆ = โ๐ for statistics comparison, as discussed in Section 2.3.2).
The comparison of the mean velocity profile ๐ = โจ๐ขโฉ normalized by the DA free-surface (or
channel-half-height) flow velocity, ๐๐ฟ , is shown in Figure 2.6(a). Excellent agreement is obtained.
1 1 1 1
0.5 0.5 0.5 0.5
0 0 0 0
-0.5 -0.5 -0.5 -0.5
0 0.5 1 0 1 2 0 0.5 1 -1 -0.5 0
Figure 2.6: Comparison of (a) mean velocity and (b-d) streamwise, wall-normal and shear
components of the Reynolds stress tensor. โฆ Experiment, DNS.
20
Figure 2.6(b-d) show the comparison of the turbulence intensities normalized by ๐ข ๐ . Very good
agreement is obtained and the slight differences in the outer layer (๐ฆ/๐ฟ โฅ 0.3) are probably due
to the fact that the flow in the experiment at the measurement station is not a fully developed
open-channel flow, but resembles a spatially-developing boundary-layer flow.
The intensities of the form-induced fluctuations normalized by ๐ข ๐ are compared in Figure 2.7.
Relatively large discrepancies are found between the present DNS and the experimental measure-
ments. This is at least partially attributed to a difference in the sampling size. The sizes of the
sediment domains are similar between the present DNS and the experiment; however, for the DNS,
the spatial averaging is carried out throughout the entire (๐ฅ, ๐ง) domain, while the spatial averaging
in the experiment was carried out over six lateral measurement frame positions at one streamwise
location only. The form-induced stresses have been found highly sensitive to the interface geom-
etry (Nikora et al., 2001; Fang et al., 2018). If we mimic the experimental sampling protocol by
applying the spanwise averaging procedure at six uncorrelated spanwise locations and then repeat
the procedure for different streamwise locations using the DNS data, we obtain a family of curves
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
-0.2 -0.2 -0.2
-0.4 -0.4 -0.4
0 1 2 0 0.5 -0.4 -0.2 0
Figure 2.7: Comparison of (a) streamwise, (b) wall-normal and (c) shear component of the
form-induced stress tensor. โฆ Experiment, DNS. Black lines represent the values spatially
averaged over the entire horizontal plane; grey lines emulate the sampling used in the experiment.
21
Interface ๐
๐/๐ฟ ๐ถ๐ ๐ฟ ๐ /๐ฟ ๐ฟโ๐ /๐ฟ ๐ฟ ๐ /๐ฟ ๐ฟโ๐ /๐ฟ ๐ln(๐พ โ ) 2
๐ln(๐พ โ)
Case 2 Regular 0.32 0.06 0.008 0.078 0.138 0.097 0.157 -6.472 0.005
Case 3 Random 0.33 0.11 0.013 0.055 0.165 0.117 0.227 -5.391 0.055
Table 2.2: Flow parameters for Cases 2 and 3. ๐ฟ ๐ and ๐ฟ ๐ are Brinkman-layer thickness and
Reynolds-stress penetration depth, both measured from ๐ฆ = โ๐; ๐ฟโ๐ and ๐ฟโ๐ are the same lengths
measured from crest; ๐พ โ is apparent permeability.
6
15
5
4
10
3
2 5
1
0 0
0 100 200 300 100 101 102
Figure 2.8: (a) Diagnostic function for fitting the logarithmic profile and (b) DA velocity profiles
for regular ( ) and random ( ) interfaces; (red) fitted logarithmic profiles.
whose envelope provides an indication of the level of uncertainty in the experimental measure-
ments. Such a family of curves are shown by the gray lines in Figure 2.7 while the black solid lines
denotes DNS results using spatial averaging over the entire domain. The VGI17 data points fall
within the scatter.
In addition to the uncertainties arising due to sampling size, subtle differences in details of the
sphere distribution near the interface may also contribute to discrepancies in form-induced velocity
fluctuations. However, since the detailed grain-distribution characteristics are not available from
VGI17, we consider the synthesized sediment sufficiently โrealisticโ as it reproduces quantitatively
the turbulence statistics and matches qualitatively the form-induced fluctuations. In comparison, a
sediment with a regular interface may yield flow statistics that are drastically different as described
in Sec. 2.3.2.
22
2.3.2 Effects of interface irregularity on flow statistics
In this section, Cases 2 and 3 are compared to understand the effects of interface characteristics on
turbulence statistics and structure.
2.3.2.1 DA velocity and friction
The location of the zero-plane displacement, ๐, and the von Kรกrmรกn constant, ๐
, of the wall-bounded
flow are determined by fitting the streamwise DA velocity profile to the logarithmic law,
1
๐+ = log(๐ฆ + ๐) + + ๐ต, (2.8)
๐
where ๐ต is the intercept of the logarithmic profile. The fitting procedure is similar to the method
used by Breugem et al. (2006) and Suga et al. (2010) and is demonstrated in Figure 2.8(a). The fitted
velocity profiles are shown in Figure 2.8(b). The fitted values of ๐
and ๐ are listed in Table 2.2; these
values are close to the results (๐
= 0.31 โ 0.32 and ๐/๐ฟ = 0.06 โ 0.1) reported by Breugem et al.
(2006), Suga et al. (2010) and Manes et al. (2011) at similar ๐
๐ ๐พ . It is well established that ๐
for
turbulence bounded by a porous wall is smaller than the corresponding value for an impermeable
0.02
0.8 0
-0.02
0.6 -0.04
-0.06
0.4 0 0.04
0.2
0
-0.2
0 0.5 1
Figure 2.9: DA velocity along ๐ฆ in linear scaling for regular ( ) and random ( )
interfaces. Inset magnifies the distribution at ๐ฆ โ โ๐.
23
wall, which is usually within 5% of 0.41 for both smooth and rough walls. Despite sharing similar
๐
values, Cases 2 and 3 differ significantly in ๐. Recalling that ๐ is measured from the sediment
crest height, the higher ๐ in Case 3 is probably due to the larger variance in local sediment heights
associated with the random interface. In the following plots we offset ๐ฆ by the amount of ๐ to
effectively collapse the logarithmic regions between the two cases.
The DA velocity profiles in linear scaling (Figure 2.9) show the differences in the velocities
near the interface. The velocity in the vicinity of the regular interface is negative with a small
magnitude. This can be related to the organized recirculation regions induced by each sphere at
the top layer (discussed in Section 2.3.2.3). This is an important observation with implications for
interfacial exchange of solutes. Also, simpler models, especially those that ignore inertial effects,
may not be able to capture these features. Such a flow pattern is also described in the conceptual
model of Pokrajac and Manes (2009). For lower ๐ฆ, the velocity decreases and eventually reaches
a constant, positive value inside the sediment. For the random interface, the velocity variation is
monotonic. Both velocity profiles exhibit an inflection point near the sediment crest, consistent
with observations of Manes et al. (2011); Voermans et al. (2017).
The depth of shear-layer penetration, ๐ฟ ๐ , is defined as the distance from the ๐ฆ = โ๐ to the ๐ฆ
location separating the constant-velocity region in the bed and the shear layer above. ๐ฟ ๐ is also
referred to as the Brinkman-layer thickness. Here, ๐ฟ ๐ is calculated in the same way as Voermans
et al. (2017), โจ๐ขโฉ๐ฆ+๐=โ๐ฟ๐ = 0.01(๐๐ โ ๐ ๐ ) + ๐ ๐ , where ๐๐ is the DA velocity at ๐ฆ = โ๐ and ๐ ๐
is the constant velocity deep in the sediment. ๐ฟ ๐ is the penetration depth measured from ๐ฆ = โ๐.
The value measured from the crest, ๐ฟโ๐ = ๐ฟ ๐ + ๐, for the two cases are shown in Table 2.2, equaling
0.138๐ฟ (regular case) and 0.165๐ฟ (random case). They are of the same order of the grain diameter,
as also observed by Goharzadeh et al. (2005).
The friction coefficient is defined as
2
๐ข๐
๐ถ๐ = 2 . (2.9)
๐๐ฟ
Table 2.2 shows that the friction coefficient on the random interface is 70% higher than the regular
one. A higher ๐ถ ๐ for the random interface is expected as the larger interface length scales (๐โ+ and
24
๐+ ) lead to higher total drag, similar to the effects of an impermeable wall with a larger roughness
length scale.
0.8 0.8
0.6 0.6
0.4
0.4
0.2
0.2
0
0
-0.2
0 2 4 6 -0.4 -0.2 0 0.2 0.4
Figure 2.10: Components of (a) the Reynolds-stress tensor and (b) the anisotropy tensor for
regular ( ) and random ( ) interfaces.
2.3.2.2 Turbulent fluctuations
Figures 2.10(a,b) compare the Reynolds stress tensor and its anisotropy between the two cases.
The Reynolds-stress anisotropy tensor is defined as
โจ๐ขโฒ๐ ๐ขโฒ๐ โฉ ๐ฟ๐ ๐
๐๐ ๐ = โ . (2.10)
โจ๐ขโฒ๐ ๐ขโฒ๐ โฉ 3
It is shown that Townsendโs wall-similarity hypothesis for a wall-bounded flow applies to the outer
layer (๐ฆ/๐ฟ > 0.2). As ๐ฆ decreases from the sediment crest, for both cases all components of the
Reynolds stress tensor are damped rapidly, together with a decrease of tensor anisotropy. In the
interface region, the random interface results in a lower Reynolds stress anisotropy with a higher
fraction of TKE residing in ๐ฃ โฒ and ๐ค โฒ motions.
The more isotropic turbulence near the random interface is linked to the augmented disturbance
from the sediment obstruction to the near-wall turbulent coherent structures. Figure 2.11 compares
25
the contours of ๐ขโฒ+ at the respective locations of the โจ๐ขโฒ2 โฉ peaks in the two cases. The low-speed
streaks near the random interface are significantly shorter in wall units than those near the regular
interface due to the physical blockage from the spheres serving as roughness elements.
Figure 2.11: Contours of instantaneous ๐ขโฒ+ in the (๐ฅ, ๐ง) plane at the respective โจ๐ขโฒ2 โฉ-peak
elevations for the (a) regular and (b) random interfaces.
The penetration depth of turbulent shear stress into the bed, ๐ฟ ๐ , is obtained from โจ๐ขโฒ๐ฃ โฒโฉ๐ฆ=โ๐ฟ ๐ =
0.01โจ๐ขโฒ๐ฃ โฒโฉ๐ , following Voermans et al. (2017), and shown in Table 2.2. โจ๐ขโฒ๐ฃ โฒโฉ๐ is the value of โจ๐ขโฒ๐ฃ โฒโฉ
at ๐ฆ = โ๐. ๐ฟโ๐ is the same depth measured from the crest. The flow over the random interface
gives significantly deeper turbulence penetration, as shown by 47% higher ๐ฟโ๐ and 20% higher in
๐ฟ ๐ compared to the regular interface.
To quantitatively compare turbulent scales, the integral length scales are shown in Figure 2.12.
300
200
100
0
0 0.1 0.2 0.3 0.4 0.5 0.6
Figure 2.12: Integral length scales of ๐ขโฒ motions along the streamwise directions for regular
( ) and random ( ) interfaces.
26
It is defined as
โซ โ
๐ฟ ๐ ๐,๐ฅ ๐ (๐ฆ) = ๐
๐ ๐ (๐ฆ, โณ๐ ๐ )๐ (โณ๐ ๐ ), (2.11)
0
where
โจ๐ขโฒ๐ (๐ฆ, ๐ฅ ๐ )๐ขโฒ๐ (๐ฆ, ๐ฅ ๐ + โณ๐ ๐ )โฉ
๐
๐ ๐ (๐ฆ, โณ๐ ๐ ) = , (2.12)
๐๐ข๐ (๐ฆ)๐๐ข ๐ (๐ฆ)
(no summation over repeated indices). โณ๐ ๐ is the separations along the ๐ฅ ๐ direction and ๐๐ข๐ (๐ฆ),
๐๐ข ๐ (๐ฆ) are the root-mean-square of turbulent fluctuations. The integration is carried out to the
value of โณ๐ ๐ at which the correlation coefficient first crosses 0.3. A threshold value between 0.3
and 0.5 is usually used to calculate the integral lengths from two-point correlation coefficients
(see, for examples, Krogstad and Antonia (1994), Christensen and Wu (2005) and Volino et al.
(2011)). It has been checked that varying this threshold from 0.2 to 0.5 would not noticeably affect
the comparison. It is shown that, for the random interface, at ๐ฆ = โ๐ the coherent motions are
more extensive in ๐ฅ due to a deeper flow penetration. At the sediment crest, however, the coherent
motions are noticeably shorter in ๐ฅ for the random interface, due to the disturbances of the larger
roughness height.
Farther away from the interface ((๐ฆ + ๐)/๐ฟ > 0.15), the difference in the ๐ฅ integral length
between the two cases systematically increases with ๐ฆ, indicating a lack of Townsendโs wall-
similarity hypothesis for this flow quantity. One reason is the relatively high roughness. Jimรฉnez
(2004) has proposed that Townsendโs hypothesis applies when the roughness height is less than
0.02๐ฟ. In the present study, the roughness height (measured from ๐ฆ = โ๐) is at least 0.06๐ฟ. Another
possible reason is the limited Reynolds number used in the present DNS. The experimental study
of Krogstad and Antonia (1994) showed a lack of similarity in the outer-layer ๐ฟ 11 in ๐ฅ between a
smooth-wall and a rough-wall turbulent boundary layer flow. However, a repeated experiment of
boundary layer flows by Krogstad and Efros (2012) at a higher Reynolds number indicated reduced
discrepancy in the outer-layer ๐ฟ 11 .
27
2.3.2.3 Form-induced fluctuations
The form-induced stresses are shown in Figure 2.13. The peak values of all components are much
smaller than the corresponding Reynolds stresses, consistent with observations of Manes et al.
(2009), Voermans et al. (2017) and Fang et al. (2018). The streamwise form-induced stress near
the regular interface region is concentrated in a narrower layer near the crest elevation with a much
higher peak value. This is due to the vertical alignment of the wake regions of the grains as shown
in Figures 2.14(b,c). Large values of โจe ๐ข 2 โฉ are found in the troughs of the top sediment surface and
in the wake regions of the surface protuberances. Since these regions are located in a narrow layer
๐ข 2 โฉ is narrower with a larger peak. In contrast, Figure 2.15 shows
near the crest, the distribution of โจe
that, for the random interface the mean recirculation regions can reach much larger sizes and that
the penetration of the time-mean flow is much deeper. In addition, โจe ๐ฃ 2 โฉ, โจe
๐ค 2 โฉ and โจe
๐ขe๐ฃ โฉ are higher
near the random interface. This may be explained by the higher spatial variation of e ๐ข associated
0.8 0.8 0.8
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
-0.2 -0.2 -0.2
-0.4 -0.4 -0.4
0 2 4 -0.1 0 0.1 0.2 0 0.2 0.4
Figure 2.13: Components of the form-induced stress tensor: (a) streamwise, (b) wall-normal and
(c) spanwise components, for regular ( ) and random ( ) interfaces.
28
with the disordered blockage, which leads to more intense e ๐ฃ and ๐คe through continuity.
Figure 2.14: (a) Time-averaged streamlines for the regular interface. (b) Streamlines on Slice A
(๐ฅ โ ๐ฆ plane across the crests of the grains). (c) Streamlines distribution on Slice B (๐ฅ โ ๐ฆ plane
across the valley between two rows of grains). Dash lines indicate ๐ฆ = โ๐.
To analyze the importance of form-induced stresses to the mean momentum balance, we
compare the magnitudes of the various shear stresses in Figure 2.16. For the regular interface, the
form-induced shear stress is negligible; the viscous shear stress contributes significantly to the total
shear stress at the sediment crest. However, for the random interface the viscous shear is weaker
due to the milder DA velocity gradient associated with the thicker shear layer; a much stronger
Figure 2.15: (a) Time-averaged streamlines for the random interface. (b) and (c) show streamlines
on two ๐ฅ โ ๐ฆ planes, Slice C and Slice D. Dash lines indicate ๐ฆ = โ๐.
29
form-induced stress is present.
0 0
-0.2 -0.2
-0.4 -0.4
-0.6 -0.6
-0.8 -0.8
-1 -1
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Figure 2.16: Shear stresses for cases with (a) regular and (b) random interfaces: Reynolds
shear stress, viscous shear stress and (red) form-induced shear stress.
It is important to note that, in the vicinity of the random interface, the form-induced shear
stress reaches a similar magnitude of the Reynolds shear stress (Figure 2.16(b)). The fact that
the form-induced stress potentially affects the DA velocity near the interface and that it strongly
depends on the interface geometry calls for a systematic study of the generation of the form-induced
stresses; however, this is beyond the scope of the present work.
2.3.2.4 Reynolds-stress budgets
The budget equations of the normal Reynolds stresses, โจ๐ขโฒ2 ๐ผ โฉ๐ , can be written as (Raupach and Shaw,
1982; Mignot et al., 2009; Yuan and Jouybari, 2018)
โฒ โฒ
๐โจ๐ข ๐ผ โฉ โฒ โฒ ๐e๐ข๐ผ โฒ โฒ ๐e๐ข๐ผ ๐ ย โฒ โฒ
0= โ2โจ๐ข ๐ผ ๐ฃ โฉ๐ โ2โจ๐ข ๐ผ ๐ข ๐ โฉ๐ โ2 ๐ข ๐ผ ๐ข ๐
ย โ ๐ข ๐ข e ๐ข๐ (2.13)
๐๐ฆ ๐๐ฅ ๐ ๐๐ฅ ๐ ๐ ๐๐ฅ ๐ ๐ผ ๐ผ ๐
| {z }| {z }| {z }| {z }
๐๐ ๐๐ ๐๐ค ๐๐ค
* + * + * +
๐ 2 ๐ขโฒ2 โฒ โฒ ๐๐ขโฒ
๐ โฒ โฒ โฒ ๐๐ ๐๐ข
โ ๐ข ๐ข ๐ข +๐ ๐ผ
โ2 ๐ขโฒ๐ผ โ2๐ ๐ผ ๐ผ
.
๐๐ฅ ๐ ๐ผ ๐ผ ๐ ๐ ๐๐ฅ ๐ ๐๐ฅ ๐ ๐๐ฅ ๐ผ ๐๐ฅ ๐ ๐๐ฅ ๐
| {z }| ๐ ๐ ๐
{z } | {z } | {z }
๐๐ก ๐๐ ฮ ๐
On the right-hand side, the first three terms are, respectively, shear production (๐๐ ) and additional
productions due to the form-induced shear (๐๐ค and ๐๐ ). The following three terms are the transport
30
terms due to the form-induced (wake) fluctuations (๐๐ค ), turbulent fluctuations (๐๐ก ), and viscous
diffusion (๐๐ ). ฮ is the pressure work and ๐ is the viscous dissipation. Summing over the equations
of the three Reynolds-stress components yields the budget equation of the TKE; the sum of the ฮ
terms yields the pressure transport, ๐๐ . The production terms and viscous dissipation of the TKE
are shown in Figure 2.17(a), whereas the transport terms are shown in Figure 2.17(b). The residuals
are around 3% of ๐๐ peak values.
Figure 2.17: Terms in the TKE budget equation. (a) Production and viscous dissipation for
regular (symbols and thin lines) and random-interface (thick lines) cases; ๐๐ ( , โณ) , ๐๐ + ๐๐ค
( , โ, red) and ๐ ( , โก). Transport terms for (b) regular and (c) random interfaces; ๐๐ค
( , red), ๐๐ก ( ), ๐๐ ( ) and ๐๐ ( ). All terms are normalized by ๐ข 3๐ /๐ฟ.
For (๐ฆ + ๐)/๐ฟ โฅ 0.2, the shear production is balanced by the viscous dissipation (not shown). ๐๐
terms reach their maximum values close to the sediment crest. The total amounts of form-induced
production (๐๐ +๐๐ค ) are significant near the interface, becoming the dominant term a short distance
below the crest. Fang et al. (2018) also observed non-negligible form-induced TKE production for
regular interface particle distribution with ๐
๐ ๐พ = 0 โ 27.
In the vicinity of the regular interface, ๐๐ , ๐๐ก and ๐๐ remove TKE from the high-TKE region
slightly above the crest and transfer it into the low-TKE region in the bed; ๐๐ dominates other TKE
transports due to high ๐๐/๐๐ฆ, while ๐๐ค is negligible. In contrast, for the random interface ๐๐ is
negligible and the augmented e ๐ฃ fluctuations (Figure 2.13) lead to a non-negligible form-induced
transport, which works against the other transport processes by moving TKE upward from low-TKE
31
region in the bed to the crest region corresponding to the TKE peak. Such process contributes to a
more equilibrium turbulence (with more energy dissipated at the location where it is generated).
The total form-induced production (๐๐ค + ๐๐ ) for individual normal Reynolds stresses are
compared in Figure 2.18(a). For both cases below the crest, it predominantly contributes to โจ๐ขโฒ2 โฉ
other than the other two components, more so for the regular interface. ๐๐ can be rewritten as
โ2โจ๐ขโฒ๐ผ ๐ฃ โฒโฉโจ๐ข ๐ผ โฉ๐๐/๐๐ฆ. Thus, ๐๐ remains positive for โจ๐ขโฒ2 โฉ and is zero for โจ๐ฃ โฒ2 โฉ and โจ๐ค โฒ2 โฉ. The
negative ๐๐ค + ๐๐ in both cases in the vicinity of the crest is then due to negative ๐๐ค , indicating
conversion of TKE to the kinetic energy of form-induced fluctuations, since this term also exists in
the form-induced stress budgets with a positive sign (Raupach and Shaw, 1982). Such conversion
is probably due to the work of particle drag against the turbulent motions of scales larger than the
particles, generating form-induced fluctuations at the scale of the particles. Deeper into the porous
bed, the scales of turbulent motions are rapidly reduced and are smaller than that of the particles
(as shown by the integral scales in Figure 2.12); here, the particle-scale form-induced shear layers
generates turbulence, represented by positive ๐๐ค .
The pressure work for a normal Reynolds stress is decomposed into the pressure transport and
60 60
40 40
20 20
0 0
-20 -20
-40 -40
-60 -60
-0.05 0 0.05 0.1 0.15 0.2 -0.05 0 0.05 0.1 0.15 0.2
Figure 2.18: (a) Form-induced production and (b) pressure-strain-rate term in regular ( ) and
random interface ( โฒ2 โฒ2 โฒ2
) cases for the budgets of โจ๐ข โฉ๐ (โ, red), โจ๐ฃ โฉ๐ (โก) and โจ๐ค โฉ๐ (โณ, blue);
normalization is done with ๐ข ๐ and ๐ฟ.
32
a pressure-strain-rate term, R = 2โจ๐โฒ ๐๐ขโฒ๐ผ /๐๐ฅ ๐ผ โฉ๐ . The latter is compared in Figure 2.18(b). R
distributes energy among different normal components of the Reynolds stress tensor; it forms the
primary source for ๐ฃ โฒ- and ๐ค โฒ- energy budgets that balances the viscous dissipation. The main
difference is that, at the crest of the regular interface the ๐ขโฒ energy converts mostly to that of ๐ค โฒ,
with R 22 โ 0, while R 22 and R 33 are of comparable magnitudes at the crest for the random case.
Liu and Katz (2018) observed experimentally that, in a shear layer formed over an open cavity,
the magnitude of R 11 is an order of magnitude higher than R 22 . The similar observation herein
may also be attributed to the local mean shear layers formed above the topmost-layer particles,
shown in Figure 2.19 using time-mean spanwise vorticity in an (๐ฅ, ๐ฆ) plane. Such shear layers are
stronger and aligned in ๐ฆ above the regular interface, while the random distribution of particle for
the random interface results in a spread in such effect on energy redistribution.
Figure 2.19: Time-averaged spanwise vorticity ๐ ๐ง normalized by ๐ข ๐ /๐ฟ at (a) Slice B of the regular
interface and (b) Slice D of the random interface.
33
2.3.3 Implications for the mechanism of flat-bed hyporheic exchange
2.3.3.1 Variation of apparent permeability
As a convention of the oil and gas industry (Edwards et al., 1990; Huang et al., 2006; Javadpour
et al., 2009), the apparent permeability, ๐พ โ , is used to establish a relationship between mean velocity
and macroscopic pressure gradient in a finite Reynolds number flow inside a porous medium. Here,
the local value of ๐พ โ is obtained based on the local volume-averaged velocity and the gradient of
the local volume-averaged pressure,
๐พ โ (๐ฅ, ๐ฆ, ๐ง) ๐โจ๐โฉ(๐ฅ, ๐ฆ, ๐ง)
โจ๐ขโฉ ๐ (๐ฅ, ๐ฆ, ๐ง) = โ , (2.14)
๐ ๐๐ฅ
where โจยทโฉ๐ and โจยทโฉ indicates, respectively, the superficial and intrinsic volume averaging performed
locally (within a volume much larger than a grain size). Normalized by ๐ข ๐ and ๐ฟ, Equation (2.14)
becomes !
๐พโ ๐โจ๐โฉ +
โจ๐ขโฉ๐ + = โ๐
๐ ๐ , (2.15)
๐ฟ2 ๐๐ฅ/๐ฟ
keeping in mind that ๐ = ๐/๐. For convenience, from here on write the dimensionless ๐พ โ /๐ฟ2 as
๐พ โ . Note that ๐พ โ is not the Darcy permeability, but a macroscopic descriptor of the complex flow
-3
10
3
6
2
1 4
0 2
0 2 4 6
Figure 2.20: Horizontal variation of the apparent permeability ๐พ โ in the Brinkman layer in the
random-interface case.
34
at the interface, including the influences of both diffusion of the mean momentum (traditionally
modelled using the Brinkman correction (Brinkman, 1949)) and microscopic inertial and drag
forces (traditionally modelled using the Forchheimerโs correction (Hassanizadeh and Gray, 1987).
)
Since the Brinkman layer (measured from crest) has a thickness (๐ฟโ๐ ) of the order of a sphere
diameter, ๐พ โ in this layer is considered constant along ๐ฆ, but varying in ๐ฅ and ๐ง, in this layer. To
calculate such ๐พ โ (๐ฅ, ๐ง), we divide the Brinkman layer horizontally into small segments of a size of
2๐ท ร 2๐ท in ๐ฅ and ๐ง. Enlarging the segment size to 4๐ท ร 4๐ท or 8๐ท ร 8๐ท does not fundamentally
change the observation. The velocity and pressure values in Equation (2.14) are spatially averaged
in each small segment respectively; the pressure gradient for each segment is obtained using
second-order central difference. The (๐ฅ, ๐ง) distribution of ๐พ โ for the random interface is shown
in Figure 2.20; the variation for the regular interface is between 1 ร 10โ3 and 2 ร 10โ3 , barely
visible using the same color bar and is thus not shown. The random interface leads to a much more
heterogeneous permeability distribution with a higher mean value. Based on its natural logarithm,
2
the mean, ๐๐๐(๐พ โ ) , and the variance, ๐๐๐(๐พ โ ) , are calculated and summarized in Table 2.2. For the
regular and random interfaces, respectively, the ๐๐๐(๐พ โ ) values are โ6.5 and โ5.4 (or twice higher
in mean apparent permeability for the random case) and the ๐๐๐(๐พ 2
โ ) values are 0.005 and 0.055.
Since the square root of permeability has a physical meaning of the effective pore size (Breugem
et al., 2006), the results here show that the random interface increases the effective pore size at
the interface associated with a larger apparent permeability. This is consistent with the larger
streamwise integral length scales of ๐ขโฒ motions at ๐ฆ โ โ๐, as shown in Figure 2.12.
The ๐พ โ value is also calculated below the Brinkman layer using equation (2.15). Here, ๐พ โ
becomes the bulk permeability, ๐พ, associated with laminar flow inside the porous medium, which
is typically predicted using a semi-empirical model such as the Kozeny-Carman equation. The
value of ๐๐(๐พ โ ) obtained using the DNS data is almost a constant of -10.9, regardless of location.
Such value indicates a ๐พ โ that is almost two order-of-magnitude smaller than the value in the
Brinkman layer. In comparison, the Kozeny-Carman equation predicts a ๐๐(๐พ) value of โ10.0
35
inside the bulk of the sediment, overestimating the bulk permeability by almost 1.5 times.
In the next section, it will be shown that the spatial variation of apparent permeability in the
Brinkman layer of the random interface is correlated with the roughness geometry. As ๐พ โ in the
Brinkman layer incorporates the effect of both macroscopic diffusion and microscopic inertia on the
macroscopic (volume-averaged) velocity, we focus on the spatial variation of these two mechanisms
as affected by the roughness geometry.
2.3.3.2 Momentum transport mechanisms and dependence on roughness geometry
The double-averaged ๐ข-momentum equation can be written as (Raupach and Shaw, 1982)
ยฏ ๐
2
ยฏ ๐ ๐โจ๐ขโฒ๐ฃ โฒโฉ๐
๐โจ๐โฉ ๐ โจ๐ขโฉ ๐โจ๐ขห ๐ฃห โฉ๐
โ + ๐ โ โ + ๐น๐ท = 0, (2.16)
๐๐ฅ ๐๐ฆ 2 ๐๐ฆ ๐๐ฆ
| {z } | {z }
D I
where ๐น๐ท is the total drag per unit mass (sum of both pressure drag and viscous drag around the
grains). According to Whitaker (1996), the viscous and Reynolds shear stress terms represent the
effects due to diffusion of the mean momentum, usually modelled using the Brinkman correction;
Figure 2.21: Terms in the DA ๐ข-momentum equation (equation (2.16)), normalized by ๐ข 2๐ /๐ฟ for (a)
regular and (b) random interfaces: Reynolds-stress term( , red), viscous stress term ( ,
red), total drag ( ), form-induced stress term ( ) and pressure-gradient term( ).
36
the sum of these two terms is denoted herein as D, representing โdiffusionโ. The terms related
to the form-induced fluctuations (i.e., the form-induced shear stress and the total drag) represent
the effects of fluid inertia at the grain scale on the mean momentum, usually modelled by the
Forchheimerโs correction. Their sum is denoted as I, representing (microscopic) โinertiaโ.
3 15 3 15
2 10 2 10
1 5 1 5
0 0 0 0
0 2 4 6 0 2 4 6
3 0.04
0.02
2
0
1 -0.02
0 -0.04
0 2 4 6
Figure 2.22: Horizontal variation of local volume-averaged D (a) and |I| (b), both normalized by
๐ข 2๐ /๐ฟ. (c) sediment surface height fluctuations โโฒ normalized by ๐ฟ for the random-interface case.
(d) Correlations of D โฒ (red) and |I|โฒ (black) with โโฒ.
Figure 2.21 shows the balance of equation (2.16) across the interface, with the spatial averaging
performed using area averaging in the ๐ฅ-๐ง plane of the whole domain. The residual is less than 2.3%
of the peak magnitude of the total force term for the two cases. For both cases, the spatial-averaged
values of D and I are both significant inside the Brinkman layer, approximately balancing each
37
other as the pressure gradient is comparatively weak. The total drag is the main contributor to I
for both cases. For the regular interface, D is contributed almost equally by both Reynolds and
viscous stresses, while it is predominantly due to the Reynolds stress for the random one; this is
consistent with the observations made for shear stress distribution (Figure 2.19). The total drag
peaks at the crest for the regular interface. This is because the main contributor of ๐น๐ท โthe pressure
dragโis caused by recirculation regions, which are distributed uniformly near the crest for this case.
For the random interface, recirculation regions of various sizes are formed due to various scales of
roughness protuberances, leading to a peak of ๐น๐ท below the crest.
We now explore potential correlation between the geometry of the random interface and the
(๐ฅ, ๐ง) variation of locally-averaged D and I magnitudes inside the Brinkman layer. The same
local volume-averaging method described in Section 2.3.3.1 is used. The results for the regular
interface are not shown, as the variations are barely visible using the same range of contour levels.
Figures 2.22(a) and (b) compare the (๐ฅ, ๐ง) distribution of locally averaged D and |I|, respectively.
The spatial average is similar for both contours, a result of relatively weak pressure gradient
discussed earlier. The spatial variation is much more intense for |I| than for D. This indicates
that, though roughness leads to local augmentation of ๐ขโฒ๐ฃ โฒ (predominantly through wake-induced
production), such spatial variation is rather weak and does not significantly affect the (๐ฅ, ๐ง) pattern
of wall-normal โจ๐ขโฉ๐ transport. In contrast, the pressure drag is generated by the small-scale wake of
each grain. Therefore, |I| is more sensitive to the local roughness height. Visual comparison with
the roughness height distribution โโฒ (๐ฅ, ๐ง) in Figure 2.22(c) shows that the distributions of |I| and โโฒ
are somewhat correlated. Indeed, the scatter plot (Figure 2.22(d)) of โโฒ and the spatial fluctuation
of |I| (departure from the ๐ฅ-๐ง mean) displays a positive correlation (with preferred distribution in
Quadrants I and III). This is consistent with the expectation that a high roughness protuberance
tends to generate a large (in ๐ฆ) recirculation region and, consequently, a higher-than-average value
of pressure drag at the corresponding (๐ฅ, ๐ง) location.
As ๐พ โ describes the global effect of both I and D, the ๐พ โ contour in Figure 2.20 also displays
correlation with โโฒ (๐ฅ, ๐ง) distribution; the correlation is negative as a higher local drag serves as a
38
stronger momentum sink.
2.3.3.3 Mechanisms producing vertical mass flux
The hyporheic flux (๐ ๐ป ) quantifies the magnitude of mass exchange across the interface during a
period of time; it is defined as the wall-normal volumetric flow rate at ๐ฆ = โ๐ directed out of the
permeable bed (Elliott and Brooks, 1997b),
โซ
1 ๐
๐ ๐ป (๐) = ๐ฃ(๐ก, ๐ฅ)| ๐ฃ>0 ๐๐ด, (2.17)
๐ด (๐ฅ,๐ง)
๐
where () denotes running time averaging within a window size ๐ and A is the fluid area of the
๐ฅ-๐ง plane at ๐ฆ = โ๐. Figure 2.23(a) shows the variation of the vertical flux, ๐ ๐ป , as a function of
๐. The error bars quantify the variance of the values obtained from the DNS data for each ๐; the
variance is zero for ๐ equal to the total simulation time and increases with decreasing ๐. For ๐
much larger than the large-eddy turn-over time, ๐ฟ/๐ข ๐ , the form-induced wall-normal velocity (e ๐ฃ ) is
the only contributor of ๐ ๐ป , while for ๐ < ๐ฟ/๐ข ๐ the turbulence fluctuation (๐ฃ โฒ) also contributes to
the exchange. Evidently, the contribution of e๐ฃ is more significant than that of ๐ฃ โฒ for both interfaces.
0.15 10
8
0.1 6
4
0.05 2
0
0 -2
0 1 2 3 -0.2 0 0.2
Figure 2.23: (a) Wall-normal flux at ๐ฆ = โ๐ normalized by ๐ข ๐ ๐ฟ2 with different averaging window
sizes ๐๐ข ๐ /๐ฟ for regular ( ) and random ( ) interfaces. (b) Pressure-work term of โจe ๐ฃ 2 โฉ๐
3
budget normalized by ๐ข ๐ /๐ฟ for regular( , โ) and random( ) interfaces.
39
Also, ๐ ๐ป is around 5 times larger for the random interface compared to the other, consistent with
the more intense Reynolds and form-induced stresses near ๐ฆ = โ๐.
The mechanisms governing ๐ฃ โฒ energy transport have been discussed in Section 2.3.2.4 using
the stress budget. For the transport of e ๐ฃ energy, studies on impermeable flows over rough surfaces
(Raupach and Shaw, 1982; Yuan and Piomelli, 2014c) showed that the process is governed by
the work done by the form-induced pressure, ฮ e = โ2โจe ๐ฃ ๐ ๐/๐๐ฆโฉ
e ๐ , as the main source, the wake
production predominantly as a sink. The wall-normal variation of ฮ e is compared in Figure 2.23(b);
the much higher overall value for the random interface is due to the more intense ๐ e fluctuations
(shown in Figure 2.24 in a (๐ฅ โ ๐ฆ) slice for each case) leading to significant ๐ ๐/๐๐ฆ
e near the interface.
Figure 2.24: Time-averaged pressure variation and 3D mean streamlines for (a) regular and (b)
random interfaces. Blue arrows indicate examples of hyporheic flow paths.
40
Thus, the main difference between the ๐ ๐ป generation in the two cases is that the random particle
arrangement leads to larger-scale and more intense time-mean pressure variations which, by working
against the time-mean velocity, generates more intense e ๐ฃ near the interface and promotes vertical
momentum exchange. Meanwhile, ๐ฃ โฒ is also intensified by the pressure work near the interface,
also contributing to ๐ ๐ป , though with a smaller magnitude.
The hyporheic flow paths are also of interests due to their impact on important macroscopic
exchange parameters such as the residence-time distribution and the penetration depth. Below the
region of turbulence penetration (1๐ท below the crest), the flow is steady and, consequently, the
mean streamlines shown in Figure 2.24 can be interpreted as the hyporheic flow paths of particles
released from the upstream.
For the regular interface, the spatial variation of ๐
e shows distributions with alternating negative
and positive values with a length scale of one sphere diameter, representing organized mean
recirculation regions that are shown in Figure 2.14 (๐) and (๐). The main mechanism of vertical
time-mean exchange is due to these organized recirculation regions at the interface. The relatively
low ๐พ โ at the uppermost layer of the regular interface limits communication between the flows
inside the sediment and in the region above; thus, in the sediment, there is little vertical exchange
and the lateral flow paths are long.
In contrast, for the random interface, the hyporheic flow paths display a multiscale pattern with
short flow paths near the interface and longer flow paths in the shape of circular arcs reaching deeper
into the bed. These paths are similar to those visualized through dye transport for a flat bed in the
experiments of Packman et al. (2004). Near the bed, the strong ๐ e gradients of larger coherence
length induce flow infiltration near the bed, reminiscent of bedform-induced advective pumping
(Packman et al., 2004). This also results in more intense e ๐ฃ , increasing form-induced mixing. The
streamwise lengths of the flow paths appear to be reduced as the wall-normal communication is
strengthened.
41
2.4 Conclusions and discussions
We reported DNS results for fully-developed turbulent flows in an open channel on grain-resolved
sediment beds with either regular or random grain positioning at the sediment-water interface. The
random interface results in a particle roughness with significantly larger root-mean-square height
and longer horizontal correlation lengths. These length scales are much smaller compared to those
typical of a bedform. Thus, the interface height variation is considered as roughness. A comparison
with experimental measurements using randomly poured spheres shows that the case with random
interface yields matching flow statistics. It is thus considered as a good approximation of a realistic
sediment.
Detailed DNS results obtained with regular and random interfaces with matching ๐
๐ ๐ = 395
and ๐
๐ ๐พ = 2.6 are compared. The differences in the effects on the turbulence and the mass and
momentum exchange across the interface are summarized as follows. (1) As a main difference
between the two cases, the larger roughness length scales brought by the random interface results
in less organized distribution of mean recirculation regions at the interface and more intense ๐ e
variations at larger scales. The result is a higher e ๐ฃ production and, consequently, more isotropic
form-induced stress tensor and significant form-induced shear stress (that is comparable to the
Reynolds shear stress). It is shown that e ๐ฃ , as opposed to ๐ฃ โฒ, is the main contributor to the
wall-normal hyporheic flux across the interface for both regular and random interface roughness.
Consequently, the hyporheic flux is significantly higher for the random interface than for the
regular one. (2) The vertical spread-out of the mean-shear-layer distribution near the random
interface leads to a more even TKE redistribution for ๐ฃ โฒ and ๐ค โฒ fluctuations and consequently
a more isotropic Reynolds stress tensor. The form-induced fluctuations also modulate the total
production of Reynolds stress through conversion between TKE and WKE, as well as introducing
an additional diffusion process that transports energy from the low-TKE region to the high-TKE
region. (3) More intense mixing is observed near the random interface due to augmented Reynolds
and form-induced shear stresses; this is reflected by a deeper turbulence penetration (44% higher
๐ฟโ๐ ) and a higher apparent permeability (twice higher ๐พ โ ) in the Brinkman layer. The local apparent
42
permeability in this layer is negatively correlated with the local height of roughness. Moreover, the
strengthened communication between the surface and subsurface flows for the random interface
leads to shortened and deeper-reaching hyporheic flow paths.
Regarding how the type of interface roughness influences the effects of wall permeability, the
results show that the penetration depths of mean shear and turbulence are significantly increased by
the random roughness. However, the roughness type does not appear to noticeably modify the von
Kรกrmรกn constant (which is lower on a permeable wall than on an impermeable one), as long as the
same type of packing is used inside the sediment. Lastly, the Kelvin-Helmholtz (K-H) instability
due to wall permeability and the K-H rollers, if present, may also be affected differently by the
two roughness types. The K-H instability is, however, expected to be very weak in this work as
๐
๐ ๐พ โ 2.5. Manes et al. (2011) systematically evaluated the effects of ๐พ with ๐
๐ ๐พ varying from 0
to 17. Within this range, the K-H instability and signature of the rollers were observed for ๐
๐ ๐พ โ 17
and not for ๐
๐ ๐พ โ 8, for example. The effects of interface roughness on K-H instability and the
resulting coherent motions for flows with high ๐
๐ ๐พ are interesting questions for future work.
It should be noted that the specific configuration of the regular packing may also affect the
drag, turbulent statistics and structure. For example, the hexagonal particle packing is expected
to give more drag and higher resistance to the streamwise velocity fluctuations. We have run an
additional DNS simulation with hexagonal arrangement at the interface (not shown herein) and
compared various flow quantities with the cubic arrangement. It was observed that, though the
hexagonal arrangement gave slightly higher ๐ถ ๐ and anisotropy of the form-induced stresses, such
differences were much smaller compared to the difference between the cubic and and the random
arrangements. The reason is probably due to the similar roughness length scales among various
regular arrangements and the organized mean-velocity pattern (e.g. mean recirculation regions)
they generate. For this reason, the cubic arrangement discussed in details herein is considered as a
representative example of regular type of packing.
The results demonstrate that subtle details of the particle roughness alone, in the absence of
bedform, can affect significantly the dynamics of the turbulence and the form-induced fluctuations.
43
Such effects translate to large differences in apparent permeability, hyporheic flux, and subsurface
flow paths of passive scalars. Future work is needed to fully understand the effects of particle
roughness on scalar and momentum transport with systematically varying roughness lengths (๐โ+ ,
๐+ ), roughness texture, as well as the effect of varying bulk permeability with random interface
roughness.
44
CHAPTER 3
QUANTIFYING THE EFFECTS OF BED ROUGHNESS ON TRANSIT TIME
DISTRIBUTIONS VIA DIRECT NUMERICAL SIMULATIONS OF TURBULENT
HYPORHEIC EXCHANGE
3.1 Introduction
The pore-resolved direct numerical simulations (DNS) of interface turbulence carried out in Chap-
ter 2 revealed that, even in the absence of bedforms, the grain-scale roughness of a macroscopically
flat bed leads to significant fluid volumetric flux into the sediment. Also, statistics and structure of
the turbulent flow were shown to depend on the texture of roughness. Similar observations were
made for a macroscopically flat bed by Chandler et al. (2016), Grant et al. (2020a) and Kim and
Kang (2020).
In this chapter, the work in Chapter 2 is expanded by documenting the effects of grain-scale
roughness and its texture on macroscopic exchange quantities including the fluid volumetric flux
into the sediment, subsurface flow paths, and transit time characteristics. The primary questions
that will be addressed in this chapter include:
1. How does bed roughness on a macroscopically flat bed influence fluid parcel transit times?
2. How do TTDs change with tracer particle entry/exit locations within a three-dimensional
flow field and with the inclusion of molecular diffusion?
To this end, turbulent open-channel flows bounded by synthesized sediment beds similar to those
studied in Chapter 2 were simulated. Computed flows for two sediment beds that share the
same bed permeabilities but have different bed roughness arrangements are compared. A forward
particle tracking method is used to calculate the transit times, considering the mechanism of
time-mean velocity advection with and without molecular diffusion. The paper is organized as
follows. Section 3.2 describes the details of the DNS simulations, the sediment bed synthesis, and
the method used to calculate transit times. Section 3.3 compares the transit time results for the
45
two cases examined and analyzes the effects of molecular diffusion, followed by conclusions in
Section 3.4.
3.2 Methodology
3.2.1 Case configuration and parameters
Figure 3.1: Simulation domain and synthesized sediment beds colored by ๐ฆ in (a) the regular and
(b) the random cases. (c) One-dimensional power spectral densities ๐ธ ๐ of the interface height
fluctuations for Case 1 (regular case, ) and Case 2 (random case, , red). (red) Slope
of power-law decay for Case 2; ๐ ๐๐๐ is the root-mean-square (rms) of interface height fluctuations.
Fully developed open-channel flows are simulated using the same DNS method mentioned
in previous Chapter. The simulation domain includes both the surface and sub-surface flows
(Figure 3.1(a,b)). ๐ฟ is the thickness of the turbulent boundary layer (i.e., open-channel height)
measured from the sediment crest; ๐ป๐ is the sediment depth measured from the crest to the bottom
boundary of the simulation domain. Symmetric boundary conditions are applied at both top and
bottom boundaries of the domain. Periodic conditions are applied at ๐ฅ and ๐ง boundaries. A
constant mean pressure gradient is used to drive the flow. The elevation of ๐ฆ = 0 is set at the
46
Interface ๐
๐ ๐พ ๐
๐ ๐ ๐ท/๐ฟ ๐ท+ ๐ป๐ /๐ฟ (๐ฟ ๐ฅ , ๐ฟ ๐ง )/๐ฟ (ฮ๐ฅ + , ฮ๐ฆ +min , ฮ๐ง+ )
Case 1 Regular 2.62 395 0.20 79 1.6 (6, 3) (1.5, 0.19, 1.5)
Case 2 Random 2.62 395 0.20 79 1.6 (6, 3) (1.5, 0.19, 1.5)
Table 3.1: Summary of simulation parameters. ๐ท is the grain diameter; ๐ฟ ๐ฅ๐ is the simulation
domain size in ๐ฅ๐ . ฮ๐ฅ + , ฮ๐ฆ + , and ฮ๐ง+ are DNS grid sizes in ๐ฅ, ๐ฆ and ๐ง, respectively, normalized
using the viscous length scale ๐/๐ข ๐ .
virtual originโan offset of the ๐ฆ axis such that the mean velocity profile in the overlap region of
the turbulent flow recovers its logarithmic form. The elevation of the virtual origin was determined
based on a fitting procedure (Breugem et al., 2006; Manes et al., 2011); it falls between 0.3๐ท and
0.5๐ท below the sediment crest, where ๐ท is the grain diameter.
The simulation parameters are listed in Table 3.1. Two cases are considered, one with a regular
interface (Case 1) in which the arrangement of grains in the uppermost layer is regular and the
other with random arrangement of grains (termed โrandom interfaceโ, Case 2). In both cases,
the same porosity of 0.41 is imposed in the sediment sufficiently far from the interface. The
permeability Reynolds number and the friction Reynolds number are the same for both cases:
โ
๐
๐ ๐พ = ๐พ๐ข ๐ /๐ = 2.62 (where ๐พ is the permeability) and ๐
๐ ๐ = ๐ฟ๐ข ๐ /๐ = 395. Following
Voermans et al. (2017), ๐ข ๐ is obtained from the maximum magnitude of the bed-normal profile of
the total shear stress, which is the sum of the viscous shear stress, the Reynolds shear stress and
the form-induced shear stress. The total numbers of grid points are 1536, 779 and 768 in ๐ฅ, ๐ฆ and
๐ง directions, respectively. The present DNS simulations are similar to the simulations reported by
Shen et al. (2020). The only difference is that the sediment thickness herein is doubled (๐ป๐ /๐ฟ = 1.6
versus 0.8 in Shen et al. (2020)), to allow for deeper subsurface flow paths, which are important for
calculation of large transit times but less important for interfacial turbulence. The total simulation
time used for data collection is around 10 large-eddy turn-over times (LETOTs), defined as ๐ฟ/๐ข ๐ ; a
LETOT represents the time scale of the evolution of the largest turbulent eddy in the open channel.
Sufficiency of the data sampling in time was validated by comparing both surface and subsurface
flow statistics, as well as TTDs, obtained with 10 LETOTs and those obtained with 20 LETOTs.
1/2
The differences were small (e.g., up to 1% in subsurface dispersive velocities, e ๐ข2 , as shown in
47
the Appendix B Figure B.4).
For both cases, the depths of the channel and the bed (both measured from the virtual origin)
are 5 times and 8 times of the grain diameter, respectively. A boundary layer thickness of 5 times
the grain size is much smaller than in many streams in reality and, consequently, the turbulence
near the interface may not be identical to those in many real streams. However, the present setup
is designed to represent streams for the regime with ๐
๐ ๐พ of ๐(1) where the form-induced velocity
is the main component of the subsurface fluid motion; the interfacial turbulence weakly affects
the majority of the parcel transit times in comparison to the form-induced velocities. In addition,
๐ท/๐ฟ of the same order of magnitude were used in existing numerical and experimental studies that
reported grain-scale physics (Manes et al., 2009; Dey and Das, 2012; Fang et al., 2018; Han et al.,
2018; Voermans et al., 2017, 2018).
The porous sediment beds are modeled as closely packed monodisperse hard spheres. In
the bulk of the sediment, the sphere positioning is determined based on molecular dynamics
simulations (Plimpton, 1995) and is random for both cases. For the uppermost-layer spheres,
two types of distributions are imposed, representing regular (Case 1) and random (Case 2) bed
roughness geometries, respectively. The random case corresponds to significantly larger root-
mean-square roughness height and longer horizontal correlation lengths. As these length scales
are much smaller compared to ๐ฟ and the scales typical of a bedform, the bed height variation
is considered as roughness. Figures 3.1(a,b) compare the sphere distributions in Cases 1 and 2,
showing the difference in grain arrangement at the top of the sediments. The sediment in Case
2 was constructed such that it yields turbulence statistics that match those observed for randomly
poured spheres in the experimental studies of Voermans et al. (2017) (case L12 therein). Case 2 is
thus considered as a good approximation of a realistic sediment.
To highlight the difference between the two roughness geometries, the one-dimensional (1D)
power spectra of height fluctuations of uppermost-layer spheres (๐ธ ๐ ) are compared in Figure 3.1(c)
as functions of the streamwise wavenumber (๐
1 = 1/๐ 1 , where ๐ 1 is the streamwise wavelength).
Both cases show a white noise regime at small wavenumbers, characterized by almost constant ๐ธ ๐ .
48
Aubeneau et al. (2015) calculated the power spectra of elevations of sand beds produced in lab
flume experiments under various discharges. They also observed a white noise region in the power
spectrum for small wavenumbers. For the random case the white noise regime applies to ๐ 1 /๐ท larger
than around 1.5. For smaller ๐1 , ๐ธ ๐ shows a power-law decay with a slope of around โ2.4. From
this power-law decay, a fractal dimension ๐ท ๐ can be identified as ๐ท ๐ = ๐ธ + (3 โ ๐ฝ)/2 = 2.3, where
๐ฝ = 2.4 is the fitted slope of the decay and ๐ธ = 2 (for a 3D surface) is the Euclidean topological
dimension. Aubeneau et al. (2015) also observed power-law decays at large wavenumbers, though
with different slopes of โ1.0 to โ1.5. In comparison, for the regular case the white noise region
covers wavenumbers up to ๐ท. At larger wavenumbers ๐ธ ๐ is characterized by higher-frequency
harmonics, which is a mathematical consequence of the slight clipping (or flattening) at the top
of the sphere surfaces that results from the spatial discretization of the sphere geometry. These
observations show that the bed roughness of Case 2 is close to a self-affine fractal (similar to real
sediment beds) while Case 1 is characterized by organized grain-wise protuberances.
3.2.2 Transit time calculation based on particle tracking
For an unsteady flow in general, the mass conservation equation for the storage of water of a certain
age in the sediment is (Botter et al., 2011; Harman, 2015)
๐๐(๐ก) ๐ ๐ (๐, ๐ก) ๐๐(๐ก) ๐ ๐ (๐, ๐ก)
= โ๐(๐ก) โ
๐โ๐ (๐, ๐ก) โ , (3.1)
๐๐ก ๐๐
where ๐(๐ก) is the instantaneous total water storage in the sediment at time ๐ก, ๐(๐ก) is the instantaneous
flux at the outflow, ๐ ๐ (๐, ๐ก) is the instantaneous RTD, defined as the probability density function
(pdf) of the age (๐) of all the water parcels in the sediment at time ๐ก, and โ ๐โ๐ (๐, ๐ก) is the backward
TTD, defined as the probability that the age of a water parcel at the outflow at time ๐ก is equal to
๐ (Harman, 2015). The only assumption used in Equation (3.1) is that there is one inlet and one
outlet of the system. In this study, the inlet is the collection of regions on the water-sediment
interface where the bed-normal velocity (๐ฃ) is negative, while the collection of regions where ๐ฃ is
positive serves as the outlet. The water age is considered zero at the inlet.
49
Note that different definitions of RTD were used in the literature. For example, the cumulative
distribution function of RTD was defined as the fraction of solute entering the sediment bed in
a short time near ๐ก = 0 and exiting the bed by time ๐ (Grant et al., 2020b). The residence time
function of EB97 is trivially different from this definition. In this article, however, the definitions
of Harman (2015) are followed.
We consider fluid parcels advected by the three-dimensional, time-averaged velocity, ๐ข๐ . The
advection by instantaneous turbulent fluctuations ๐ขโฒ๐ is not accounted for. This is because, at the
present ๐
๐ ๐พ = 2.6, the magnitudes of normal Reynolds stress components are significantly smaller
than those of the form-induced stresses in the sediment below the Brinkman layer, as shown by
the DNS data of Shen et al. (2020) for the present cases. At a higher ๐
๐ ๐พ range, however, the
effect of ๐ขโฒ๐ on the transit times is likely to be more important inside the sediment (Voermans et al.,
2017). Kim et al. (2020) showed that at ๐
๐ ๐พ = 50 and ๐
๐ ๐ โ 2, 000 significant turbulent motions
in the bedโincluding downwelling and upwelling flowsโoccur and are modulated by large-scale
turbulent motions in the surface flow. The present analyses, therefore, are relevant for the regime
of ๐
๐ ๐พ = ๐(1) only, where the interfacial turbulence is present but weak compared to dispersive
velocities.
The consideration of the water parcels advected by the time mean velocity implies a steady state
assumption for Equation (3.1), which becomes
๐ ๐ ๐ (๐)
๐โ ๐โ๐ (๐) + ๐ = 0. (3.2)
๐๐
Equation (3.2) shows that the backward TTD and the RTD are related; one can be calculated from
the other given known values of the outflow flux and water storage. In this work, we calculate โ ๐โ๐ (๐)
based on its definition. The discussion is focused on โ ๐โ๐ (๐) and the corresponding cumulative
โ
โ
distribution function, ๐ ๐ (๐). The RTD can be calculated based on Equation 3.2; details are
presented in the Appendix B.
In the present particle tracking approach, a fluid parcel is considered as a tracked particle.
To calculate the transit time (or the maximum water age achieved by an individual parcel in the
sediment), fluid parcels are traced starting from a ๐ฆ elevation near the interface throughout the
50
(๐ฅ, ๐ง) plane. Only the fluid parcels entering the sediment (i.e., those released at locations where
๐ฃ < 0) are tracked. The tracked paths are in general three-dimensional. The parcel position ๐ฅ๐ (๐ก) at
a time ๐ก (with ๐ก = 0 at the parcel release time) is tracked based on
๐๐ฅ๐
= ๐ข๐ [๐ฅ๐ (๐ก)], (3.3)
๐๐ก
where ๐ก is a fictitious time used for particle tracking, not to be confused with the simulation time in
DNS. Equation (3.3) is discretized and solved using the explicit Euler scheme. Once a subsurface
travel path is determined, the transit time for the corresponding fluid parcel is calculated as the total
time it spends along the subsurface path.
All tracked fluid parcels share the same release elevation (๐ฆ 1 ) and the same โexitโ elevation (๐ฆ 2 ).
If ๐ฆ 1 โค ๐ฆ 2 , a parcel is considered to have left the sediment once it passes ๐ฆ 2 . If ๐ฆ 1 > ๐ฆ 2 , the tracked
parcel first goes downward, passing the ๐ฆ 2 elevation, and then upward toward the ๐ฆ 2 elevation; in
this case, the exit is considered to occur when the parcel passes ๐ฆ 2 the second time. The periodic
boundary conditions in ๐ฅ and ๐ง are used: parcels leaving the domain at boundary ๐ฅ/๐ฟ = 6 reenters
the domain at boundary ๐ฅ/๐ฟ = 0, and vice versa, for up to twice of the domain lengths. Longer
tracking lengths up to 4 times of the domain lengths in ๐ฅ and ๐ง were also tested, yielding similar
transit time distributions. One particle is released at each computational grid point in the fluid
domain at ๐ฆ = ๐ฆ 1 . In total, around ๐(106 ) particles are released. To test the result convergence,
we used twice the number of released particles, which led to differences of less than 1% in the
โ
โ
fitted TTD slope and the median transit time value. The cumulative backward TTD, ๐ ๐ (๐), is
calculated as the fraction of the number of tracked subsurface paths that correspond to transit times
โโ
shorter than ๐. โ ๐โ๐ (๐) is then obtained as ๐ ๐ ๐ /๐๐. The calculated โ
๐โ๐ (๐) is smoothed based on a
moving averaging procedure. The size of the averaging window widens with the transit time value
as approximately 0.1๐.
To validate the implementation of the particle tracking method, the computed residence time
function as defined in EB97 is compared to the analytical solution for a subsurface flow induced
by a 2D sinusoidal interface pressure distribution. Details of the validation are included in the
Appendix B. The comparison shows that the particle-tracked residence time function matches
51
well with the analytical solution given a sufficiently fine spatial resolution and a sufficiently deep
sediment domain (Figure B.2).
Figure 3.2: (a) Sketch of typical tracked flow paths for the random interface (Case 2). The
cumulative backward transit time distribution (b) and its pdf (c) are shown for two sets of entry
and exit locations: ๐ฆ 1 = ๐ฆ 2 = 0 ( ) and ๐ฆ 1 = ๐ฆ 2 = โ3๐ท ( ) for the regular interface (black
lines) and the random interface (red lines). (d) Modified pdf showing the probability density in
evenly spaced logarithmic increments. In (b), the horizontal dashed lines ( , blue) mark the 50th
(median) and 95th percentiles. In (c), the dashed blue lines ( , blue) show the fitted power law
with slopes indicated.
In the analyses of EB97, detailed flow information around sediment grains is not available. The
water column and the subsurface flow were delineated by a zero-thickness interface separating the
surface flow and the subsurface flow where the Darcy model applies. In that context, a particle
52
tracking method for residence time calculation assumes that the release elevation (๐ฆ 1 ) and exit
elevation (๐ฆ 2 ) of the tracked fluid parcels are both at the exact water-sediment interface. In the
present DNS, there is a transition layer between the surface flow and the Darcy-flow region. This
layer includes the roughness sublayer (Raupach et al., 1991)โa near-bed layer where form-induced
stresses are importantโand the Brinkman layerโthe region of mean shear layer penetration,
๐โจ๐ขโฉ/๐๐ฆ > 0. Varying entry and exit elevations parcels are tested. Examples of tracked flow paths
are sketched in Figure 3.2(a).
3.3 Results
3.3.1 Effect of bed roughness geometry on transit times
โโ
Figures 3.2(b) and 3.2(c) compare ๐ ๐ (๐) and โ ๐โ๐ (๐) with the same entry and exit elevations for
two sets of values: ๐ฆ 1 = ๐ฆ 2 = 0 (near the sediment crest) and ๐ฆ 1 = ๐ฆ 2 = โ3๐ท (a sufficiently deep
elevation beyond which the median ๐ of the โ ๐โ๐ (๐) distribution becomes largely independent of
the entry and exit elevations for both cases, as shown later in Figure 3.5). The elevation โ3๐ท is
equivalent to โ0.6๐ฟ. The transit times are normalized by the LETOT (๐ฟ/๐ข ๐ ) of the open-channel
turbulence, representative of the time scale of largest turbulent coherent motions in the water
column. For large transit times stronger fluctuations of the TTD curves are seen in Figure 3.2(c),
since the sampling of these long paths available from the DNS domain is extremely limited. The
โ๐โ๐ (๐) distributions display power-law decays with a slope of around โ1.3 at large transit times
for both cases. The existence of a power-law tail indicates that the sediment domain in the DNS
simulations is sufficiently large (in both ๐ฅ and ๐ฆ) to capture deeper flow paths, and is consistent
with experimental observations by Aubeneau et al. (2015) and results of numerical simulations
obtained by Lee et al. (2020) of exchanges induced by multiscale bed geometries. Note that,
according to Equation (3.2), a power-law decay of the TTD (โ ๐โ๐ ) is associated with a power-law
decay of the RTD (๐ ๐ ), which is shown in Figure B.1 in the Appendix B. For both cases, a deeper
entry elevation excludes shorter flow paths close to the surface and, consequently, increases the
probability of longer paths. This is reflected by the shifting of the TTDs toward larger times in
53
(๐ฆ 1 , ๐ฆ 2 ) Characteristic Case 1 (A) Case 2 (A) Case 1 (MD) Case 2 (MD)
(0, 0) ๐50 0.294 0.368 0.467 0.574
(0, 0) ๐95 10.6 10.6 110 80.0
(โ3๐ท, โ3๐ท) ๐50 7.33 8.67 13.6 13.8
(โ3๐ท, โ3๐ท) ๐95 205 393 8.35 ร 105 2.17 ร 106
Table 3.2: Various characteristics of transit times calculated based on mean advection only and
both mean advection and molecular diffusion, at two elevations of particle entry (๐ฆ 1 ) and exit (๐ฆ 2 ).
Transit times are normalized by ๐ฟ/๐ข ๐ . Note. โAโ denotes calculation based on time-mean
advection alone; โMDโ denotes calculation accounting for both time-mean advection and
molecular diffusion (discussed in Section 3.3.3).
Figure 3.2(c) as ๐ฆ 1 and ๐ฆ 2 change from 0 to โ3๐ท. The change of entry and exit elevations does not
significantly modify the slope of the power-law tail.
When the particle release location at 0 was used, the โ ๐โ๐ distribution in the case with regular
roughness displays a relatively sharp drop of probability for dimensionless ๐ at around 400. This
is thought to be related to the organized fluid motions around the regularly arranged grains at the
surface, leading to preferential values of ๐ associated with these short paths. A similar drop in
probability is also evident at dimensionless ๐ of around 200 to 400 in both cases, when a particle
release elevation of โ3๐ท was used. This is because the algorithm tracks undulating segments of the
large number of paths that reach not far beyond ๐ฆ = โ3๐ท (visualized in Figures 3.3(e,f)). Lee et al.
(2020) studied bed geometries synthesized with systematically varying spectral characteristics.
They also observed variations in the slope of the RTD and suggested that they are connected to the
flow recirculation regions induced by bed morphology.
An alternative approach to plot the transit time distribution, different from โ ๐โ๐ (๐), is to divide
the unit area under the curve into evenly spaced logarithmic increments of ๐. This allows the actual
probability density to be shown in a plot with a logarithmic ๐ axis. This distribution, shown in
โโ
Figure 3.2(d), is defined as ๐ ๐ ๐ /๐ log10 ๐ and is equivalent to 2.303๐ โ ๐โ๐ (Grant et al., 2020b).
To quantitatively compare the TTDs, the 50th and 95th percentiles (denoted as ๐50 and ๐95 ,
respectively) of the โ๐โ๐ distributions are calculated. Specifically, ๐50 represents the central tendency
of the distribution, while ๐95 is considered as a characteristic ๐ associated with the rare, long-transit-
54
time subsurface paths. Their values are tabulated in Table 3.2 under columns โCase 1 (A)โ and
โCase 2 (A)โ. The comparison shows that the two cases give rather similar ๐50 values (up to 10-20%
different) for both elevations. For ๐ฆ 1 = ๐ฆ 2 = โ3๐ท, ๐95 is significantly higher in Case 2 than in Case
1 (393 versus 205), due to longer and deeper reaching paths, which will be shown in Figure 3.3(g,h).
An important conclusion that follows from this result is that the difference in bed roughness texture
extends deep into the sediment and significantly modifies the transit time values.
Figure 3.3: Time-mean subsurface flow paths corresponding to the median (50th-percentile)
transit time (a,b,e,f) and the 95th-percentile transit time (c,d,g,h) for entry/exit elevations at
๐ฆ 1 = ๐ฆ 2 = 0 (a,b,c,d) and ๐ฆ 1 = ๐ฆ 2 = โ3๐ท (e,f,g,h), for the regular (a,c,e,g) and random (b,d,f,h)
interfaces.
55
To explain the differences in the TTD characteristics, Figure 3.3 shows tracked subsurface flow
paths for ๐ฆ 1 = ๐ฆ 2 = 0 (a-d) and for ๐ฆ 1 = ๐ฆ 2 = โ3๐ท (e-h) for both cases. Only the paths with
travel times equal to ๐50 and ๐95 are plotted. In each subplot, around 200 pathlines are visualized.
For both cases the paths associated with ๐50 are short and do not reach beyond one grain diameter
below ๐ฆ 1 . Those associated with ๐95 , however, reach far deeper and are described as follows. For
the regular interface, relatively uniform small-scale arc-shaped paths are observed (due to flows
around regularly packed grains at the bed surface). For the random case the tracked paths are more
complex, with short paths near the interface and also longer, deeper paths. As the transit time of a
water parcel depends on the depth and the length of the flow path, ๐95 in Case 2 is longer than that
in Case 1. For paths at ๐ฆ 1 = ๐ฆ 2 = โ3๐ท, Case 1 also yields multiscale flow paths similar to those
for Case 2, as the short paths near the interface for Case 1 are now excluded by the deeper entry
elevation. The paths in Case 2 reach deeper into the sediment than for Case 1, leading to longer
transit times. These observations indicate that the flow pattern throughout the bed is dependent on
the roughness characteristics.
In Figure 3.3, nearly all of the flow paths are oriented downstream in ๐ฅ. Exceptions are the
flow paths located in the lower portion of the recirculation regions near the interface, where the
flow is directed upstream. The very high fraction of subsurface flow paths oriented downstream
is a major difference from flow paths induced by bedforms, such as those induced by triangular
bedforms in the experimental study of Elliott and Brooks (1997a), where a significant fraction
of large-scale upstream-oriented subsurface flow paths were found. This difference is probably
due to the negligible hydrostatic pressure variations along the interface in the present cases as the
roughness height is much lower than the open-channel height.
To characterize the velocity along the subsurface flow paths, Figure 3.4 compares the wall-
normal velocity magnitudes (|๐ฃ|) that are conditionally averaged along the down- and up-going
portions of the flow paths, respectively. The result is a function of ๐ฆ and the bed-normal flow path
direction. Flow paths tracked with ๐ฆ 1 = ๐ฆ 2 = 0 are used. At a given ๐ฆ near the interface, |๐ฃ| for the
random case is much larger than the value for the regular case. This is consistent with the deeper
56
flow paths for the former (Figure 3.3(d)) than for the latter (Figure 3.3(c)). In addition, for both
cases the speed of downward flow is higher compared to that of upward flow. This difference has
an impact on the transit time values when the entry elevation differs from the exit one.
y/
|v|/uโง
Figure 3.4: Wall-normal variation of ๐ฃ magnitude conditionally averaged along downward (red)
or upward (black) flow paths tracked at ๐ฆ 1 = ๐ฆ 2 = 0, for regular ( ) and random ( ) cases.
Next, independent variations of entry (๐ฆ 1 ) and exit (๐ฆ 2 ) elevations are studied. Both ๐ฆ 1 and ๐ฆ 2 are
varied between the sediment crest and ๐ฆ/๐ฟ = โ0.9. An ensemble of particle-tracking simulations
based on 42 values each of ๐ฆ 1 and ๐ฆ 2 (total 1764) was carried out. For each (๐ฆ 1 , ๐ฆ 2 ) pair, ๐50 is
calculated and plotted in a map shown in Figure 3.5. The variation of ๐95 was also evaluated and a
similar comparison was observed; the ๐95 plots are shown in the Appendix B Figure B.5.
First, focus on the diagonal line represented by ๐ฆ 1 = ๐ฆ 2 . The short transit times of ๐(1)๐ฟ/๐ข ๐ at
๐ฆ 1 = ๐ฆ 2 โ 0 are mostly due to mild undulations of subsurface paths (Figure 3.3a,b). For sufficiently
deep entry/exit elevations, ๐50 plateaus to values of ๐(10)๐ฟ/๐ข ๐ . The dependence of transit times on
the entry and exit elevations of particles is not restricted to elevation variation within the Brinkman
layer, which corresponds approximately to ๐ฆ/๐ฟ โ (โ0.1, 0.1). Beyond this layer, deeper entry or
exit elevations will still result in progressively longer transit times due to progressive exclusion of
57
shorter flow paths, until only long flow paths remain.
Next, small values of |๐ฆ 1 โ ๐ฆ 2 | produce short transit times due to a large fraction of tracked
paths being associated with short flow paths that do not extend far beyond the elevation of ๐ฆ 1 , while
with large |๐ฆ 1 โ ๐ฆ 2 | differences a larger fraction of the tracked paths are longer and deep-reaching
ones, which yields much longer transit times. The regions above the diagonal line (i.e., ๐ฆ 1 > ๐ฆ 2 )
and below the line (i.e., ๐ฆ 1 < ๐ฆ 2 ) are not symmetric. The median transit times with ๐ฆ 1 below ๐ฆ 2 are
larger than those with ๐ฆ 1 above ๐ฆ 2 at the same |๐ฆ 1 โ ๐ฆ 2 | difference. This is because the downward
motion of a fluid parcel in the bed is, on average, faster than the upward motion, as shown in
Figure 3.4. One also observes that the median transit times with ๐ฆ 1 > ๐ฆ 2 are less well-converged
statistically than those with ๐ฆ 1 < ๐ฆ 2 , as shown by the more fluctuating contour lines. This is due to
the fewer flow paths used for averaging when ๐ฆ 1 > ๐ฆ 2 in the current computation procedure (see
Appendix B Figure B.3). At large ๐ฆ 1 -๐ฆ 2 separations, much longer transit times are observed for the
regular case compared to the random one. This difference is probably a result of the significantly
smaller magnitude of bed-normal velocity along flow paths in the regular case for ๐ฆ/๐ฟ > โ0.5
(a) (b)
Figure 3.5: Median transit time normalized by ๐ฟ/๐ข ๐ with independent variations of ๐ฆ 1 (entry) and
๐ฆ 2 (exit) elevations, for regular (a) and random (b) cases. Similar variations of ๐95 are shown in
the Appendix B Figure B.5.
58
(Figure 3.4).
These results highlight the fact that subsurface flow paths are multiscale in nature and typically
extend deeper beyond the Brinkman layer. They also show that both the entry and exit elevations
are important parameters in measuring and reporting transit times in the presence of bed roughness;
however, these details are not often presented in the literature. The roughness geometry affects not
only the transit time when ๐ฆ 1 is set equal to ๐ฆ 2 , but also the dependence of transit time on ๐ฆ 1 and
๐ฆ 2 when ๐ฆ 1 โ ๐ฆ 2 .
In the next two sections we address the following questions: (1) What factors determine the (๐ฅ, ๐ง)
locations of the entries and exits of subsurface flow paths (and consequently TTD characteristics)?
(2) How important is the molecular diffusion for solute transit time?
3.3.2 Interfacial ๐ฃ distribution
Figure 3.6 shows the typical flow pattern around a roughness protuberance formed as a grain
cluster. A stagnation point appears upstream of the protuberance, leading to a region with local ๐
higher than the spatial average (indicated by ๐ e > 0). In addition, a recirculation region appears
downstream, leading to a local low-๐ region (indicated by ๐ e < 0). Such interfacial pressure
Figure 3.6: Sketch showing how local high-pressure (stagnation point) and low-pressure
(recirculation) regions lead to entry and exit of fluid parcels (red spheres and red lines),
specifically, in Case 2. Roughness geometry. ๐(๐ฅ, ๐ฆ, ๐ง) and ๐ฃ(๐ฅ, ๐ฆ, ๐ง) are the time-averaged
e ๐ฆ, ๐ง) = ๐(๐ฅ, ๐ฆ, ๐ง) โ โจ๐โฉ(๐ฆ) is the form-induced
values of pressure and wall-normal velocity; ๐(๐ฅ,
pressure, i.e., deviations from the plane-averaged value.
59
variation is in line with the description of local pressure variation on a bedform (triangular dune)
by Savant et al. (1987). These local pressure maxima and minima lead to bed-normal pressure
gradients at the bed surface: ๐๐/๐๐ฆ > 0 upstream and ๐๐/๐๐ฆ < 0 downstream.
Figure 3.7: Bed-parallel contours for the random interface (Case 2) at ๐ฆ = 0: (a) bed-normal mean
velocity, (b) form-induced pressure, (c) bed-normal gradient of mean pressure, and (d) transit time
of flow paths entering at each location. Only 1/4 of the domain is shown. All quantities in (a-d)
are normalized by ๐ฟ and ๐ข ๐ . Joint probability density functions of ๐๐/๐๐ฆ and ๐ฃ in the interface
region, for Case 1 (e) and Case 2 (f).
The bed-normal mean momentum equation near the interface connects the pressure gradients
60
to the entry and exit of subsurface flow paths:
๐ท๐ฃ ๐๐ฃ ๐๐ฃ 1 ๐๐ ๐๐ขโฒ๐ ๐ฃ โฒ ๐2๐ฃ
โก + ๐ข๐ =โ โ +๐ , (3.4)
๐ท๐ก ๐๐ก ๐๐ฅ๐ ๐ ๐๐ฆ ๐๐ฅ๐ ๐๐ฅ๐ ๐๐ฅ๐
where ๐ท๐ฃ/๐ท๐ก is the material derivative of ๐ฃ advected by the time-mean flow, and the three terms
on the right-hand-side represent forces due to, respectively, time-mean pressure, Reynolds stress
and viscous stress. Equation (3.4) describes how various factors (on the right-hand-side) affect the
time-mean ๐ฃ velocity of a fluid parcel. It is hypothesized that the local pressure gradient (the first
term on the right-hand-side of (3.4)) is the main factor in determining the entry (๐ฃ < 0) and exit
(๐ฃ > 0) locations of subsurface flow paths. Specifically, a positive ๐๐/๐๐ฆ upstream of the obstacle
would lead to negative ๐ฃ values (entries of flow paths) for a fluid parcel initially with ๐ฃ = 0, while
a negative ๐๐/๐๐ฆ downstream would lead to exits of flow paths. Such mechanism is sketched in
Figure 3.6. This hypothesis is to be validated next.
Figures 3.7(a-d) compare contours of ๐ฃ, ๐, e and ๐๐/๐๐ฆ in a bed-parallel plane at ๐ฆ = 0 for
the random case. The correlation between positive-๐ e regions, positive ๐ฆ-gradients of pressure,
and negative ๐ฃ values (entries), is apparent. These regions are typically seen upstream of clusters
of grains serving as roughness protuberances. Similarly, downstream of these protuberances,
one typically observes low pressure values, negative pressure gradients, and positive ๐ฃ (exits).
Figure 3.7(d) shows the contour of transit time of flow paths entering the subsurface from this
plane. Regions with the highest pressure values are seen to generate the longest transit times.
The local correlation between ๐ฃ and the bed-normal pressure gradients is better demonstrated
using contours of their joint probability density function (jpdf) in Figures 3.7(e,f) for the two cases,
respectively. Since the local ๐ฆ-gradient of ๐ varies significantly in ๐ฆ between the crest elevation
and the bottom of the Brinkman layer, here the jpdf value for each ๐๐/๐๐ฆ and ๐ฃ combination is
calculated at every ๐ฆ grid point in the region bounded by these two elevations and then averaged.
The jpdf values shown in Figure 3.7(e,f) are normalized such that the area integral in each of
these two plots is one. Both quantities are normalized by ๐ฟ and the bulk velocity (๐ข ๐ ) in the
โซ๐ฟ
water column; ๐ข ๐ is calculated as 1/๐ฟ 0 โจ๐ขโฉ๐๐ฆ. Here ๐ข ๐ instead of ๐ข ๐ is used for normalization
to highlight the effect of roughness geometry. For both cases, the vertical pressure gradients are
61
indeed negatively correlated with the vertical flux. Comparing the two cases, the random case
displays much wider distribution of ๐ฃ and narrower distribution of ๐ฆ-gradient of pressure. This
means that ๐ฃ is more sensitive in this case to changes in ๐๐/๐๐ฆ, probably because of the higher
effective permeability near the interface than for the regular case (Shen et al., 2020). Notice also
that in both cases the magnitude of ๐ฃ < 0 reaches much higher values than the magnitude of ๐ฃ > 0,
again indicating faster fluid parcel entry into the bed than exit from the bed.
Figure 3.7 indicates that, for both cases, the interfacial pressure distribution is closely correlated
with local ๐ฃ and therefore is probably the driving force of ๐ฃ transport near the interface. Forces of
the Reynolds and viscous stresses play secondary roles. In other words, the differences in transit
time characteristics between Cases 1 and 2 are probably due to the different interfacial pressure
distributions caused by the roughness geometries obstructing the surface flow in different ways.
3.3.3 Effect of molecular diffusion on TTD
The transit time calculation using Equation (3.3) considers the movement of a fluid parcel. On the
other hand, the transit time of a solute is affected also by molecular diffusion. An exact estimation
of the solute transit time would require tracking either the instantaneous location of a tracer particle
(during the DNS simulation) or the concentration variation which yields a breakthrough curve.
These tasks are beyond the scope of the present work. An alternative approach is to use the random
walk method (Kinzelbach, 1988) to incorporate the effect of molecular diffusion into Equation (3.3)
and this approach is used here.
The random walk particle tracking is based on a set of stochastic differential equations. Follow-
ing Sun et al. (2015) who carried out pore-scale random walk particle tracking for a 3D domain,
we solve
โ๏ธ
๐ฅ๐ (๐ก + ฮ๐ก) = ๐ฅ๐ (๐ก) + ๐ข๐ (๐ฅ๐ , ๐ก)ฮ๐ก + ๐๐ 2๐ท ๐ ฮ๐ก, (3.5)
where ฮ๐ก is a small time differential, ๐๐ is a random number of standard Gaussian distribution
with unit variance and zero mean, and ๐ท ๐ is the molecular diffusion coefficient. Equation (3.5)
is a generalization of Equation (3.3). Compared to the approach of Kinzelbach (1988), here the
62
dispersion terms are ignored by setting the longitudinal and transverse dispersivity lengths to zeros.
This is because the DNS resolves the time-mean flow details. Note also that the diffusion due to the
turbulent velocities ๐ขโฒ๐ is assumed to be much smaller than that due to e ๐ข๐ because of much weaker
Reynolds stresses compared to the dispersive stresses below the interface. The smoothness of
the results of Equation (3.5) requires the step of particle advection to be a fraction of the grid cell
size; Kinzelbach (1988) recommended five to ten steps per cell. Here, the ฮ๐ก value corresponds to
five steps per grid cell at a convection velocity of โจ๐ขโฉ below the Brinkman layer. ๐ท ๐ /๐ข ๐ ๐ฟ takes a
value of 2.5 ร 10โ6 to achieve a Schmidt number (๐๐ = ๐/๐ท ๐ ) of around 1000, representative of
dissolved oxygen in water.
Figure 3.8: Cumulative backward transit time distribution (a) and its pdf (b) based on mean
advection only ( ) and based on both mean advection and molecular diffusion ( ), for Case
1 (black) and Case 2 (red) at ๐ฆ 1 = ๐ฆ 2 = 0. In (b), (blue) shows the fitted power law with slope
indicated.
The backward TTD and its cumulative distribution are shown in Figure 3.8 for ๐ฆ 1 = ๐ฆ 2 = 0.
Various TTD characteristics are tabulated in Table 3.2 under columns โCase 1 (MD)โ and โCase
2 (MD)โ. The TTD power-law decay slope is milder for both cases, with the slope of around
โ1.1, compared to โ1.3 considering advection only. Changing to a deeper elevation, such as
๐ฆ 1 = ๐ฆ 2 = โ3๐ท, does not affect the power-law slope. For both cases, the addition of molecular
63
diffusion increases significantly the probability of transit times being longer than ๐(102 )๐ฟ/๐ข ๐ , while
early time values are only slightly affected. This increase in transit times is attributed to the fact that
solute particles are displaced from the time-mean paths by molecular diffusion, resulting in paths
with greater lengths. The addition of molecular diffusion does not remove the difference brought
by a change in the roughness geometry. For a particle release elevation below the interface region
(i.e. ๐ฆ 1 = ๐ฆ 2 = โ3๐ท) which tracks a large fraction of deep subsurface paths, ๐95 normalized by
๐ฟ/๐ข ๐ is significantly longer in Case 2 (2.17 ร 106 ) than in Case 1 (8.35 ร 105 ).
For a 2D sinusoidal bedform analyzed using an advection-dispersion model, Bottacin-Busolin
and Marion (2010) also observed a milder slope (โ0.5) of the residence time function when
both mean advection and mechanical dispersion were accounted for, compared to the slope with
advection only (โ1.0). According to Equation (3.2) and the as shown in the Appendix B, the
present slopes of โ ๐โ๐ shown in Figure 3.8(b) indicate power-law decay slopes of โ0.1 and โ0.3 of
the same residence time function, with or without molecular diffusion, respectively. The difference
between the slopes observed herein and those in the study of Bottacin-Busolin and Marion (2010) is
possibly due to the three-dimensionality of the present flow and the capability of DNS in capturing
the actual velocities near the interface and within the bed.
3.4 Conclusions
Results of transit times are reported for turbulent open-channel flows over flat sediment beds at
๐
๐ ๐พ = 2.6 with two different bed roughnesses consisting of either regularly (Case 1) or randomly
(Case 2) positioned sediment grains at the bed surface. The transit time calculations are based on
DNS simulations similar to the ones reported in Shen et al. (2020). The transit times are calculated
using the particle tracking method, considering (i) advection by time-mean velocity only and (ii)
both time-mean advection and molecular diffusion.
Despite the absence of bedforms, both bed roughnesses generated a large spectrum of transit
times, whose characteristics are similar to those generated by bedforms. The TTDs show power-law
tails at late times for both cases; the time range and slope of the tail vary only weakly with the
64
roughness geometry. For both cases, the power-law slopes range from โ1.3 considering the mean
advection only, and around โ1.1 considering both the mean advection and the molecular diffusion.
We also show that, with pore-resolved data, the definition of a surface-subsurface interface becomes
obscure; the calculated transit times vary with the elevations chosen to mark the โentryโ or โexitโ of
particles.
Near the interface, roughness leads to regions with local pressure maxima or minima due to its
effect on the near-bed turbulence. These regions are found to correlate closely with the distribution
of interfacial bed-normal velocity. This suggests that locations of subsurface flow path entries and
exits are mainly determined by the interfacial pressure, while being less sensitive to forces due
to the Reynolds and viscous stresses. In other words, the advective pumping mechanism may be
used to describe the subsurface flow induced by roughness on a macroscopically flat bed, similar
to those induced by bedforms.
Compared to the regular roughness, the random roughness leads to longer transit times in the
TTD power-law tail associated with deep subsurface pathlines. Specifically, the pure-advection
particle tracking in Case 2 yields 18% and 91% larger transit times of the 50th and the 95th
percentiles, respectively, than in Case 1 when ๐ฆ 1 = ๐ฆ 2 = โ3๐ท. This difference occurs because
the random roughness induces more intense interfacial pressure variation (due to grain clusters
forming protuberances larger than a single grain), with a wider range of horizontal length scales.
Consequently, the subsurface flow paths in Case 2 can reach farther and deeper than in Case 1.
With the addition of molecular diffusion in the transit time calculation, the transit time increases
(especially at late times) for both cases, as expected. The differences between the two cases remain
significant.
This work demonstrates that grain-scale roughness and its geometry are important for hyporheic
exchange, as shown by the subsurface flow characteristics and transit times. The effect appears to be
fundamentally similar to how a larger-scale bedform affects the exchange. These findings would not
have been possible if a pore-unresolved simulation model (such as the Darcy flow model) was used
instead. In addition, DNS data similar to the present data can be used for closure developments,
65
such as the parameterization of advection due to e ๐ข๐ and ๐ขโฒ๐ in TTD or RTD prediction based on a
Darcy flow. Additional factors that are likely to affect the overall transit time of solutes include (i)
grain-scale solute concentration variation near the bed surface and (ii) turbulent diffusion, which
is strong in the Brinkman layer and is likely to affect ๐ values calculated with ๐ฆ 1 or ๐ฆ 2 close to 0.
The importance of these factors and their dependence on roughness details need to be addressed in
future work. In addition, scenarios that include both roughness and bedforms need to be explored
for possible interactions between the effects of the two.
66
CHAPTER 4
PORE-RESOLVED DIRECT NUMERICAL SIMULATIONS OF HYPORHEIC
EXCHANGE INDUCED BY BEDFORMS AND BED ROUGHNESS
4.1 Introduction
Bedforms are the most important source of flow resistance at the local scale in river channels.
The word โBedformโ are used herein refers to ripples or dunes as described by Vanoni (2006). In
sand-bedded rivers, form resistance is primarily generated by bedforms. The turbulent motions and
flow separation generated by bedforms are responsible for transfer of mean flow energy to turbulent
kinetic energy (Venditti, 2013).
The objective of this Chapter is to identify to what extent the grain-scale bed roughness affects
the hyporheic exchange with the presence of bedform. To this end, pore-resolved DNS simulations
of turbulent half-channel flows bounded by the same permeable dune-shape bedform with two
different bed roughnesses are carried out and compared. The effects of grain-scale roughness and
its texture on macroscopic exchange quantities including the subsurface flow paths and transit time
characteristics are documented. A forward particle tracking method is used to calculate the transit
times, considering the advection by time-mean velocity only, same as in Chapter 3.
The chapter is organized as follows. Section 4.2 describes the numerical techniques used for
DNS and transit time calculation, the problem set-up, and the syntheses of bedforms with two
different roughnesses superimposed on them. Section 4.3 compares the effects of bed roughness
geometry on velocities, pressure field and transit times for the two cases examined and analyzes the
effects of roughness with the presence of bedform, followed by conclusions in Section 4.4.
67
4.2 Methodology
4.2.1 Configurations and parameters
Fully developed open-channel flows are simulated using the same DNS method mentioned in
previous chapters. The simulation domain includes both the surface and sub-surface flows (Fig-
ure 4.1(a,b)). ๐ฟ is the thickness of the turbulent boundary layer (i.e., open-channel height) measured
from the trough of the bedform; ๐ป๐ is the sediment depth measured from the bedform trough to
the bottom boundary of the simulation domain. ๐ป is the bedform height measured from its trough
to its crest. The bedform wavelength is denoted by ๐. The ๐ฅ location of the crest is ๐ฟ ๐ = ๐/2.
Symmetric boundary conditions are applied at both top and bottom boundaries of the simulation
domain. Periodic conditions are applied at ๐ฅ and ๐ง boundaries. A constant mean pressure gradient
is used to drive the flow. The elevation of ๐ฆ = 0 is set at the bedform trough elevation. The same
double-averaging (DA) decomposition of an instantaneous flow variable as that used in previous
Figure 4.1: Simulation domain and synthesized bedforms with (a) regular and (b) random
roughnesses, colored by ๐ฆ.
68
Roughness ๐ท/๐ฟ ๐ป๐ /๐ฟ ๐ป/๐ฟ ๐/๐ฟ ๐ฟ ๐ /๐ (๐ฟ ๐ฅ , ๐ฟ ๐ง )/๐ฟ (ฮ๐ฅ + , ฮ๐ฆ +min , ฮ๐ง + )
Case 1 Regular 0.05 2 0.15 3 0.5 (3, 1) (3.7, 0.3, 4.1)
Case 2 Random 0.05 2 0.15 3 0.5 (3, 1) (3.7, 0.3, 4.1)
Table 4.1: Summary of parameters. ๐ท is the grain diameter; ๐ฟ ๐ฅ๐ is the simulation domain size in
๐ฅ๐ ; ๐
๐ ๐พ = 2.5, ๐
๐ ๐ = 1580, ๐ท + = 79. ฮ๐ฅ + , ฮ๐ฆ + , and ฮ๐ง + are DNS grid sizes in ๐ฅ, ๐ฆ and ๐ง,
respectively, normalized using the viscous length scale ๐/๐ข ๐ .
chapters are applied. Given the streamwise dependence of this flow, a modified double-averaging
decomposition is used, with the spatial averaging performed in ๐ง only, instead of averaging in both
๐ฅ and ๐ง as in the previous chapters.
The simulation and case parameters are listed in Table 4.1. Two cases are considered, one
with regular roughness (Case 1) in which the arrangement of grains in the uppermost layer is
regular and the other with random arrangements of uppermost-layer grains (Case 2). In both
cases, the same porosity of 0.4 is imposed in the bulk of the sediment. The friction Reynolds
number and permeability Reynolds number are the same for both cases: ๐
๐ ๐ = ๐ฟ๐ข ๐ /๐ = 1580 and
โ
๐
๐ ๐พ = ๐พ๐ข ๐ /๐ = 2.5, where ๐พ is the permeability. The ๐
๐ ๐พ is the same as that used in previous
chapters. The grain diameter in wall units is ๐ท + = 79. ๐ข ๐ is obtained as the average of the local
๐ข ๐ , calculated from the maximum magnitude of the local ๐ฆ profile of the total shear stress, which
is the sum of the viscous shear stress, the Reynolds shear stress and the form-induced shear stress.
The number of grid points are 1280, 1147 and 384 in ๐ฅ, ๐ฆ and ๐ง, respectively. The total simulation
time used for data collection is around 10 large-eddy turn-over times (LETOTs), defined as ๐ฟ/๐ข ๐ .
The porous bedforms are immobile and modeled as closely packed mono-disperse hard spheres.
The locations of the grains are determined based on molecular dynamics simulations (Plimpton,
1995; Shen et al., 2020). Two bedforms with the same macroscopic dune geometry but different
roughnesses at the uppermost layer are synthesized and shown in Figure 4.1. One is the โregularโ
case (Case 1) formed by regular distribution in (๐ฅ, ๐ง) of grains and the other is the โrandomโ case
(Case 2) formed by random grain distribution of grains in both (๐ฅ, ๐ง) and ๐ฆ. Both roughnesses are
statistically similar to those imposed on a flat bed in the previous chapters. Below the bed surface,
the grain distribution is random in both cases.
69
Figures 4.2 compares the bed surfaces in Cases 1 and 2, showing the difference in grain
arrangement at the top of the bedform.
Figure 4.2: Roughness surfaces colored by ๐ฆ (a,b) and local height fluctuations measured from
the bedform surface (c,d) in Cases 1 (a,c) and Case 2 (b,d). (e) One-dimensional power spectral
densities ๐ธ ๐ of the local height fluctuations for Case 1 (regular roughness, ) and Case 2
(random roughness, , red). ๐ ๐๐๐ is the root-mean-square (rms) of interface height
fluctuations.
To quantify the difference between the two roughness geometries, the one dimensional (1D)
70
power spectra (๐ธ ๐ ) of local height fluctuations measured from the bedform surface are compared
in Figure 4.2(e) as functions of the streamwise wavenumber (๐
1 = 1/๐ 1 , where ๐ 1 is the streamwise
wavelength). Similar spectral distributions as those seen in Chapter 3 are observed here, expectedly.
For the random case the white noise regime applies to ๐ 1 /๐ท larger than around 2.0. For smaller ๐ 1 ,
๐ธ ๐ shows a power-law decay with a slope of around โ2.4, which is the same as shown in Chapter
2. In comparison, for the regular roughness the white noise region covers wavenumbers up to ๐ท.
4.2.2 Transit time calculation based on particle tracking
The same definition and calculation of transit times are used as in Chapter 3 and are not repeated
here. The only difference is the particle release and exit locations. Different from a constant height
elevation for particle enter and exit used in Chapter 3, here we release particles on the streamwise-
varying bedform interface and track particles until they exit from the interface, as shown in Fig. 4.3.
Particles residing in the sediment without exiting or particles exiting the sediment elsewhere other
than the interface are not accounted for in the transit time calculation. The particles are seeded
uniformly on the bedform surface. In total, around ๐(106 ) particles are released. To test the result
Figure 4.3: Conceptual sketch to show the enter and exit locations of tracked particles.
71
convergence, we used twice the number of released particles, which led to less than 1% difference
in the fitted TTD slope or the median transit time value.
Figure 4.4: (a) Time- and spanwise-averaged streamwise velocity โจ๐ขโฉ +๐ง for cases with regular
( ) and random ( , red) roughnesses. โจ๐ขโฉ +๐ง contours for (b) regular and (c) random
roughness.
4.3 Results
In this section, Cases 1 and 2 are compared to understand the effects of roughness on the hyporheic
exchange with the presence of bedform, specifically the effects on turbulent flow and pressure
distributions along the SWI and how they influence the mass transport in the sediment.
72
4.3.1 Effects of roughness geometry on velocity
The ๐ฆ profiles of the time- and spanwise-averaged streamwise velocity normalized by average ๐ข ๐
are compared in Fig. 4.4 (a). These profiles are qualitatively similar with previous observations
of Chen et al. (2015) and Kirca et al. (2020), which show that the flow is accelerated towards the
bedform crest, due to the streamline convergence, and decelerated after the crest. The random
roughness leads to a lower bulk velocity compared to the regular one and, consequently, results in
a 17% higher friction coefficient,
2
๐ข๐
๐ถ๐ = 2 , (4.1)
๐๐ฟ
of 0.0109 compared to 0.0093 in the regular-roughness case. Here, ๐๐ฟ is the DA velocity at half
channel height. The higher ๐ถ ๐ is due to the larger effective length scale of the random roughness
enhancing the local wall friction on the bedform surface, consistent with the observations made on
a macroscopically flat bed for the same roughnesses in Chapter 2.
Profiles of the viscous shear stress normalized by local wall units are shown in Fig. 4.5(a).
Despite the adverse-pressure-gradient (APG) region on the lee side and a sharp edge of the bedform
at its crest, no mean flow separation (or change of sign of the ๐ฆ-gradient of streamwise mean
velocity) is observed. More interestingly, the inflection point of the mean velocity profile stays
at the permeable wall, as opposed to the presence of an inflection point inside an APG boundary
layer bounded by an impermeable wall (as shown by the similarity solution of the Falkner-Skan
equation or by extensive experimental and numerical evidences in turbulence boundary layers).
This indicates that a permeable wall promotes inner-outer interaction in a turbulent boundary layer
and, consequently, reduces the likelihood of APG-induced flow separation.
The Reynolds shear stress profiles (Fig. 4.5(b)) show deeper penetration of turbulence into the
bed in the random case, which indicates stronger wall-normal mean momentum transfer at the
interface. Different from the viscous shear stress reaching the maximum values on the bed surface,
the maximum magnitudes of the Reynolds shear stress are reached above the wall, farther from
the wall in the APG region, consistent with observations made for impermeable-wall APG flows.
In addition, the streamwise maximum of the Reynolds shear stress occurs on the lee side, not
73
corresponding to the ๐ฅ location of the maximum of the ๐ฆ-gradient of the streamwise mean velocity.
This indicates that the turbulence are not generated by the mean shear alone; other mechanisms,
such as the shear production due to mean strain rate at the roughness scale (i.e. wake production
(Raupach and Shaw, 1982; Yuan and Piomelli, 2014c; Yuan and Jouybari, 2018; Mangavelli et al.,
2021)), may be important.
ยฏ ๐ง+ /๐๐ฆ + and (b) Reynolds shear stresses
Figure 4.5: Profiles of (a) viscous shear stresses โ๐โจ๐ขโฉ
โโจ๐ขโฒ๐ฃ โฒโฉ๐ง+ for cases with regular ( ) and random ( , red) roughnesses.
4.3.2 Effects of roughness geometry on pressure field
The distribution of the time- and spanwise-averaged pressure across the half channel and the
permeable bed is compared in Fig. 4.6 (a,b). In both cases the bedform introduces favorable
streamwise pressure gradients on the stoss side and adverse gradients on the lee side. Similar to
74
what was observed by Cardenas and Wilson (2007) using RANS simulation of surface flows, the
minimum pressure is located at the crest of bedform. The pressure distribution along the SWI
drives the flow in the sediment and controls the surface-subsurface dynamics exchange and mass
transport.
Figure 4.6: Time- and spanwise-averaged pressure field for cases with (a) regular and (b) random
roughnesses. Pressure distributions along the SWI normalized using bulk velocity (c) or friction
velocity (d) for both cases.
The interfacial pressure non-dimensionalized using either ๐ข ๐ or average ๐ข ๐ is shown in Fig. 4.6
(c) and (d). When normalized using ๐ข ๐ (Fig. 4.6 (d)) the pressure variation induced by the random
roughness is significantly weaker (by almost a half in magnitude) than that generated by the regular
one. However, the profiles almost collapse when the pressure is non-dimensionalized using ๐ข ๐
instead (Fig. 4.6 (c)). This implies that the difference of pressure variations at the SWI between
the two cases is mostly due to the higher ๐ถ ๐ in the random case. The interfacial pressure variation
reflects the pressure drag induced by the dune, which is a quadratic function of the speed of flow
75
passing it. As the random roughness provides stronger flow resistance than the regular one, it
reduces the speed of flow passing the dune and thus reduces the pressure variation. Although the
size of roughness is much smaller compared to the size of bedform, it significantly influences the
pressure distribution along the interface.
Figure 4.7: Distribution of the Clauser parameter ๐ฝ in the two cases. The black dash line shows
the bedform interface geometry.
The strength of pressure gradient can be quantified by the Clauser parameter (Clauser, 1956),
๐ฝ(๐ฅ) = (๐ฟโ /๐๐ค ) (๐๐โ /๐๐ฅ), where ๐๐โ /๐๐ฅ is the streamwise pressure gradient at the free-surface
(or half channel height), ๐ฟโ is the displacement thickness, and ๐๐ค is the wall shear stress. Here, ๐๐ค
ยฏ ๐ง
is set to the maximum value of local total shear stress ๐ = ๐ ๐โจ๐๐ฆ
๐ขโฉ
โ โจ๐ขโฒ๐ฃ โฒโฉ๐ง โ โจ๐ขห ๐ฃห โฉ๐ง . The distributions
of ๐ฝ in both cases are shown in Figure 4.7, which shows a sinusoidal variation with negative values
(i.e. favorable pressure gradients (FPG) and flow acceleration) on the stoss side and positive ones
(i.e. APG and flow deceleration) on the lee side. The regular roughness case yields larger ๐ฝ
magnitudes with both signs than does the random case, consistent with the discussions above. The
๐ฝ values on the lee side indicates that the turbulence is subjected to mild APGs (Aubertine and
Eaton, 2005). A different bedform shape that induces stronger APG is needed to test the effect of
wall permeability on flow separation.
Compared to that on a flat bed, the effect of roughness geometry on the pressure distribution on a
bedform is significantly different. On a flat bed, the roughness protuberance creates local stagnation
76
Interface ๐50 ๐95 Slope
Case 1 Regular 3.00 385 โ1.52
Case 2 Random 0.0854 7.47 โ1.54
Table 4.2: Various characteristics of transit times calculated based on mean advection only.
Transit times are normalized by ๐ฟ/๐ข ๐ . Note. ๐50 and ๐95 represent 50th (median) and 95th
percentiles of the TTD distributions.
point upstream of the protuberance and low pressure region downstream, inducing roughness-scale
pressure heterogeneity at the interface. As shown earlier in Figure 2.24, the more intense pressure
variations are observed at the random interface of flat bed due to larger roughness scales. However,
with the presence of bedform, the higher wall friction due to the random roughness reduces the
bulk velocity and the macroscopic pressure variation on the bedform. Meanwhile, the macroscopic
pressure variation induced by bedform masks the roughness-scale pressure variation induced by
the roughness. In other words, the effect of roughness on the pressure field appears global when a
bedform is present, while, on a flat bed, the roughness effect is local.
4.3.3 Effects of roughness geometry on transit times
The effect of bed roughness on modulating the interfacial pressure has significant implications
in the exchange of mass (e.g. water and solutes) across the interface, which can be quantified,
for examples, using transit time distribution. Table 4.2 compares characteristics of transit times
calculated based on mean advection only. Particles are released at the bedform interface and tracked
until they exit from the interface, as shown in Figure 4.3. Here, ๐50 represents the central tendency
of the distribution, while ๐95 is considered as a characteristic ๐ associated with the rare, long-transit-
time subsurface paths. The comparison shows that the ๐50 and ๐95 of the transit time distribution
in the regular case are significantly higher than in the random case. However, a constant slope of
the power-law tail of around โ1.5 is observed for both cases. This is similar with the observations
of Sawyer and Cardenas (2009) that the tails of residence time distributions follow a power law
regardless of permeability structure near the SWI. An important conclusion that follows is that the
77
difference in grain-scale roughness texture on top of a bedform extends deep into the sediment and
significantly modifies the transit time distributions.
Figure 4.8: Cumulative backward transit time distributions (a) and their pdf distributions (b) for
regular ( ) and random ( , red) cases. In (a), the horizontal dashed lines ( , blue) mark
the 50th (median) and 95th percentiles. In (b), the dashed blue lines ( , blue) show the fitted
power law with slopes indicated. Subsurface flow paths in the (c) regular and (d) random cases.
To explain such difference in the TTD characteristics, Figure 4.8(c,d) compares the subsurface
paths obtained using the particle tracking approach. For both cases, multiscale subsurface flow
paths are induced, with the longest ones characterized by roughly a half of the bedform length scale
๐. The general pattern is consistent with observations in existing numerical and experimental studies
of subsurface flows induced by bedforms based on advective pumping (Cardenas and Wilson, 2007;
Elliott and Brooks, 1997b). Here, the regular roughness is shown to induce a larger fraction of
78
deep-reaching paths and, consequently, overall longer transit times. This is probably due to the
higher-magnitude interfacial pressure variation observed in Fig. 4.6 for this case. Figure 4.9 shows
the histogram of the number of pathlines that reaches a depth of each given ๐ฆ value. The area of
the integral equals one. It is evident that a larger fraction of pathlines in the regular case reaches
deep in the bed than in the random case.
Figure 4.9: Histogram of the number of flow paths reaching a depth of a given ๐ฆ value for the two
cases. Values are normalized such that the area integral is one.
The comparison with TTD results calculated for flat beds (Table 3.2) shows that the roughness
influences the mass exchange differently depending on whether a bedform is present. On a flat
bed, the random roughness induces slightly larger ๐50 measured at the interface than the regular
roughness. However, with the presence of bedform, the regular roughness case induces a much
longer transit times. This is due to the different roles played by roughness. The pressure difference
here is much larger than on a flat bed. The results show that grain-scale roughness on top of a
79
bedform modulates the pumping effects of bedform in a way that significantly modifies the transit
times.
4.4 Conclusions
Pore resolved DNSs are carried out over immobile bedforms at ๐
๐ ๐พ = 2.5 and ๐
๐ ๐ = 1580 in order
to investigate the effects of two roughnesses, regular and random, on the flow and transit times with
the presence of bedforms. The time-averaged velocity profiles show that the roughness affects the
flow velocity near the SWI, resulting in different wall friction depending on the roughness texture
(similar to what happens in impermeable-bed flows), viscous shear stress and shear penetration
depth.
The two roughness textures also induces different pressure variation along the bedform surface.
On the bedform with regular roughness case, a larger pressure variation is observed than in the
random roughness case. This is contrary to what has been observed on a flat bed (in Chapter 3)
where the random roughness induces more intense pressure variation of the roughness scale. This
discrepancy between flat-bed and bedform cases occurs because of the different roles played by
roughness. On flat bed, roughness protuberance creates roughness-scale pressure heterogeneity,
while on a bedform, the roughness also modulates the pressure variation of the bedform scale.
In addition, contrary to the flat-bed flows where the random roughness augments mass storage
in the bed by inducing longer and deeper flow paths and longer transit times, with the presence of
bedform the opposite is observed: the random roughness damps the pumping effect of the bedform
and reduces storage. The observations show that grain-scale roughness superimposed on a bedform
does not simply act as a smaller-scale โbedformโ with a lesser effect on pumping. Instead, roughness
modulates significantly the bedfrom-induced pumping exchange, by changing the macroscopic
pressure distribution along SWI. This work demonstrates that grain-scale roughness on top of a
bedform should not be ignored when studying hyporheic exchange, as shown by its influence on
the pressure variation and transit times.
80
CHAPTER 5
CONCLUSIONS
5.1 Summary
In this dissertation, three-dimensional pore-resolved direct numerical simulation are carried out to
elucidate multiscale flow physics regarding mass and momentum transport across water-sediment
interfaces. River-bed sediments representative of realistic ones are synthesized based on molecular
dynamics simulations. An immersed boundary method is used to impose the boundary conditions
of fluid regime in a complex sediment geometry. The methods are validated by comparing with
experimental results. A particle tracking method is implemented to calculate the transit time
distributions in the sediment, with or without the contribution of molecular diffusion. These
numerical methods have been applied to understand the effects of roughness on flows bounded by
flat beds and bedforms. The most important findings and their implications are summarized below.
Bed roughness plays a key role in inducing hyporheic exchange on both flat beds and bedforms.
On a flat bed, small bed morphological features on hyporheic exchange should not be ignored,
as they significantly affect the dynamics of turbulence, time-mean flow, subsurface pathlines and
transit time distributions. Specifically, the larger roughness length scales brought by the random
interface results in more intense mixing, pressure variations, more isotropic form-induced stress
tensor, significantly higher wall-normal hyporheic flux across SWI and deeper-reaching hyporheic
flow paths. These differences explain the dependence of transit times on the roughness geometry.
The results also demonstrate that the grain roughness on a flat bed induces hyporheic exchange
predominantly through the pumping mechanism, similar to the effect of a bedform.
With the presence of bedform, grain roughness still plays an important role in modifying the
wall friction, interfacial pressure distribution and velocity profiles, as well as penetration depth.
The large-scale variation of the interfacial pressure is induced by the bedform geometry, while the
magnitude of its variation is significantly modulated by the texture of roughness superimposed on
81
the bedform. Roughness still leads to local pumping of interfacial flow; but a more significant
part of its effect is the augmentation of hydrodynamic drag, which slows down the core flow and
consequently reduces the pressure variation on the bedform. As a result, the transit times are more
sensitive to the roughness texture when a bedform is present than when it is not. These results
demonstrate that the effect of bed roughnesses interacts nonlinearly with that of larger bed scales
such as bedforms.
5.2 Future work
Pore-resolved simulations such as the DNS herein offer an unprecedented amount of information
useful in direct characterization of pore-scale drivers of the flow. In future work, other contributing
factors of hyporheic exchange can be further studied with the present methodologies. These factors
include Reynolds numbers, bulk permeability, spatial heterogeneity of bed permeability, grain-
diameter distributions, bedform configurations such as asymmetric dunes or those with high slope,
and moving or deformable beds. DNS data of surface-subsurface exchange within this parameter
space are needed to fully understand the hyporheic exchange and validate pore-unresolved models
of the exchange. In addition, direct simulations of reactions and microbial biomass grow/decay
may be used to understand pore-scale physics in an array of increasingly important topics, such as
transport of nutrients/pollutants and greenhouse-gas production in river beds.
82
APPENDICES
83
APPENDIX A
DETAILS OF POROUS-BED SYNTHESES
Both the porous beds with random and regular interface are generated by simulating the pouring
of hard spheres of the same diameter, ๐ท, into the MD simulation domain. The particles are subject
to the gravitational force in the ๐ฆ direction and the forces between particles determined by the
Hookean-style granular potential with history effects. The top and bottom boundaries are fixed,
where the repulsive boundary condition is imposed; the periodic boundary condition is imposed
on the ๐ฅ and ๐ง boundaries.
For the bed with the random interface packing, the particles are released at each time step from
the top boundary and move towards the bottom boundary due to imposed gravity (Figure A.1). The
simulation ends as the sediment-bed domain is filled with randomly packed spheres. For the bed
with the regular interface, special treatment is required at the bottom boundary. Specifically, a layer
of regularly packed spheres are initially positioned on the bottom boundary (Figure A.2a); a second
Repulsive
Periodic
Periodic
Periodic
Repulsive
Periodic
(a) (b)
Figure A.1: Bed synthesis with random interface packing.
(a) (b) (c) (d) (e)
Figure A.2: Bed synthesis with regular interface packing.
84
layer of spheres are also manually positioned following hexagonal close packing (Figure A.2b).
However, as is, the first two layers of spheres lead to locally low value of ๐ at such elevation.
To maintain the running average of ๐ (๐ฆ) (with an averaging window size of ๐ท) as almost a
constant value of ๐ ๐๐ฃ๐ , the second layer of spheres are shifted away from the wall by a small
amount determined by the targeted ๐ ๐๐ฃ๐ value (of an order of ๐(0.1๐ท) with 3% fluctuations added,
Figure A.2(c)). Then, the similar procedure as in Figure A.1 is employed to fill the rest of the porous
bed with randomly packed spheres. In the end, the MD simulation domain is flipped upside-down
to yield a regularly-packed uppermost layer. Once the particle positions are determined using the
aforementioned procedure, the particles are overlaid with the DNS grid and the volume of fluid for
each grid cell is determined for the immersed boundary method of the fluid solver.
85
APPENDIX B
SUPPORT INFORMATION FOR CHAPTER 3
Calculation of the residence time distribution
Integrating the master equation (Equation (3.2)) for a steady system with the boundary condition
of zero age for water parcels entering via inflow, i.e.,
๐ฝ
๐ ๐ (0) = , (B.1)
๐
and noting that ๐ and ๐ฝ = ๐ are constants for a steady flow, one obtains
โซ ๐ โซ ๐
๐๐ ๐ (๐) = โ
๐ โ
๐โ๐ (๐)๐๐. (B.2)
๐ ๐ 0
Figure B.1: Residence time distributions (RTD) normalized by ๐/๐ for ๐ฆ 1 = ๐ฆ 2 = 0 ( ) and
โ3๐ท ( 6 7
). For ๐๐ข ๐ /๐ฟ > ๐(10 ) when ๐ฆ 1 = ๐ฆ 2 = 0 or ๐(10 ) when ๐ฆ 1 = ๐ฆ 2 = โ3๐ท, the RTD
rapidly decreases; this is not physical but an artifact associated with the finite domain size in ๐ฅ
(limited to 2๐ฟ ๐ฅ ) imposed for particle tracking as the tracking procedure cannot be indefinitely long.
86
โโ โโ
Since ๐ ๐ ๐ = โ ๐โ๐ ๐๐ and ๐ ๐ (0) = 0 (as the minimal age is 0), this yields
๐ h โ
โ i
๐ ๐ (๐) = 1 โ ๐ ๐ (๐) . (B.3)
๐
Equation (B.3) states that the total volume of water storage with age ๐, which equal to the sum
of the volume of stored water that is aging to ๐ + ๐๐ and the volume of water leaving the bed at
โโ
age ๐, is the same as the total water volume leaving the bed at age ๐ or higher (i.e., ๐(1 โ ๐ ๐ )).
This is a result of the time-invariant water age distribution in a steady system. The calculated ๐ ๐
distributions for the regular and random cases, with ๐ฆ 1 = ๐ฆ 2 = 0 or โ3๐ท, are shown in Figure B.1.
Validation of particle tracking method
To verify the implementation of the particle tracking method, the subsurface flow induced by
a 2D sinusoidal pressure field applied at the interface is considered, which corresponds to the
stationary bedform case studied in EB97. Note that we do not discuss the results of the residence
time function from the DNS data. This quantity is used for validation purpose in this section only.
An analytical solution of the Darcy velocity was derived by EB97 as
โ โ
๐ข โ = โ cos ๐ฅ โ ๐ ๐ฆ , ๐๐๐ ๐ฃ โ = โ sin ๐ฅ โ ๐ ๐ฆ , (B.4)
where ๐ is the wavelength of the imposed interface pressure distribution and โ denotes normalization
using the maximum velocity amplitude at the interface and ๐/(2๐). Following EB97, consider a
pulse of solute entering the bed at a location ๐ฅยฎ. Assume that the flow is steady. Denote the fraction
of solute molecules that remains in the bed after an elapsed time ๐ as ๐
(ยฎ ๐ฅ , ๐) = 1 before
๐ฅ , ๐). ๐
(ยฎ
the pulse of solute leaves the bed and ๐
(ยฎ ๐ฅ , ๐) = 0 after it leaves the bed. The overall residence time
function averaged across the bed surface, ๐
(๐), is obtained by averaging ๐
at all surface locations,
weighted by the local bed-normal fluid volumetric flux at the bed surface, ๐(ยฎ ๐ฅ ),
โซ
(๐ฅ,๐ง)
๐(๐ฅ, ๐ง)๐
(๐ฅ, ๐ง, ๐) ๐๐ฅ๐๐ง
๐
(๐) = โซ , (B.5)
(๐ฅ,๐ง)
๐(๐ฅ, ๐ง) ๐๐ฅ๐๐ง
โซ
where (๐ฅ,๐ง) ๐๐ฅ๐๐ง denotes integration throughout the fluid area at the bed surface. According to
EB97, the analytical solution of ๐
, given in implicit form, is
2 cosโ1 ๐
๐โ = . (B.6)
๐
87
Equation (B.4) is used as the advection velocity for particle tracking. The ๐ฅ โ domain size
covers two wavelengths, ๐ฅ โ โ [โ๐/2, 3๐/2], while the ๐ฆ โ domain size is varied from [โ๐/2, 0]
to [โ2๐, 0] to study the convergence of the particle tracking method. Uniform grids in ๐ฅ and ๐ฆ
are used to discretize the velocity fields in Equation (B.4). The spatial resolution is also varied,
from 10 points to 50 points per wavelength. The particle-tracked ๐
is compared with the analytical
solution provided by Equation (B.6) in Figure B.2.
10 0 10 0
Slope -1
-1 -1
10 10
-2 -2
10 10
0 1 2
10 10 10 10 0 10 1 10 2
โค โค
AAAB6nicbZDLSsNAFIZP6q3WW9Slm8Ei1E1Jili7kYIblxXtBdpQJtNJO3QyCTMToYQ+ghsXirj1idz5Nk7TIN5+GPj4zzmcM78fc6a043xYhZXVtfWN4mZpa3tnd8/eP+ioKJGEtknEI9nzsaKcCdrWTHPaiyXFoc9p159eLerdeyoVi8SdnsXUC/FYsIARrI11W8GnQ7vsVJ1M6C+4OZQhV2tovw9GEUlCKjThWKm+68TaS7HUjHA6Lw0SRWNMpnhM+wYFDqny0uzUOToxzggFkTRPaJS53ydSHCo1C33TGWI9Ub9rC/O/Wj/RwYWXMhEnmgqyXBQkHOkILf6NRkxSovnMACaSmVsRmWCJiTbplLIQGpnQEupnOTTcrxA6tap7Xq3d1MrNyzyOIhzBMVTAhTo04Rpa0AYCY3iAJ3i2uPVovVivy9aClc8cwg9Zb5/HF43b
(a)
AAAB7XicbZBLSwMxFIXv1Fetr6pLN8EiiIsyI8XaXcWNywr2Ae1YMmmmjWaSIckIZSj4E9y4UMSt/8ed/8bpzCC+DgQ+zrkhN8cLOdPGtj+swsLi0vJKcbW0tr6xuVXe3uloGSlC20RyqXoe1pQzQduGGU57oaI48Djterfn87x7R5VmUlyZaUjdAI8F8xnBJrE6A4Oj66NhuWJX7VToLzg5VCBXa1h+H4wkiQIqDOFY675jh8aNsTKMcDorDSJNQ0xu8Zj2ExQ4oNqN021n6CBxRsiXKjnCoNT9fiPGgdbTwEsmA2wm+nc2N//L+pHxT92YiTAyVJDsIT/iyEg0/zoaMUWJ4dMEMFEs2RWRCVaYmKSgUlpCIxXKoF7LoeF8ldA5rjonVeeyVmme3Wd1FGEP9uEQHKhDEy6gBW0gcAMP8ATPlrQerRfrNRstWHmFu/BD1tsnon2P+Q== AAAB7XicbZBLSwMxFIXv1Fetr6pLN8EiiIsyI8XaXcWNywr2Ae1YMmmmjWaSIckIZSj4E9y4UMSt/8ed/8bpzCC+DgQ+zrkhN8cLOdPGtj+swsLi0vJKcbW0tr6xuVXe3uloGSlC20RyqXoe1pQzQduGGU57oaI48Djterfn87x7R5VmUlyZaUjdAI8F8xnBJrE6A4Oj66NhuWJX7VToLzg5VCBXa1h+H4wkiQIqDOFY675jh8aNsTKMcDorDSJNQ0xu8Zj2ExQ4oNqN021n6CBxRsiXKjnCoNT9fiPGgdbTwEsmA2wm+nc2N//L+pHxT92YiTAyVJDsIT/iyEg0/zoaMUWJ4dMEMFEs2RWRCVaYmKSgUlpCIxXKoF7LoeF8ldA5rjonVeeyVmme3Wd1FGEP9uEQHKhDEy6gBW0gcAMP8ATPlrQerRfrNRstWHmFu/BD1tsnon2P+Q==
โง โง
AAAB6nicbZDLSsNAFIZP6q3WW9Slm8Ei1E1Jili7kYIblxXtBdpQJtNJO3QyCTMToYQ+ghsXirj1idz5Nk7TIN5+GPj4zzmcM78fc6a043xYhZXVtfWN4mZpa3tnd8/eP+ioKJGEtknEI9nzsaKcCdrWTHPaiyXFoc9p159eLerdeyoVi8SdnsXUC/FYsIARrI11W/FPh3bZqTqZ0F9wcyhDrtbQfh+MIpKEVGjCsVJ914m1l2KpGeF0XhokisaYTPGY9g0KHFLlpdmpc3RinBEKImme0Chzv0+kOFRqFvqmM8R6on7XFuZ/tX6igwsvZSJONBVkuShIONIRWvwbjZikRPOZAUwkM7ciMsESE23SKWUhNDKhJdTPcmi4XyF0alX3vFq7qZWbl3kcRTiCY6iAC3VowjW0oA0ExvAAT/BscevRerFel60FK585hB+y3j4ByJyN3A==
(b)
Figure B.2: Validation of particle tracking against EB97, with various spatial resolution and a
sediment domain of ๐ฆ โ [โ๐/2, 0] (a) or ๐ฆ โ [โ๐, 0] (b).
For sufficiently deep sediment (with depth
ONR 331 of at least
Hydrodynamics one Review,
Program bedform wavelength)
16-23 September 2020 and sufficiently
fine spatial resolution (at least around 50 grid points per wavelength), the particle-tracked residence
time distribution matches well with the analytical solution. Figure B.2(a) shows that the values of
๐
at large times are significantly underestimated with a small sediment depth due to the absence
of deeper flow paths, while Figure B.2(b) shows that the large-time slope is slightly different at a
coarse resolution.
Some details of the particle tracking implementation
The following provides explanation as to why a better statistical convergence is obtained for the
mean residence times calculated with ๐ฆ 1 < ๐ฆ 2 than those calculated with ๐ฆ 1 > ๐ฆ 2 , with the same
|๐ฆ 1 โ ๐ฆ 2 |. To calculate ๐1 at each (๐ฆ 1 , ๐ฆ 2 ) combination, the same number of particles are released
88
in the fluid domain of the bed-parallel plane at ๐ฆ 1 . As illustrated in Figure B.5, among the particles
that are released at a ๐ฆ 1 higher than ๐ฆ 2 and eventually leaving the sediment (e.g., released at ๐ฆโฒ1 ,
exiting at ๐ฆโฒโฒ2 ), only a fraction of them reach a depth of ๐ฆโฒโฒ2 , leading to fewer averaging samples and
more undulatory ๐1 contour lines in Figure 5. On the other hand, for particles released at a ๐ฆ 1
lower than ๐ฆ 2 (e.g., released at ๐ฆโฒโฒ1 , exiting at ๐ฆโฒ2 ), all of them pass ๐ฆโฒ2 to exit the sediment.
Figure B.3: Sketch to demonstrate that a particle entry elevation (๐ฆ 1 ) lower than the exit elevation
(๐ฆ 2 ) tracks a larger fraction of subsurface flow paths than in the case with ๐ฆ 1 > ๐ฆ 2 , considering the
same |๐ฆ 1 โ ๐ฆ 2 |. Here, (๐ฆโฒ1 , ๐ฆโฒโฒ2 ) tracks 2 out of the 4 paths (50%) entering from ๐ฆโฒ1 , while (๐ฆโฒโฒ1 , ๐ฆโฒ2 )
tracks 2 out of the 2 paths (100%) entering from ๐ฆโฒโฒ1 .
89
Figure B.4: Streamwise dispersive velocity intensities in the sediment bed: comparison between
results obtained from DNS data taken for simulation times of 10 ( , red) and 20 ( )
large-eddy turn-over times, for the regular- (a) and random-interface (b) cases.
Figure B.5: Transit time of the 95th percentile (๐95 ) normalized by ๐ฟ/๐ข ๐ with independent
variations of ๐ฆ 1 (entry) and ๐ฆ 2 (exit) elevations, for the regular- (a) and random-interface (b) cases.
90
BIBLIOGRAPHY
91
BIBLIOGRAPHY
Amadio, G. (2014). Particle packings and microstructure modeling of energetic materials. Uni-
versity of Illinois at Urbana-Champaign.
Aubeneau, A. F., Martin, R. L., Bolster, D., Schumer, R., Jerolmack, D., and Packman, A.
(2015). Fractal patterns in riverbed morphology produce fractal scaling of water storage times.
Geophysical Research Letters, 42(13):5309โ5315.
Aubertine, C. D. and Eaton, J. K. (2005). Turbulence development in a non-equilibrium turbulent
boundary layer with mild adverse pressure gradient. Journal of Fluid Mechanics, 532:345โ364.
Baranov, V., Lewandowski, J., and Krause, S. (2016). Bioturbation enhances the aerobic respiration
of lake sediments in warming lakes. Biology Letters, 12:20160448.
Battin, T. J., Kaplan, L. A., Findlay, S., Hopkinson, C. S., Marti, E., Packman, A. I., Newbold,
J. D., and Sabater, F. (2008). Biophysical controls on organic carbon fluxes in fluvial networks.
Nature Geoscience, 1:95โ100.
Bencala, K. E. (2000). Hyporheic zone hydrological processes. Hydrological Processes,
14(15):2797โ2798.
Benettin, P., Kirchner, J. W., Rinaldo, A., and Botter, G. (2015). Modeling chloride transport using
travel time distributions at Plynlimon, Wales. Water Resources Research, 51(5):3259โ3276.
Benettin, P., Rinaldo, A., and Botter, G. (2013). Kinematics of age mixing in advection-dispersion
models. Water Resources Research, 49(12):8539โ8551.
Boano, F., Harvey, J. W., Marion, A., Packman, A. I., Revelli, R., Ridolfi, L., and Wรถrman, A.
(2014). Hyporheic flow and transport processes: Mechanisms, models, and biogeochemical
implications. Reviews of Geophysics, 52(4):603โ679.
Bottacin-Busolin, A. and Marion, A. (2010). Combined role of advective pumping and mechanical
dispersion on time scales of bed formโinduced hyporheic exchange. Water resources research,
46(8).
Botter, G., Bertuzzo, E., and Rinaldo, A. (2011). Catchment residence and travel time distributions:
The master equation. Geophysical Research Letters, 38(11).
Breugem, W. and Boersma, B. (2005). Direct numerical simulations of turbulent flow over a
permeable wall using a direct and a continuum approach. Physics of fluids, 17(2):025103.
Breugem, W., Boersma, B., and Uittenbogaard, R. (2006). The influence of wall permeability on
92
turbulent channel flow. Journal of Fluid Mechanics, 562:35โ72.
Brinkman, H. (1949). A calculation of the viscous force exerted by a flowing fluid on a dense
swarm of particles. Flow, Turbulence and Combustion, 1(1):27.
C. Carman, P. (1937). Fluid flow through granular beds. Chemical Engineering Research & Design
- CHEM ENG RES DES, 75.
Cardenas, M. B. (2015). Hyporheic zone hydrologic science: A historical account of its emergence
and a prospectus. Water Resources Research, 51(5):3601โ3616.
Cardenas, M. B. and Wilson, J. L. (2007). Dunes, turbulent eddies, and interfacial exchange with
permeable sediments. Water Resources Research, 43(8).
Cardenas, M. B., Wilson, J. L., and Haggerty, R. (2008). Residence time of bedform-driven
hyporheic exchange. Advances in Water Resources, 31(10):1382โ1386.
Chandler, I., Guymer, I., Pearson, J., and Van Egmond, R. (2016). Vertical variation of mixing
within porous sediment beds below turbulent flows. Water resources research, 52(5):3493โ3509.
Chen, X., Cardenas, M. B., and Chen, L. (2015). Three-dimensional versus two-dimensional bed
form-induced hyporheic exchange. Water Resources Research, 51(4):2923โ2936.
Christensen, K. T. and Wu, Y. (2005). Characteristics of vortex organization in the outer layer of
wall turbulence. In TSFP DIGITAL LIBRARY ONLINE. Begel House Inc.
Clauser, F. H. (1956). The turbulent boundary layer. In Advances in applied mechanics, volume 4,
pages 1โ51. Elsevier.
Cossu, C. and Hwang, Y. (2017). Self-sustaining processes at all scales in wall-bounded turbulent
shear flows. Philosophical Transactions of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 375(2089):20160088.
Dairain, A., Maire, O., Meynard, G., Richard, A., Rodolfo-Damiano, T., and Orvain, F. (2020).
Sediment stability: can we disentangle the effect of bioturbating species on sediment erodibility
from their impact on sediment roughness? Marine Environmental Research, 162:105147.
Danesh-Yazdi, M., Foufoula-Georgiou, E., Karwan, D. L., and Botter, G. (2016). Inferring changes
in water cycle dynamics of intensively managed landscapes via the theory of time-variant travel
time distributions. Water Resources Research, 52(10):7593โ7614.
Dey, S. and Das, R. (2012). Gravel-bed hydrodynamics: double-averaging approach. Journal of
Hydraulic Engineering, 138(8):707โ725.
Dutta, R., Nicolle, J., Giroux, A., and Piomelli, U. (2016). Evaluation of turbulence models
93
on roughened turbine blades. In IOP Conference Series: Earth and Environmental Science,
volume 49, page 062007. IOP Publishing.
Edwards, D., Shapiro, M., Bar-Yoseph, P., and Shapira, M. (1990). The influence of reynolds
number upon the apparent permeability of spatially periodic arrays of cylinders. Physics of
Fluids A: Fluid Dynamics, 2(1):45โ55.
Elliott, A. (1991). Transfer of solutes into and out of streambeds. PhD thesis, California Institute
of Technology.
Elliott, A. and Brooks, N. (1997a). Transfer of nonsorbing solutes to a streambed with bed forms:
Laboratory experiments. Water Resources Research, 33(1):137โ151.
Elliott, A. and Brooks, N. (1997b). Transfer of nonsorbing solutes to a streambed with bed forms:
Theory. Water Resources Research, 33(1):123โ136.
Fang, H., Han, X., He, G., and Dey, S. (2018). Influence of permeable beds on hydraulically
macro-rough flow. Journal of Fluid Mechanics, 847:552โ590.
Flack, K. A. and Schultz, M. P. (2014). Roughness effects on wall-bounded turbulent flows. Phys.
Fluids, 26:101305.
Ghodke, C. D. and Apte, S. V. (2016). Dns study of particle-bedโturbulence interactions in an
oscillatory wall-bounded flow. Journal of Fluid Mechanics, 792:232โ251.
Goharzadeh, A., Khalili, A., and Jรธrgensen, B. B. (2005). Transition layer thickness at a fluid-porous
interface. Physics of Fluids, 17(5):057102.
Grant, S. B., Gomez-Velez, J. D., Ghisalberti, M., Guymer, I., Boano, F., Roche, K., and Harvey,
J. (2020a). A one-dimensional model for turbulent mixing in the benthic biolayer of stream and
coastal sediments. Water Resources Research, 56(12):e2019WR026822.
Grant, S. B., Monofy, A., Boano, F., Gomez-Velez, J. D., Guymer, I., Harvey, J., and Ghisalberti,
M. (2020b). Unifying advective and diffusive descriptions of bedform pumping in the benthic
biolayer of streams. Water Resources Research, 56(11):e2020WR027967.
Han, X., Fang, H., He, G., and Reible, D. (2018). Effects of roughness and permeability on solute
transfer at the sediment water interface. Water research, 129:39โ50.
Harman, C. J. (2015). Time-variable transit time distributions and transport: Theory and application
to storage-dependent transport of chloride in a watershed. Water Resources Research, 51(1):1โ30.
Harvey, J. W., Bรถhlke, J. K., Voytek, M. A., Scott, D., and Tobias, C. R. (2013). Hyporheic
zone denitrification: Controls on effective reaction depth and contribution to whole-stream mass
balance. Water Resources Research, 49(10):6298โ6316.
94
Hassanizadeh, S. M. and Gray, W. G. (1987). High velocity flow in porous media. Transport in
porous media, 2(6):521โ531.
He, G., Han, X., Fang, H., Reible, D., and Huang, L. (2019). Effects of roughness reynolds number
on scalar transfer mechanisms at the sediment-water interface. Water Resources Research,
55(8):6811โ6824.
Huang, H., Ayoub, J. A., et al. (2006). Applicability of the forchheimer equation for non-darcy flow
in porous media. In SPE Annual Technical Conference and Exhibition. Society of Petroleum
Engineers.
Javadpour, F. et al. (2009). Nanopores and apparent permeability of gas flow in mudrocks (shales
and siltstone). Journal of Canadian Petroleum Technology, 48(08):16โ21.
Jimรฉnez, J. (2004). Turbulent flows over rough walls. Annu. Rev. Fluid Mech., 36:173โ196.
Kalbus, E., Schmidt, C., Molson, J., Reinstorf, F., and Schirmer, M. (2009). Influence of aquifer
and streambed heterogeneity on the distribution of groundwater discharge. Hydrology and Earth
System Sciences, 13(1):69โ77.
Keating, A., Piomelli, U., Bremhorst, K., and Nesic, S. (2004). Large-eddy simulation of heat
transfer downstream of a backward-facing step. Journal of Turbulence, 5(20):1โ3.
Khirevich, S. (2011). High-performance computing of flow, diffusion, and hydrodynamic dispersion
in random sphere packings.
Kim, J. S. and Kang, P. K. (2020). Anomalous transport through free-flow-porous media interface:
Pore-scale simulation and predictive modeling. Advances in Water Resources, 135:103467.
Kim, T., Blois, G., Best, J., and Christensen, K. T. (2017). Experimental study of turbulent structure
over permeable walls with a refractive-index-matching technique. In TSFP DIGITAL LIBRARY
ONLINE. Begel House Inc.
Kim, T., Blois, G., Best, J. L., and Christensen, K. T. (2018). Experimental study of turbulent flow
over and within cubically packed walls of spheres: Effects of topography, permeability and wall
thickness. International Journal of Heat and Fluid Flow, 73:16โ29.
Kim, T., Blois, G., Best, J. L., and Christensen, K. T. (2020). Experimental evidence of amplitude
modulation in permeable-wall turbulence. Journal of Fluid Mechanics, 887.
Kinzelbach, W. (1988). The random walk method in pollutant transport simulation. In Groundwater
flow and quality modelling, pages 227โ245. Springer.
Kirca, V. S. O., Saghebian, S. M., Roushangar, K., and Yagci, O. (2020). Influence of surface
roughness of dune bedforms on flow and turbulence characteristics. International Journal of
95
Sediment Research, 35(6):666โ678.
Knapp, J. L. A., Gonzalez-Pinzon, R., Drummond, J. D., Larsen, L. G., Cirpka, O. A., and Harvey,
J. W. (2017). Tracer-based characterization of hyporheic exchange and benthic biolayers in
streams. Water Resour. Res., 53:1575โ1594.
Kong, F. and Schetz, J. (1982). Turbulent boundary layer over porous surfaces with different
surfacegeometries. In 20th Aerospace Sciences Meeting, page 30.
Kozeny, J. (1927). On capillary movement of water in soil. German, Sitzber Akad. Wiss. Wien,
136.
Krogstad, P.-ร
. and Antonia, R. (1994). Structure of turbulent boundary layers on smooth and
rough walls. Journal of Fluid Mechanics, 277:1โ21.
Krogstad, P.-ร
. and Efros, V. (2012). About turbulence statistics in the outer part of a boundary
layer developing over two-dimensional surface roughness. Physics of Fluids, 24(7):075112.
Kuwata, Y. and Suga, K. (2016a). Lattice boltzmann direct numerical simulation of interface
turbulence over porous and rough walls. International Journal of Heat and Fluid Flow, 61:145โ
157.
Kuwata, Y. and Suga, K. (2016b). Transport mechanism of interface turbulence over porous and
rough walls. Flow, Turbulence and Combustion, 97(4):1071โ1093.
Laube, G., Schmidt, C., and Fleckenstein, J. H. (2018). The systematic effect of streambed
conductivity heterogeneity on hyporheic flux and residence time. Advances in Water Resources,
122:60 โ 69.
Lee, A., Aubeneau, A. F., and Cardenas, M. B. (2020). The sensitivity of hyporheic exchange to
fractal properties of riverbeds. Water Resources Research.
Liu, X. and Katz, J. (2018). Pressure rate of strain, pressure diffusion, and velocity pressure gradient
tensor measurements in a cavity flow. AIAA Journal, 56(10):3897โ3914.
Malard, F., Tockner, K., Dole Oliver, M. J., and Ward, J. (2002). A landscape perspective of
surfaceโsubsurface hydrological exchanges in river corridors. Freshwater Biology, 47(4):621โ
640.
Manes, C., Poggi, D., and Ridolfi, L. (2011). Turbulent boundary layers over permeable walls:
scaling and near-wall structure. Journal of Fluid Mechanics, 687:141โ170.
Manes, C., Pokrajac, D., McEwan, I., and Nikora, V. (2009). Turbulence structure of open
channel flows over permeable and impermeable beds: A comparative study. Physics of Fluids,
21(12):125109.
96
Mangavelli, S. C., Yuan, J., and Brereton, G. J. (2021). Effects of surface roughness topography in
transient channel flows. J. Turbul., 22:434โ460.
Mignot, E., Barthรฉlemy, E., and Hurther, D. (2009). Double-averaging analysis and local flow
characterization of near-bed turbulence in gravel-bed channel flows. Journal of Fluid Mechanics,
618:279โ303.
Mittal, R., Dong, H., Bozkurttas, M., Najjar, F., Vargas, A., and Von Loebbecke, A. (2008). A
versatile sharp interface immersed boundary method for incompressible flows with complex
boundaries. Journal of computational physics, 227(10):4825โ4852.
Nikora, V., Goring, D., McEwan, I., and Griffiths, G. (2001). Spatially averaged open-channel flow
over rough bed. Journal of Hydraulic Engineering, 127(2):123โ133.
Nikora, V. I., Goring, D. G., and Biggs, B. J. (1998). On gravel-bed roughness characterization.
Water Resources Research, 34(3):517โ527.
Packman, A. I., Brooks, N. H., and Morgan, J. J. (2000). A physicochemical model for colloid
exchange between a stream and a sand streambed with bed forms. Water Resources Research,
36(8):2351โ2361.
Packman, A. I., Salehin, M., and Zaramella, M. (2004). Hyporheic exchange with gravel beds:
basic hydrodynamic interactions and bedform-induced advective flows. Journal of Hydraulic
Engineering, 130(7):647โ656.
Padhi, E., Penna, N., Dey, S., and Gaudio, R. (2018a). Hydrodynamics of water-worked and
screeded gravel beds: A comparative study. Physics of Fluids, 30(8):085105.
Padhi, E., Penna, N., Dey, S., and Gaudio, R. (2018b). Spatially averaged dissipation rate in flows
over water-worked and screeded gravel beds. Physics of Fluids, 30(12):125106.
Persson, N. J., Albohr, O., Tartaglino, U., Volokitin, A. I., and Tosatti, E. (2005). On the nature of
surface roughness with application to contact mechanics, sealing, rubber friction and adhesion.
J. Phys.: Condens. Matter, 17:R1โR62.
Plimpton, S. (1995). Fast parallel algorithms for short-range molecular dynamics. Journal of
computational physics, 117(1):1โ19.
Pokrajac, D. and Manes, C. (2009). Velocity measurements of a free-surface turbulent flow
penetrating a porous medium composed of uniform-size spheres. Transport in porous media,
78(3):367.
Raupach, M., Antonia, R., and Rajagopalan, S. (1991). Rough-wall turbulent boundary layers.
Applied mechanics reviews, 44(1):1โ25.
97
Raupach, M. R. and Shaw, R. (1982). Averaging procedures for flow within vegetation canopies.
Boundary-layer meteorology, 22(1):79โ90.
Roche, K., Blois, G., Best, J., Christensen, K., Aubeneau, A., and Packman, A. (2018). Turbulence
links momentum and solute exchange in coarse-grained streambeds. Water Resources Research.
Salehin, M., Packman, A. I., and Paradis, M. (2004). Hyporheic exchange with heterogeneous
streambeds: Laboratory experiments and modeling. Water Resources Research, 40(11).
Savant, A. S., Reible, D. O., and Thibodeaux, L. (1987). Convective transport within stable river
sediments. Water Resour. Res., 23:1763โ1768.
Sawyer, A. H. and Cardenas, M. B. (2009). Hyporheic flow and residence time distributions in
heterogeneous cross-bedded sediment. Water Resources Research, 45(8).
Scotti, A. (2006). Direct numerical simulation of turbulent channel flows with boundary roughened
with virtual sandpaper. Physics of Fluids, 18(3):031701.
Shen, G., Yuan, J., and Phanikumar, M. S. (2020). Direct numerical simulations of turbulence and
hyporheic mixing near sedimentโwater interfaces. Journal of Fluid Mechanics, 892:A20.
Song, X., Chen, X., Zachara, J. M., Gomez-Velez, J. D., Shuai, P., Ren, H., and Hammond, G. E.
(2020). River dynamics control transit time distributions and biogeochemical reactions in a
dam-regulated river corridor. Water Resources Research, 56(9):e2019WR026470.
Stoesser, T., Frohlich, J., and Rodi, W. (2007). Turbulent open-channel flow over a permeable bed.
In Proceedings of the Congress-International Association for Hydraulic Research, volume 32,
page 189.
Stonedahl, S. H., Harvey, J. W., Wรถrman, A., Salehin, M., and Packman, A. I. (2010). A multiscale
model for integrating hyporheic exchange from ripples to meanders. Water Resources Research,
46(12).
Stubbs, A., Stoesser, T., and Bockelmann-Evans, B. (2018). Developing an approximation of a
natural, rough gravel riverbed both physically and numerically. Geosciences, 8(12):449.
Suga, K., Matsumura, Y., Ashitaka, Y., Tominaga, S., and Kaneda, M. (2010). Effects of wall
permeability on turbulence. International Journal of Heat and Fluid Flow, 31(6):974โ984.
Suga, K., Mori, M., and Kaneda, M. (2011). Vortex structure of turbulence over permeable walls.
International Journal of Heat and Fluid Flow, 32(3):586โ595.
Sun, Y., Park, C.-H., Pichot, G., and Taron, J. (2015). Random walk particle tracking. In
Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and
Benchmarking, pages 153โ184. Springer.
98
Taneda, S. (1956). Experimental investigation of the wake behind a sphere at low reynolds numbers.
Journal of the Physical Society of Japan, 11(10):1104โ1108.
Vanoni, V. A. (2006). Sedimentation engineering. American Society of Civil Engineers.
Venditti, J. (2013). 9.10 Bedforms in Sand-Bedded Rivers, pages 137โ162.
Voermans, J., Ghisalberti, M., and Ivey, G. (2017). The variation of flow and turbulence across the
sedimentโwater interface. Journal of Fluid Mechanics, 824:413โ437.
Voermans, J. J., Ghisalberti, M., and Ivey, G. N. (2018). A model for mass transport across the
sediment-water interface. Water Resources Research, 54(4):2799โ2812.
Volino, R. J., Schultz, M. P., and Flack, K. A. (2011). Turbulence structure in boundary layers over
periodic two-and three-dimensional roughness. Journal of Fluid Mechanics, 676:172โ190.
Whitaker, S. (1996). The forchheimer equation: a theoretical development. Transport in Porous
media, 25(1):27โ61.
Wondzell, S. M. (2006). Effect of morphology and discharge on hyporheic exchange flows in
two small streams in the Cascade Mountains of Oregon, USA. Hydrological Processes: An
International Journal, 20(2):267โ287.
Wรถrman, A., Packman, A. I., Johansson, H., and Jonsson, K. (2002). Effect of flow-induced
exchange in hyporheic zones on longitudinal transport of solutes in streams and rivers. Water
Resources Research, 38(1):2โ1.
Wรถrman, A., Packman, A. I., Marklund, L., Harvey, J. W., and Stone, S. H. (2006). Exact
three-dimensional spectral solution to surface-groundwater interactions with arbitrary surface
topography. Geophysical research letters, 33(7).
Yuan, J. and Jouybari, M. A. (2018). Topographical effects of roughness on turbulence statistics in
roughness sublayer. Physical Review Fluids, 3(11):114603.
Yuan, J., Nicolle, J., Piomelli, U., and Giroux, A. (2014). Modelling roughness and acceleration
effects with application to the flow in a hydraulic turbine. In IOP Conference Series: Earth and
Environmental Science, volume 22, page 022016. IOP Publishing.
Yuan, J. and Piomelli, U. (2014a). Estimation and prediction of the roughness function on realistic
surfaces. Journal of Turbulence, 15(6):350โ365.
Yuan, J. and Piomelli, U. (2014b). Numerical simulations of sink-flow boundary layers over rough
surfaces. Physics of Fluids, 26(1):015113.
Yuan, J. and Piomelli, U. (2014c). Roughness effects on the reynolds stress budgets in near-wall
99
turbulence. Journal of Fluid Mechanics, 760.
Yuan, J. and Piomelli, U. (2015). Numerical simulation of a spatially developing accelerating
boundary layer over roughness. Journal of Fluid Mechanics, 780:192โ214.
Zagni, A. F. and Smith, K. V. (1976). Channel flow over permeable beds of graded spheres. Journal
of the hydraulics division, 102(2):207โ222.
Zippe, H. J. and Graf, W. H. (1983). Turbulent boundary-layer flow over permeable and non-
permeable rough surfaces. Journal of Hydraulic research, 21(1):51โ65.
100