EXPLORING MICROBIAL AND HYDRAULIC CONTROLS ON REACTIVE SOLUTE
TRANSPORT AND REACTION MICROZONE FORMATION IN SURFACE-GROUND

WATER INTERFACES

By

Sinchan Roy Chowdhury

A THESIS

Michigan State University

in partial fulﬁllment of the requirements

Submitted to

for the degree of

Geological Sciences – Master of Science

2018

ABSTRACT

EXPLORING MICROBIAL AND HYDRAULIC CONTROLS ON REACTIVE SOLUTE
TRANSPORT AND REACTION MICROZONE FORMATION IN SURFACE-GROUND

WATER INTERFACES

By

Sinchan Roy Chowdhury

, Org C = 20 mgl−1

Several experimental studies have detected the presence of microzones in hyporheic sed-
iments. These microzones are small scale anoxic pores, embedded within oxygen rich porous
media and can act as anaerobic reaction sites producing reduction compounds such as nitrous
oxide, which is a greenhouse gas. Microbes are one of the key controls on nutrient transforma-
tion in hyporheic sediment and microbial biomass growth is also capable of altering hydraulic
ﬂux in sediment, causing ‘bioclogging’.
In this thesis, I developed one of the ﬁrst computa-
tional modeling approaches that combined hydraulics and microbial kinetics to explore the pres-
ence of microzones in stream sediments. The model was used to explore stream and sediment
conditions with diﬀerent hydraulic ﬂux (0.1-1.0 mday−1 Darcy ﬂux), nutrient concentrations
(O2 = 8mgl−1
, NH3 = 1− 0.5 mgl−1), and bioclog-
ging scenarios (with and without). The model domain was constructed from a pore network model
with random sized pore throat radii resulting in a heterogeneous and anisotropic ﬂow domain that
resembled a streambed, comprising of medium sand with a hydraulic conductivity of 0.8mday−1.
Results indicate that microzone formation is controlled by the hydraulic ﬂux, the nutrient concen-
tration and bioclogging strongly inhibiting microzone formation. Bioclogging scenarios typically
produced unstable microzones, which perished a few days after formation. Overall, results from
the modeling show that anoxic microzones are likely to form under many hyporheic zone condi-
tions, but their distribution and biogeochemical function will be dynamic and diﬃcult to measure
in the ﬁeld. Future investigations, will need to develop ﬁeld investigations capable of observing
microzones to aid in further model development and validation.

, NO−

3 = 1.5− 3 mgl−1

ACKNOWLEDGEMENTS

I would like to thank all memebers of my comittee: Jay Zarnetske, Phani Kumar Mantha,
Wei Zhang and Matthew Schrenk for their invaluable guidance during the course of dvelopment of
this thesis. I also thank the National Science Foundation (NSF), Environmental Science and Policy
Program (ESPP) for providing funding for the research work. I would also like to thank Joe Lee
Cullin, Ariel Shogran and Tyler Hampton for their help with feedback for my research work. Last
but not the least, I would like to thank my parents for providing emotional support and being there
for me during testing times.

iii

TABLE OF CONTENTS

v

.
.

.
.

.
.

vi

.
.
.

2.1 Model Geometry .
2.2 Flow physics .
.

LIST OF TABLES .
.
.
LIST OF FIGURES .
.
.
KEY TO ABBREVIATIONS .
CHAPTER 1

.
.
.
.
.
.
INTRODUCTION .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
CHAPTER 2 MODEL DEVELOPMENT AND CODE IMPLEMENTATION . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
.
7
.
.
.
9
2.2.1 Mass Balance .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Momentum Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.3
Fortran coding for the ﬂow physics . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Biogeochemical Kinetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Implemention of microbial kinetics in Fortran . . . . . . . . . . . . . . . . 20
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.1 Transport equation for a solute . . . . . . . . . . . . . . . . . . . . . . . . 22
Implementing solute transport in Fortran . . . . . . . . . . . . . . . . . . . 23
2.4.2

2.3.1 Monod Kinetics .
2.3.2

2.4 Solute Transport .

. .

.
.
.

.

.

.

.

.
.

.

.

CHAPTER 3 COMPUTATIONAL MODELING AND ASSESSMENT OF THE FOR-

.

.
.

. .

Introduction .

3.1
.
.
3.2 Model Formulation .
.
3.3 Numerical Simulations
.
3.4 Results and Discussions .

MATION CRITERIA FOR HYPORHEIC ANOXIC MICROZONES . . . . 24
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
.
3.4.1 Bioclogging eﬀects are signiﬁcant control on microzones . . . . . . . . . . 28
Temporal conditions with and without bioclogging . . . . . . . . 28
Spatial conditions with and without bioclogging . . . . . . . . . 31
Stable Microzones are modulated by transport(hydraulic ﬂux)
and reaction (nutrients) limitations . . . . . . . . . . . . . . . . . 34
3.4.2 Transferable Findings and Implications for Future Research . . . . . . . . . 36
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.5 Conclusion . .

3.4.1.1
3.4.1.2
3.4.1.3

. .

.

.

.

BIBLIOGRAPHY .

.

.

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

iv

LIST OF TABLES

Table 2.1: Variables in Figure 2.8 .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Table 2.2: Variables in Figure 2.9 .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Table 2.3: Variables in Figure 2.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Table 3.1: Parameter values used for the microbial kinetics modeling . . . . . . . . . . . . 27

Table 3.2: Solute and biomass parameter values for simulations . . . . . . . . . . . . . . . 27

v

LIST OF FIGURES

Figure 1.1: Schematic of the Nitrogen Cycle (Bothe et al 2007)

. . . . . . . . . . . . . . .

Figure 1.2: Nitrous Oxide production pathways(Zhu et al 2013)

. . . . . . . . . . . . . . .

Figure 1.3: Hydrologic exchange denoted by arrows in a section of stream and hyporheic

zone(Stonedahl et al 2010)

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.4: Field data from 15NO−

enrichment in oxic zones. (shown as blue in the ﬁgure from Briggs et al., 2015)

3 tracer experiments clearly shows anomolous 15N2

Figure 2.1: Model Flowchart .

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.2: Schematic of the 2D pore network . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.3: Rayleigh Distributions used for generating the porous domain . . . . . . . . . .

Figure 2.4: MATLAB Code for generating the Pore Throat Distribution . . . . . . . . . . .

1

2

2

4

6

8

8

9

Figure 2.5: Constriction of ﬂow through a pore throat. R is the unclogged radius while r

represents the the shrinking radius due to bioclogging . . . . . . . . . . . . . . 12

Figure 2.6: Bioﬁlm growing in a pore throat

. . . . . . . . . . . . . . . . . . . . . . . . . 12

Figure 2.7: Convention for mapping variables in source code: Hydraulic conductivity . . . 14

Figure 2.8: Code for calculating porosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Figure 2.9: Code for calculating the pore throat conductivities . . . . . . . . . . . . . . . . 16

Figure 2.10: Code for formulating the coeﬃcient matrix for the linear equations . . . . . . . 17

Figure 2.11: Code for calculating microbial growth and nutrient uptake rates, from Equa-

tion 2.21 to Equation 2.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 20

Figure 2.12: Code for calculating O2 concentration in a pore throat

. . . . . . . . . . . . . . 23

Figure 3.1: Hydraulic conductivity change over a 30d period. Each curve represents an

initial Darcy ﬂux value (0.1,0.5,1.0 md−)

. . . . . . . . . . . . . . . . . . . . 28

vi

Figure 3.2: a. Comparison of anoxic microzone population numbers under diﬀerent
nutrient and ﬂux conditions, but with no biomass growth occurring. b.
Comparison of microzone population numbers under diﬀerent nutrient and
ﬂux conditions when biomass growth is occurring. . . . . . . . . . . . . . . . . 29

Figure 3.3: Example spatial map of O2 concentrations for simulation of 0.5mday−1 ﬂux
with high nutrient concentrations at 400 h, noting where the demarcation
of bulk oxic and anoxic conditions is in this example, plus six examples of
anoxic (< 2mgl−1) microzones are circled in red to highlight their location.
Figure 3.4: a. Spatial oxygen proﬁles at 180 h, and b. for 240h for same hydraulic ﬂux
and nutrient concentrations (0.5md-1 ﬂux and high nutrient concentrations).

. . 30

. . 31

Figure 3.5: Temporal evolution of system Peclet number for each of the bioclogging sim-
ulations. Each row is labeled with the corresponding hydraulic ﬂux condition,
while the ﬁrst column corresponds to high nutrient ﬂux and the second col-
umn refers to low nutrient ﬂux conditions. The values are log10 transformed
and the y=0 line corresponds to a Peclet number of 1.

. . . . . . . . . . . . . . 32

Figure 3.6: Biomass proﬁle for 1.0md−1 ﬂux with high nutrient concentrations at 160
h. with a. biomass concentration in pore throats aligned along the direction
ﬂow, and b. biomass concentration in pore throats aligned transverse to the
direction of ﬂow .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

.

.

.

Figure 3.7: a Oxygen proﬁle at 200 h for 0.1md−1 ﬂux and high nutrient concentrations

and b. Oxygen proﬁle at 400 h for 0.1md−1 ﬂux and low nutrient concnetrations. 34

Figure 3.8: a Oxygen proﬁle at 400 h for 0.5md−1 ﬂux and high nutrient concentrations

and b. Oxygen proﬁle at 400 h for 0.5md−1 ﬂux and low nutrient concnetrations. 35

Figure 3.9: a Oxygen proﬁle at 400 h for 1.0md−1 ﬂux and high nutrient concentrations

and b. Oxygen proﬁle at 700 h for 1.0md−1 ﬂux and low nutrient concnetrations. 36

vii

KEY TO ABBREVIATIONS

O2

N2

Oxygen

Nitrogen

Nitrous Oxide

Nitrate

Day

Carbon

Ammonia

N2O
N O−
3
N H3
d
C
mgl−1 Milligrams per litre
Damkohler Number
Da
Pe
h,hr,hrs Hours

Peclet Number

viii

CHAPTER 1

INTRODUCTION

The element nitrogen is a key part of all ecosystems and hence undergoes multiple inorganic
and organic N forms and transformation pathways in the (Figure 1.1) [6]. However, contamination
of groundwater due to nitrogen compounds such as excess NO−
3 has become a serious worldwide
issue, and it is mostly due to anthropogenic loading caused by agricultural activities as well as
human waste [28][31]. For example, this excess of nitrogen in the ecosystem is the leading cause
for lake eutrophication and algal blooms which pose serious environmental hazards [9].

Figure 1.1: Schematic of the Nitrogen Cycle (Bothe et al 2007)

Recent research eﬀorts, however, have turned their focus on another member of the nitrogen
cycle, N2O, and the various pathways that lead up to its formation(Figure 1.2). This gas is a
greenhouse gas, that is 300 times more potent than CO2 [43][19][1][11]. The nitrogen cycling in
nature is a complex process which not only is controlled by biogeochemistry, but is also profoundly
impacted by the surrounding abiotic reactions mitigated by hydrology whether be it in lakes [32]
or streams [26]. The hyporheic zone in watersheds is a known hotspot for the transformation of
multiple nitrogen species, including the processes that create N2O (Figure 1.2). The hyporheic

1

Figure 1.2: Nitrous Oxide production pathways(Zhu et al 2013)

zone is deﬁned spatially as the subsurface domain along the wetted perimeter of a stream, surface
water-groundwater exchange occurs[39](Figure 1.3).

Figure 1.3: Hydrologic exchange denoted by arrows in a section of stream and hyporheic
zone(Stonedahl et al 2010)

Typically, this physically-driven surface water exchange brings oxic water and dissolved
solutes into the surrounding subsurface, and they interact with anoxic groundwater conditions
before reemerging downgradient back into the surface ﬂow. This kind of hyporheic exchange is

2

distinct from typical groundwater due to its bi-directional nature and that it occurs over relatively
small spatial (1mm-1km) and temporal (1s-1y) scales [5] compared to purely subsurface ﬂow. This
region is known as a biological hotspot with microbially mediated nutrient transformations, which
can have a major impact on stream ecology and water quality[22]. The hyporheic exchange forces
a ﬂux of oxygen, carbon and nitrogen compounds into and out of the stream sediment where they
are subjected to microbial transformations. The general understanding of this phenomenon is as
follows:

1. Highly oxic water ﬂows into the stream sediments as a result of down welling or pressure

gradient created due to streambed morphology.

2. Energetically favorable aerobic respiration occurs within the sediment.

3. The oxygen content gradually diminishes along the ﬂow pathway and the water turns anoxic.

This would generally be the case for homogeneous sediments[29]. However natural sediments are
rarely homogenous. Studies looking at denitriﬁcation(Figure 1.4), an anaerobic respiration pathway,
in natural heterogeneous hyporheic sediments have detected the presence of reduced nitrogen prod-
ucts in bulk oxic sediments[41].When it comes to hyporheic modeling studies regarding this issue,
simple ﬂow and reactive transport models cannot seem to explain in this phenomenon. Several so-
phisticated modeling attempts have been made to capture the eﬀects of nitrogen dynamics[13][21].

A popular method employed in understanding of reactive solute transport is the residence
time approach[8]. In this method a parcel of water is visualized as moving through a porous domain
and the time it spends moving is measured as a function of the length and ﬂow velocity. And then
the residence time is used to calculate nutrient transformation rates between any two locations in
the domain by measuring the diﬀerence in residence time and nutrient concentration between those
two locations[18].

3

Figure 1.4: Field data from 15NO−
in oxic zones. (shown as blue in the ﬁgure from Briggs et al., 2015)

3 tracer experiments clearly shows anomolous 15N2 enrichment

When it comes to any nitrogen nutrient cycling in nature, microbes are often at the heart of
the transformation process. They can be classiﬁed as aerobes, anaerobes or facultative. Strains
such as Nitrosomonas europaea play a crucial role in N2O production as it is capable of oxidizing
ammonia both under aerobic and anaerobic conditions[19].

Microbes can also impact the sediment hydraulics itself through the phenomenon of bio-
clogging. Microbial biomass growth in soil systems is not a simple growth phenomenon as
with microbes growing in colonies but is accompanied by generation of extracellular polymeric
substances (EPS). This EPS is composed mainly of polysaccharides and proteins, and it allows
microbes to better latch on to soil particle surfaces as well as acts as a buﬀer against dehydrat-
ing circumstances [38].The actual microbial biomass in the bioﬁlm consists of various species
including bacteria fungi, yeast, protozoa etc. The EPS plays a major role in governing the hydro-
dynamics of the system by physical constricting the pore space and hence reducing the overall ﬂow
through the system with time[36]. Consequently, microbial activity is also directly impacting the
hydraulics of the system and this aspect is rarely considered in hyporheic investigations and models.

The hyporheic zone is where the microbes, the hydraulics and the porous domain heterogene-
ity come together in streams and river corridors, to form an eﬃcient bio-reactor system for nutrient

4

processing[24]. A popular approach of modeling this transport phenomenon in the hyporheic zone
is the dual-domain theory where you divide up your porous domain into mobile and less-mobile
porosity to simulate a heterogeneous environment[14]. Recent work by Briggs et al 2013[7] out-
lines a method where by analyzing breakthrough curves, one can quantify the relative proportions
of these two kinds of porosity. The interaction in terms of mass transfer between the less-mobile
and mobile porosity leads to the formation of small scale anoxic regions or ‘microzones.’ So, these
microzones are relatively isolated reactive sites, embedded in a porous domain which is primarily
composed of mobile pore spaces.

In this thesis, I develop a new modeling appraoach that can quantify the impact of dynamic
bioclogging on the formation of these reaction microzones as well as the generation of N2O in
the hyporheic zone. When appropriate, to help guide the research, speciﬁc objectives and testable
hypotheses are presented in the subsequent chapters.

5

CHAPTER 2

MODEL DEVELOPMENT AND CODE IMPLEMENTATION

In this chapter, I deﬁne the geometry and formulation strategies used to construct the pore
network grid, plus the physics and chemistry implemented to model the ﬂow and transport of solute
through the network as well as the bioclogging in the system. A conceptual overview of the entire
model framework is presented in Figure 2.1. The framework shows the sequence of calculations in
the modeling process. Calculations start with the hydraulic head distribution, which is followed by
the biogeochemical kinetics consisting of calculating the nutrient uptake rate and the bioclogging
eﬀects due to biomass growth. Following that, the uptake rates and the calculated hydraulics are
used to compute the nutrient transport, the results of which are used to calculate the biogeochemical
parameters all the while adjusting the head distribution caused due to bioclogging. Then snippets
from the actual source code are presented as part of model description to demostrate how diﬀerent
aspects of the model are computed using a combination of Fortran and MATLAB.

Figure 2.1: Model Flowchart

6

2.1 Model Geometry

The modeling framework uses a 2-dimensional pore network model, comprising of spherical
pore bodies and cylindrical pore throats, with a coordination number of 4 to represent the porous
medium[27]. Each element of the pore network model (i.e. pore throats and pore bodies) has a
ﬁnite 3-dimensional volume, which when summed up, comprises the entire pore volume of the
domain. Porosity is a volume property, so the porosity of the domain is calculated as :

η =

Vpore
LB (cid:52) z

(2.1)

where η is the porosity, Vpore is the pore volume, L, B are the length and breadth of the domain
respectively and (cid:52)z is a small length in the direction orthogonal to the plane of Land B, which is
slightly larger than the pore body diameter.

Figure 2.2 presents a schematic of the structure of the pore network. Although the network
is depicted as consisting of similar sized pore throats, the pore throats are in fact variable in size
throughout the domain, simulating a heterogeneous porous media.

In order to achieve the heterogeneity, the pore radii were randomly drawn from two diﬀerent
Rayleigh distributions, diﬀering by 2 orders of magnitude in their distribution values. Figure 2.3a
reprsents the larger of the two distributions and forms the distribution for the bulk of the pore throats
in the network. The remaining pore throats ( 20% for the simulation domain used) are drawn from
the smaller distrubution shown in Figure 2.3b. This produces a heterogeneous porous media with
high spatial variability in ﬂux movement.

The source code for generating the porous domain is provided in Figure 2.4. In lines 1 and 2
of Figure 2.4, I generate the total grid for the pore network and in lines 16 and 17, two sets of arrays
are drawn from two Rayleigh distributions. Each arrays has two colums corresponding to one of the
dimensions of the pore network. In line 20 a random permutation of the indices of the two arrays
are made so that the substitutions made are completely random. The actual substitution happens
in line 26 as 20% of the values in the arrays with the larger distribution values are substituted by

7

Figure 2.2: Schematic of the 2D pore network

(a) Larger Distribution

(b) Smaller Distribution

Figure 2.3: Rayleigh Distributions used for generating the porous domain

values from the array with the smaller distribution values.

Then the next step in analysis of the model geometry is computation of the eﬀective pore
throat length and and the pore volume of the generated pore network. The eﬀective length needs to
be calculated as the pore throats are shortened by the pore bodies at both ends.

lj =

Length o f

the domain in direction j

Number o f pore bodies − 1

− 2 ∗ Rporebody

(2.2)

where lj is the eﬀective length of the pore throat in j direction (either along the length or breadth

8

Figure 2.4: MATLAB Code for generating the Pore Throat Distribution

of the domain) and Rporebody is the pore body radius. This represents the actual throat length
traversed by the hydraulic ﬂux while traveling from one pore body to an adjacent one. The total
pore volume of the domain is calculated as a summation of the total volume of the pore throats and
pore bodies in the system.
Vpore = (N ∗ M)(4

3 πR3) + (N − 1) ∗ (M − 2) ∗ (πr2

(2.3)

iyly) + (M − 1) ∗ N ∗ (πr2

ixlx)

where Vpore = total pore volume, R = pore body radius, ly/lx = eﬀective pore throat length in y or
x direction, riy/rix = pore throat radius at the ith location in y or x direction
2.2 Flow physics

The ﬂow in the pore network model is assumed to be a quasi-steady ﬂow without any source
or sink present anywhere within the network. The ﬂow in each pore throat is computed using

9

a modiﬁed Hagen-Poiseuille equation which can take into account the eﬀects of bioclogging,
microbial mass grows in the system [35].

2.2.1 Mass Balance

Equations were ﬁrst derived to determine the hydraulic head distribution through out the

netwrk. The mass continuity equation for a steady state, incompressible ﬂow can be written as:

Þ
Þ

(cid:174)(cid:53) · (cid:174)v = 0
((cid:174)(cid:53) · (cid:174)v)dv = 0
((cid:174)v · (cid:174)n)da = 0


AiVi = 0

(2.4)

(2.5)

(2.6)

(2.7)

i

where v= the instantaneous velocity at a point, A = area of the control surface, V = average velocity
of ﬂow Therefore, the hydraulic head distribution can be achieved simply by solving Equation 2.7
at every single pore node in the network.

2.2.2 Momentum Balance

In this section, I am deriving the Hagen-Poiseuille equation for laminar ﬂow through a cirular
pipe of radius R. The momentum balance equation for a ﬂuid is expressed by the Navier-Stokes
equation (for a Newtonian ﬂuid):

∂(ρ(cid:174)v)
∂t

+ (cid:174)v · (cid:174)(cid:53)(ρ(cid:174)v) = −(cid:174)(cid:53)p + µ (cid:53)2 (cid:174)v +

1
3 µ(cid:174)(cid:53)((cid:174)(cid:53) · (cid:174)v) + ρ(cid:174)g

(2.8)

where ρ is the density of the ﬂuid, v is the ﬂow velocity, p is the pressure, µ is the dynamic viscosity
and g is the acceleration due to gravity.
In order to simplify Equation 2.8 certain assumptions are made:

1. The ﬂow is incompressible

10

2. The pore throat is a cylinder, so equations are solved in the cylindrical coordinate system

3. The ﬂow is in steady state

4. Gravity forces are neglected

5. Velocity only in longitudinal direction is considered

These assumptions allows Equation 2.8 to be simpliﬁed to

1
µ
Integrating equation 2.9 twice with respect to r :

(r ∂vx
∂r

) =

1
r

∂
∂r

∂p
∂x

vx =

1
4µ

∂p
∂x

r2 + C1 ln r + C2

(2.9)

(2.10)

vx has a ﬁnite value at r = 0 and since ln r is undeﬁned at r = 0, C1 = 0. vx is also 0 at

r = R. From the second boundary condition C2 = − 1
4µ
Thus, the ﬁnal Hagen-Poiseuille equation has the form:

∂p

∂x R2

vx = − 1
4µ

∂p
∂x

(R2 − r2)

(2.11)

Based on Equation 2.12 , thte toal mass ﬂux through the cross-section of a cylinder will be equal
to:

Þ R

0

(2.12)

(2.13)

q = −

(R2 − r2) 2πrdr

1
∂p
4µ
∂x
q = − π
8µ

R4

∂p
∂x

In the pore throats of the network, the bioﬁlm is assumed to grow symmetrically from the surface
of the throat towards the center (Figure 2.6) and this way over time, it clogs the pore throat thereby
reducing the hydraulic ﬂux through the section (Figure 2.5) [35][12].

Applying Equation 2.13 seperately to both the open section of the pore throat as well as the
section covered by the bioﬁlm, the total hydraulic mass ﬂux through the pore throat comes out to

11

Figure 2.5: Constriction of ﬂow through a pore throat. R is the unclogged radius while r represents
the the shrinking radius due to bioclogging

Figure 2.6: Bioﬁlm growing in a pore throat

be:

Substituting − ∂p
∂x

qtotal = − π
8µ
by ρg(cid:52)h
l

∂p
∂x

r4 + (− π
8µbio

∂p
∂x

(R4 − r4))

the ﬁnal form of the equation is obtained as:

qtotal =

(cid:52)h
l

π ρg
8µ

(r4 +

(R4 − r4)

X

)

(2.14)

(2.15)

qtotal is the total ﬂux through a bioclogged pore throat, ρ is the density of water, µ is the
viscosity of water, g is the acceleration due to gravity, l is the length of the pore the pore throat, X
is the ratio of the viscosity of the bioﬁlm to water, R is the total radius of the pore throat, r is the
radius of the unclogged section of the pore throat.

12

Referring back to Equation 2.7, the mass balance equation required to be solved at every pore

node in order to obtain the hydraulic head distribution is given by:

i=4

i=1

(cid:52)hi
li

π ρg
8µ

(r4
i +

i − r4
i )
(R4
X

) = 0

(2.16)

The mapping between r and R is a function of the volume of biomass in a pore throat. Let
the conentration of bioﬁlm (biomass per unit pore volume) be B, the density of the bioﬁlm be ρbio
and the length of the pore throat be l. Then the relation between r and R can be written as:

BπR2l = ρbioπ(R2 − r2)l

r = R(1 − B
ρbio

)0.5

(2.17)

(2.18)

13

2.2.3 Fortran coding for the ﬂow physics

In this section, I present some snippets of the fortran code used for computing the hydraulic

head distribution.

H3

Ke2(i-1,j)

H1

H

Ke1(i,j)

H2

Ke1(i,j+1)

Ke2(i,j)

H4

10

9

8

7

6

5

4

3

2

1

1

2

3

4

5

6

7

8

Figure 2.7: Convention for mapping variables in source code: Hydraulic conductivity

Some variables in the code have a speciﬁc convention with a number included either in the

name or in the index. The numbers read as follows:

• 1 refers to the x-direction

• 2 refers to the y-direction

• 3 refers to the node itself

For example, Figure 2.7 shows the conductivity variable Ke, which holds the pore throat
conductivities,referenced around a node (i, j). The way this type of referenced variables are read in

14

the the code is, in the x-direction, the value held by the variable is always to the left of the node and
in the y-direction the value referenced is below the node. So in Figure 2.7, Ke1(i, j) refers to the
conductivity of pore throat connecting H and H1 and Ke2(i, j) refers to the conductivity of pore
throat connecting H and H4.

Figure 2.8: Code for calculating porosity

Variable Deﬁnition
Vx
Vy
por
Ra
Lx
Ly
Re
n,m

Pore volume of all pore throats aligned in x-direction
Pore volume of all pore throats aligned in y-direction
Porosity of the domain
Radius of pore bodies
Length of the domain in x-direction
Length of the domain in y-drection
Radius of pore throats
Grid size parameters

Table 2.1: Variables in Figure 2.8

Figure 2.8 shows the code for calculating the porosity of the domain. The commented lines

contain OpenMP directives in order to carry out the task using parallel computing.

15

Figure 2.9: Code for calculating the pore throat conductivities

Variable Deﬁnition
lenx
leny
rho
g
meu
Ke
re11
re22

Eﬀective length of pore throats aligned in x-direction
Eﬀective length of pore throats aligned in y-direction
Density of water
Acceleration due to gravity
Dynamic viscosity of water
Conductivity of pore throats
Constricted radius of pore throats aligned in x-direction
Constricted radius of pore throats aligned in y-drection

Table 2.2: Variables in Figure 2.9

The conductivity of a pore throat is caluclated as follows:

Ke =

qtotal(cid:52)h

(2.19)

where (cid:52)h is the hydraulic head diﬀerence across a pore throat. Therefore from Equation 2.15, the
analytical relation for Ke is

Ke =

(r4 +

π ρg
8µl

(R4 − r4)

X

)

(2.20)

Figure 2.9 shows the formulation of Equation 2.20 in the code.

16

Figure 2.10: Code for formulating the coeﬃcient matrix for the linear equations

The 2 dimensional ﬂow ﬁeld results in a sparse pentadiagonal coeﬃcient matrix A for the
equation Ax = b, where x represents the vector of unknown hydraulic head values and b is the
vector containing the boundary conditions. Figure 2.10 shows the formulation of the coeﬀcient
matrix where the vectors A1, A2, A3, A4, A5 represents the diagonals of the matrix A. The solver
used for determining the vector x is the PARDISO solver available as a part of the Intel Math Kernal
Library.

17

2.3 Biogeochemical Kinetics

In this section, I discuss the methodologies implemented to calculated the biomass numbers

as well as the nutrient values used in the model.

2.3.1 Monod Kinetics

s

a

dM
dt

= (µo[

s

a

][

n

][

][

o

][

Ko + o

Kao + a

]−ko)+(µn[

3 respiration) as summarized below.

One of the most popular and robust methods of modeling microbial growth throughout the
decades has been the Monod model [23][37][42]. In this study, I implemented the multiple Monod
Kinetic model outlined in Widdowson et al 1988[37], which can track microbial biomass growth
both through aerobic and anaerobic conditions (NO−
1
Kso + s
M
M= microbial biomass
µo= growth rate for oxic conditions
µn growth rate for anoxic conditions
s, o, n, a= cocnetration for organic carbon substrate, O2, NO−
Kn, Ksn, Kan= respective substrate saturation constants under NO3 respiration
Ko, Kso, Kao= respective substrate saturation constants under O2 respiration
In= an inhibitor factor = (1 + O
Kc
ko, kn= microbial decay coeﬀcients under oxic and anoxic conditions

)−1 where Kc is a paramter

]−kn)In (2.21)

Ksn + s

Kn + n

Kan + a

3 and NH3

Equation 2.21 allows the tracking of microbial biomass under both oxic and anoxic conditions.
The inhibition factor In modulates the NO−
3 respiration in Equation 2.21 by assuming small values
when the O2 concentrations are high and gradually increases in value as the O2 gets depleted over
time. Now when it comes to other available equations that can be used to track biomass growth,
there are models which take into consideration only the nutrient which has the lowest limiting
constraint instead of all the solutes, but Equation 2.21 has been employed here as it is one of
the most frequently occuring and mechanistically accurate models in environmental engineering
literatures.

18

Once the mass is calculated, the nutrients consumed in that process is then estimated using

stoichiometric coeﬃcients.

rs =

[

µo
Yo

o

Ko + o

][

n

][

s

][

a

Kan + a

]In

s

][

Kso + s
rn = ηµn[
ro = γ µo[
s

][

] +

[

µn
Yn

Kao + a

a

n

Kn + n

o

Ko + o
a

s

s

Kn + n

Ksn + s
a
]In
]
Kao + a
s
][

][
Kan + a
][
n

][
Ksn + s
][
] +  µn[
3 , O2, NH3 respectively

Kso + s

a

(2.22)

(2.23)

(2.24)

(2.25)

][

a

Kan + a

]In

ra = ψ µo[

o

][

Kn + n

Ko + o

Ksn + s

Kso + s

Kao + a
rs, rn, ro, ra = uprake rates for substrate, NO−
Yo, Yn= yield coeﬃcients under oxic and anoxic conditions respectively
γ = O2 use coeﬃcients for heterotrophic biomass
η = NO−
, ψ = NH3 use coeﬃcients for NO−

3 use coeﬃcient for heterotrophic biomass growth
3 and O2 respectively

When computing for the biomass using Equation 2.21 the following assumptions are made:

1. Planktonic microbial populations are neglected.

2. The microbial biomass growth includes the entire mass of the bioﬁlm.

3. Biomass attachment and detachment mechanisms are not considered.

4. All biogeochemical reactions take place within only the bioﬁlm.

Moreover, the bioﬁlm in the system is considered to be permeable, allowing transport through

the bioﬁlm by means of both advection and diﬀusion.

19

2.3.2

Implemention of microbial kinetics in Fortran

The following ﬁgure shows the imlementation of the Monod Kinetics in Fortran.

Figure 2.11: Code for calculating microbial growth and nutrient uptake rates, from Equation 2.21
to Equation 2.25

20

Variable Deﬁnition
Inhibition factor
Ii
Bnold
Biomass concentration at the old time step
BnNew Biomass concentration at the new time step
dtreac
rnos
rhob
Anold
Onold
Snold
Nnold
mo
mn
Kdo
Kdn

Reaction time step
Nitrous oxide production rate
Density of bioﬁlm
Ammonia concentration at old time step
Oxygen concentration at old time step
Substrate concentration at old time step
Nitrate concentration at old time step
Growth rate under toxic condition
Growth rate under nitrate respiration
Decay coeﬃcient under oxic conditions
Decay coeﬃcient under anoxic conditions

Table 2.3: Variables in Figure 2.11

Table 2.3 lists out the variables for Figure 2.11 and any variable not listed in the table were
deﬁned prior to this. Figure 2.11 shows the computation for the microbial biomass and nutrient
uptake rates, from Equation 2.21 to Equation 2.25.

21

2.4 Solute Transport

In the ﬁnal section for this chapter, I discuss the numerical techniques implemented to bring
together the biogeochemical parameters and the advection relations developed in previous sections
in order to formulate a multispecies transport model.

2.4.1 Transport equation for a solute

Solute transport is modeled using the standard advection-diﬀusion equation [4]. The for-
mulation consists of a 1 dimensional transport equation at every pore throat and a 2 dimensional
equation at every pore body.

= D((cid:53) · ((cid:53)C)) − v · (cid:53)C + S

∂C
∂t

(2.26)

where C is the solute concentration, v is the ﬂuid velocity, D is the diﬀusion coeﬃcient, S is
a source/sink term for the solute. Integrating Equation 2.26 with respect to the control volume, the

equation becomesÞ ∂C

∂t

(D((cid:53) · ((cid:53)C)) − v · (cid:53)C − ((cid:53) · v)C + S)dV

Þ
V = DA((cid:53)C) +
DA((cid:53)C) + AvaverageC

dV =

∂C
∂t

∂C
∂t

=

AvaverageC + SV

V

+ S

where vaverage is the average ﬂux velocity across the control surface A.

(2.27)

(2.28)

(2.29)

Applying a ﬁnite diﬀerence explicit scheme to discretize the derivative terms and incorporat-
ing the modiﬁed Hagen-Poiseuille equation (Equation 2.15) to calculate the advective ﬂux across
the control surfaces:

i + (

i

Cn+1
i

= Cn

(cid:52)hi
li

π ρg
8µ

(r4
i +

i − r4
(R4
i )
X

i − Cn
Cn
i−1

li

(cid:52) t) 1
V

DAi

+ S (cid:52) t

(2.30)

i (cid:52) t +

)Cn

i

As mentioned earlier, the bioﬁlm is treated as permeable, therefore it allows solute transport not
only by advection, but also by diﬀusion. The diﬀuive ﬂux is treated as two seperate quantities with

22

(cid:52)hi
li

π ρg
8µ

(r4
i +

Cn+1
i

= Cn

i +(

i (cid:52)t+

one part ﬂowing across the open section and othe other ﬂowing across the bioclogged section. So
the ﬁnal form of the solute transport equation for a single specie is given by:
i ))Cn

+S(cid:52)t
(2.31)
n + 1 refers to the next time step while n refers to the current time step. Dw refers to the diﬀusion
coeﬃcient for solute ﬂowing through the open section and Db refers to the diﬀusion coeﬃcient for
solutes ﬂowing through the bioclogged section.

(Dwπr2

i +Dbπ(R2

i −r2

i )
i − r4
(R4
X

)Cn

i − Cn
i−1

li

(cid:52)t) 1
V

i

i

Equation 2.31 is solved at every pore throat and pore body in the network individually for
every single solute being tracked. Since the biogeochemical reactions are occuring only in the pore
throats, S is zero for every solute at the pore bodies.

2.4.2

Implementing solute transport in Fortran

This ﬁnal code section and end of modeling chapter, shows the implementation of Equation
2.31 in fortran. The code is for solute transport of O2 through a pore throat and hence it is a one
dimensional equation.

Figure 2.12: Code for calculating O2 concentration in a pore throat

Dir and aDir are multidimensional arrays of binary values 0 or 1 which corrects for the

direction of the advective solute ﬂux depending upon hydraulic head gradient.

23

CHAPTER 3

COMPUTATIONAL MODELING AND ASSESSMENT OF THE FORMATION

CRITERIA FOR HYPORHEIC ANOXIC MICROZONES

3.1 Introduction

In a valley bottom, along the river corridor, the hyporheic zone is a known control point for
biogeochemical conditions in the watershed. Hyporheic zones, are regions of bi-directional ﬂow
in the river corridor, where surface ﬂux enters the sediment and mixes with groundwater and then
exits the sediment over relatively short spatial scales[5]. This ﬂux of water also exchanges critical
nutrients like NO−
3 ,dissolved organic carbon (DOC) along with O2 in the sediment matrix, leading
to enhanced microbial activity in the hyporheic zone relative to the surface water [13][40][41].
Oxygen is thermodynamically favored for microbial respiration in the hyporheic zone, leading to
many hyporheic pore waters transitioning from bulk oxic conditions to bulk an anoxic conditions
as the O2 is respired. However, recent ﬁeld studies, including isotopic tracers of redox sensitive
solutes, indicate the presence of anaerobic denitriﬁcation products, namely N2O and N2 forming
in bulk oxic porous domains of hyporheic zones[15][41][8]. The bulk oxic conditions should
inhibit anaerobic processes, and yet anaerobic processes products exist. To explain this apparent
paradox, the theory of anaerobic ‘microzones’ was created, or pore spaces zones depleted in oxygen
embedded within oxygen-rich porous domain. These anaerobic microzones can arise from spatial
heterogeneities of physical and biogeochemical conditions in the stream sediments (e.g., variable
pore geometries, [8] variable carbon availability [29]). As indicated by the isotope tracer studies,
these microzones, if present and stable, could facilitate multiple functions, not currently accounted
for in hyporheic modeling studies, like enhanced formation of N2O, a greenhouse gas over 300
times more potent than CO2 [11][19] [43]. So, accounting for the eﬀect of microzones in any
model trying to predict the production of N2O may be crucial.

The microbial metabolic activity that likely plays a key role in the development of these

24

microzones, can also lead to the expansion of microbial biomass in the sediment matrix over
time [23][37]. This gradual accumulation of biomass, which primarily appears as the formation
of bioﬁlms, attached to sediment particles, can alter the hydraulic ﬂux through the sediment.
This biomass accumulation can lead to a phenomenon known as ‘bioclogging,’ which causes a
reduction in overall hydraulic conductivity and porosity [36] [34] [35]. Microbial biomass growth
is comprised of microbial cells growth as well as the accumulation of extra cellular polymeric
substances, which together comprise the structure of the bioﬁlm[3]. A multitude of studies have
observed the phenomenon of bioclogging and some have even showed an increase in hydraulic
dispersivity values due to bioclogging [30][33]. Despite its potential importance, bioclogging has
received little attention as a physical phenomenon in hyporheic hydraulic analysis and it should be
a part of studies and models, especially when reactive transport due to microbial activity is under
investigation [10]. However, several modeling eﬀorts for porous sediments, unrelated to hyporheic
zones, have tried to capture the mechanisms of bioclogging, either through colony formation [36]
or by means of bioﬁlm formations in pore network models [12][35].

As microzones in stream sediments is an emerging topic, modeling eﬀorts on microzones
formations in the hyporheic zone have been limited, with the most notable work being by Briggs
et al. (2015) [8] and Sawyer (2015) [29], and most of the other research being experimental and
inconclusive [16][41]. In the Briggs et al. (2015) simulation study, they did not account for the for
the role of variable stream chemistry or bioclogging, as it was predicated solely on the assumption
that the key controlling factors in determination of the microzone formation are residence time and
Damkohler number for O2 (i.e., DaO2, a number indicating if anoxia is transport or reaction rate
limited), as applied in previous hyporheic studies [42]. While the residence time was treated at the
pore scale time of travel, the DaO2 was derived with a constant O2 uptake rate for each simulation
throughout the model domain, thereby ignoring the heterogeneity in the distribution of biomass in
a porous domain and the eﬀect of availability of other nutrients on biomass growth. Consequently,
the aim of the present study is to bring together the various process-based models available for
hydraulics, nutrient transformations, and biomass development, to create more realistic mechanistic

25

simulations that can track dynamic hydraulic, microbial, and chemical conditions in a hyporheic
zone that is experiencing diﬀerent stream chemistry and microbial biomass growth conditions.
Through the development and use of this new model, our objective was to assess the interacting
controls on hyporheic pore water chemistry in space and time, with emphasis on revealing the
understudied role of bioclogging. Our initial guiding hypothesis for the eﬀect of bioclogging in
the hyporheic zone was: the conditions that promote microbial biomass growth (e.g., nutrient and
substrate availability) in the sediments would, in turn, create more heterogeneity in the ﬂowpaths
and residence times of the solutes and increase the abundance of anoxic microzones. To assess our
objective and test this bioclogging hypothesis, we explored a range of realistic stream and hyporheic
conditions that lead to or inhibited the formation and persistence of anoxic microzones occurring
in bulk oxic pore water.

3.2 Model Formulation

Please refer to Chapter 2 for details on the model formualation

3.3 Numerical Simulations

The numerical experiments explored diﬀerent hyporheic scenarios where realistic hydraulic
ﬂux, nutrient concentrations, and bioclogging conditions were varied. In total, 12 sets of simulations
were performed comprising of 3 diﬀerent hydraulic ﬂux conditions run for two diﬀerent states of
nutrient conditions (i.e., “high” and “low” concentration) and ﬁnally with and without bioclogging
occurring. The high and low concentrations were selected to represent typical more polluted “high”
and pristine “low” conditions, while reﬂecting the general stoichiometry observed in freshwaters
[20]. Each simulation utilized parallel processing with 40 Intel Xeon cores with a turbo clock
speeds of 3.7 Ghz and 27 MB cache memory each and total time for 1 simulation set, depending on
hydraulic parameters, ranged between 5 to 50 hours.The time step for numerical simulation is taken
as the value satisfying both the Courant and diﬀusion number criteria. To reduce the computational
load, the microbial kinetic step is computed every 30 seconds.

26

Parameters Values
µo
µn
ko
kn
Yo
Yn
γ
η

ψ
Kn
Ksn
Kan
Ko
Kso
Kao
Dw
Db

Source
1.3e − 4s−1
Values assumed
1.2e − 4s−1
Values assumed
0.23e − 6s−1
Values assumed
0.23e − 6s−1
Values assumed
Values assumed
0.34
Values assumed
0.38
Widdowson et al 1988
1.4
Widdowson et al 1988
2.2
0.122
Widdowson et al 1988
Widdowson et al 1988
0.122
0.0026 mgcm−1 Widdowson et al 1988
0.04 mgcm−1
Widdowson et al 1988
0.04 mgcm−1
Widdowson et al 1988
0.00077 mgcm−1 Widdowson et al 1988
0.04 mgcm−1
Widdowson et al 1988
0.001 mgcm−1
Widdowson et al 1988
2e − 9
Briggs et al 2015
2e − 10
Ezeuko et al 2011

Table 3.1: Parameter values used for the microbial kinetics modeling

Variables
Initial and inﬂux O2 cocentration
Initial and inﬂux organic C concentration
Initial biomass concentration(mass per unit pore volume)
Density of Bioﬁlm
Viscosity ratio (X)
Total simulation time
Initial and inﬂux NO−
3 concentration (HIGH)
Initial and inﬂux NH3 concentration (HIGH)
Initial and inﬂux NO−
3 concentration (LOW)
Initial and inﬂux NH3 concentration (LOW)

Values
8mgl−
20mgl−
1e − 4Kgm−
20Kgm−
1e9
30days
3mgl−
1mgl−
1.5mgl−
0.5mgl−
Table 3.2: Solute and biomass parameter values for simulations

27

3.4 Results and Discussions

3.4.1 Bioclogging eﬀects are signiﬁcant control on microzones

3.4.1.1 Temporal conditions with and without bioclogging

Numerical experiments where all hydraulic, biogeochemical, and microbial biomass were
dynamic had similar outcomes. Namely, most scenarios start displaying signiﬁcant bioclogging
eﬀects after 4 days with the simulations for higher nutrient ﬂux displaying faster clogging than
their low nutrient ﬂux counter parts (Figure 3.1). Hence, bioclogging fundamentally shifts the

Figure 3.1: Hydraulic conductivity change over a 30d period. Each curve represents an initial
Darcy ﬂux value (0.1,0.5,1.0 md−)

ability of water and solutes to move through the simulated hyporheic sediments. The pattern of
reduction in hydraulic conductivity due to bioclogging shows similar behavior as seen in the 2
dimensional simulations for pore network models performed by Thullner and Baveye, 2008 [35] for
rock type porous media. The general thresholding behavior for hydraulic conductivity (Figure3.1)

28

(a)

(b)

Figure 3.2: a. Comparison of anoxic microzone population numbers under diﬀerent nutrient and
ﬂux conditions, but with no biomass growth occurring. b. Comparison of microzone population
numbers under diﬀerent nutrient and ﬂux conditions when biomass growth is occurring.

shows that that bioclogging initially proceeds slowly, but after reaching a threshold biomass volume
level, results in a rapid decrease of the hydraulic conductivity. This decrease ceases at a very low,
stable hydraulic conductivity that is determined and bounded by the hydraulic ﬂux rate through the
bioﬁlm.

When we document and enumerate the overall anoxic microzone conditions through time,
the microzone population that occurs under a given simulation set is provided in Figure 3.2. The
percentage of anoxic pores which are formed in the bulk oxic region at any given time, where
the oxic region is deﬁned as the majority of pores (i.e > 50%) having O2 above 2mgl−1. We
demarcated the border between the oxic and anoxic portions of the porous domain by selecting the
column of the spatial grid where the number of oxic pores was less than or equal to 50% of the total
pores along the grid column. The results clearly show that anoxic microzones are formed under
most simulation conditions. However, stable microzones (Figure 3.3) (i.e., microzones that persist
throughout the total simulation period) are only generated under circumstances where bioclogging
is absent in the hyporheic zone. Figure 3.2 also shows that changing the hydraulic and nutrient
ﬂux conditions also determine the abundance and stability of microzones in the hyporheic zone.

29

Figure 3.3: Example spatial map of O2 concentrations for simulation of 0.5mday−1 ﬂux with high
nutrient concentrations at 400 h, noting where the demarcation of bulk oxic and anoxic conditions is
in this example, plus six examples of anoxic (< 2mgl−1) microzones are circled in red to highlight
their location.

Under conditions without the presence of bioclogging, ﬂuxes with higher nutrient concentrations
start developing microzones much earlier than their low nutrient ﬂux counterparts, but all the
simulations tend to move towards an microzone equilibrium condition. While under the inﬂuence
of bioclogging, the microzones do start forming but they are not sustained and disappear only hours
after forming.

Previous simulation studies (Briggs et al., 2015) linked the formation of microzones to
residence time with respect to a steady DaO2. But our results show that residence time alone
cannot predict the total number of microzones as the heterogeneity in the distribution of biomass
and its subsequent behaviors (i.e., growth, respiration) under diﬀerent nutrient ﬂux conditions,
plays a key role in the emergence of the microzones.

30

3.4.1.2 Spatial conditions with and without bioclogging

None of the bioclogging simulations were able to produce stable microzones across the
hyporheic zone. This may be due to the fact that each simulation starts oﬀ as an advection dominated
system, but due to increasing eﬀects of bioclogging, slowly transitions into a diﬀusion dominated
system. And, then, if given enough time, diﬀusion eliminates any concentration gradients in the
domain leading to the loss of the anoxic microzones. For example, as the concentration proﬁles
evolve through time (Figure 3.4), the bioclogged system approaches more or less uniform O2
concentration distribution.

(a)

(b)

Figure 3.4: a. Spatial oxygen proﬁles at 180 h, and b. for 240h for same hydraulic ﬂux and nutrient
concentrations (0.5md-1 ﬂux and high nutrient concentrations).

To further illustrate the impact of bioclogging and the prevalence of advection dominating
the outcome of microzone formation, we determined the evolution of the Peclet number (ratio of
advective and diﬀusive ﬂux) in the porous domain under the inﬂuence of bioclogging (Figure 3.5).
The resulting Peclet patterns show that every simulation starts with an advection dominated
system (i.e., Peclet number > 1; shown as > 0 in Figure 3.5 due to log transformation in plot),
but gradually over a period of days transitions into a diﬀusion dominated system. This transition
occurred at roughly the same timeframe of 6-9 d when nutrients were high, but the transition was
more gradual and took longer to achieve when nutrients were low.

31

Figure 3.5: Temporal evolution of system Peclet number for each of the bioclogging simulations.
Each row is labeled with the corresponding hydraulic ﬂux condition, while the ﬁrst column cor-
responds to high nutrient ﬂux and the second column refers to low nutrient ﬂux conditions. The
values are log10 transformed and the y=0 line corresponds to a Peclet number of 1.

In summary, the simulation results contradict our original hypothesis where we stated that

bioclogging would enhance the microzone population.

This is because we presumed that the biomass growth would take place more or less throughout
the domain, impacting the smallest pore throats preferentially. In fact, as Figure 3.6 demonstrates
clearly, the major biomass growth takes place at the entrance of the hyporheic zone, creating a
contiguous bioﬁlm layer perpendicular to and eﬀecting all ﬂow through the domain. Recent results
from an independent and more simpliﬁed modeling framework of hyporheic bioclogging [10] also
observed this focused and eﬀective bioclogging occuring directly at the sediment interface. These
consistent ﬁndings suggest that bioclogging is likely to be a localized phenomena that is occurring
near the inﬂuent (i.e., downwelling) region of hyporheic zones composed of sandy sediments as

32

(a)

(b)

Figure 3.6: Biomass proﬁle for 1.0md−1 ﬂux with high nutrient concentrations at 160 h. with a.
biomass concentration in pore throats aligned along the direction ﬂow, and b. biomass concentration
in pore throats aligned transverse to the direction of ﬂow

in our study. We are also able to show that the biomass at the distal end of the ﬂow paths does
not attain the same growth levels because they are deprived of the essential nutrients required for
their growth. Our ﬁndings are corroborated by ﬁeld observations,[25][2] where microbial biomass
was found to be accumulated at shallowest depths of streams, indication that bioclogging is indeed
a localized phenomenon, controlled by biomass accruing at the inﬁltration sites of the sediment
domain.

In the future, detailed ﬁeld and laboratory investigations should be able to validate our
simulated distribution of bioﬁlm abundance. It should be acknowledged that this study is constrained
by the fact the simulation conditions represents more of a controlled environment in a laboratory,
rather than a natural setting in the ﬁeld. In the ﬁeld, access to studying these small scale issues can
be challenging due to water depth, nondestructive sampling of sediments, and that stream beds are
dynamic. This last point of bed dynamics is worth discussing because beds composed of medium
sand, like that our simulations, are typically mobile (eroding and depositing) under most ﬂow
conditions [17]. This mobile bed results in a continuously changing streambed morphology and
hence it is likely that the bioclogged layers forming at the inﬂuent of the sediment can be eroded
and transported away resulting in the regular declogging of the system. However, given more stable

33

bed conditions, bioclogged sediments are likely to form in stream and river hyporheic zones [2].

3.4.1.3 Stable Microzones are modulated by transport(hydraulic ﬂux) and reaction (nutri-

ents) limitations

While bioclogging is critical to inhibiting microzones formation, the experiments without
bioclogging reveal other criteria needed for stable microzone formation. First, stable microzone
formation is not guaranteed for a hyporheic zone where bioclogging does not occur. For example, no
microzones were stable at 0.1md−1 Darcy ﬂux and high nutrient concentrations without bioclogging
(Figure 3.7a).The highest percentage of microzone obtained for this simulation was 16.67%, but
the microzones were only ephemeral and disappeared by 300 h. This is mainly due to the fact
that the high nutrient availability increased the oxygen respiration rates, while the hydraulic ﬂux
was not suﬃcient enough to replenish the oxygen consumed. Hence, the hyporheic domain rapidly

(a)

(b)

Figure 3.7: a Oxygen proﬁle at 200 h for 0.1md−1 ﬂux and high nutrient concentrations and b.
Oxygen proﬁle at 400 h for 0.1md−1 ﬂux and low nutrient concnetrations.

becomes bulk anoxic and no microzones can exist in such situations. However, at lower nutrient
concentrations, the oxygen uptake rate is more limited (Figure 3.7b), allowing stable microzones
(8.33% o f oxic domain) to form even under the low hydraulic ﬂux conditions.

34

The highest fraction of stable microzones formed under two types of conditions when neither
transport or reaction limitations occurred. First, for intermediate hydraulic ﬂux 0.5md−1 and high
nutrient conditions, (7.9% o f oxic domain) and second for low hydraulic ﬂux 0.1md−1 and low
nutrient conditions (8.33% o f oxic domain). However, for the 0.5md−1 ﬂux with high nutrient
concentrations (Figure 3.8a), the bulk oxic domain covered a substantially larger fraction of the
domain when compared to the low ﬂux scenario (Figure 3.7b). For low nutrient concentrations with
0.5md−1 ﬂux, (Figure 3.8b) the overall domain remains much more oxygenated and microzones
formed are much lower (3.65% o f oxic domain) due to less oxygen uptake, caused by the lack of
nutrients. Figure 3.9a shows the spatial distribution of the microzones for a ﬂux rate of 1.0md−1.

(a)

(b)

Figure 3.8: a Oxygen proﬁle at 400 h for 0.5md−1 ﬂux and high nutrient concentrations and b.
Oxygen proﬁle at 400 h for 0.5md−1 ﬂux and low nutrient concnetrations.

These results are in line with expectations based on previous simulation studies [42][8], where
the fraction of microzones in the bulk oxic domain (5.76% o f oxic domain) decreased when
compared to results obtained from lower ﬂux conditions, primarily because more O2 is moved
into the hyporheic pore waters. Figure 3.9b shows the dramatic impact of nutrients on reaction
limitations, where, not only does the fraction of microzones change (1% o f oxic domain), but
the spatial distribution of microzones is also altered. This result suggest that while domain
heterogeneity and hydraulic ﬂux are key parameters in microzone formation, unlike previous

35

modeling of hyporheic microzones [8], a dynamic DaO2 ought to also be considered in order to
completely analyze a hyporheic system.

(a)

(b)

Figure 3.9: a Oxygen proﬁle at 400 h for 1.0md−1 ﬂux and high nutrient concentrations and b.
Oxygen proﬁle at 700 h for 1.0md−1 ﬂux and low nutrient concnetrations.

3.4.2 Transferable Findings and Implications for Future Research

In summarizing our observations based on the simulation results, it becomes clear that there
are three main criteria that must be met in order for sandy hyporheic zones to develop and maintain
stable microzones, and, with them, their ecological functions. These criteria are as follows:

1. It must be an advection dominated hyporheic ﬂow systems.

2. There must be suﬃcient nutrient supplies to the hyporheic zone to stimulate oxygen uptake.

3. There cannot be stable bioclogging that forms in the hyporheic pores.

If these criteria are met, the actual quantum of stable microzones is likely to be a function of
the transport and reaction limitations that control the nutrient ﬂux into and distribution within the
hyporheic zone. However, in order to maximize the amount of anoxic microzones formed in bulk
oxic domain there must be a balance between the hydraulic ﬂux and inﬂux nutrient concentrations

36

such that there is also a sizable bulk oxic domain to harbor these anoxic microzones. Under low
ﬂux conditions, the formation of microzones is more reaction limited, where only under a low
nutrient concentration can a diﬀerence in oxygen concentration be obtained between a microzone
and the surrounding oxic domain. On the other hand, in the case of very high ﬂux conditions, the
microzones are ﬂux limited, where microzone numbers are diminished due to forcing of oxygen
into the anoxic pore spaces.

3.5 Conclusion

The formation of microzones in hyporheic sediments is a complex and highly non-linear
process, involving the interaction of hydraulic ﬂux, nutrient availability, and biomass growth
in the sediment matrix. This also results in variation in microzone abundance being dynamic
through space and time, so this indicates that the potential ecological functions of microzones, such
as greenhouse gas production, are also going to be dynamic through space and time. Although
bioclogging turns out to be detrimental to stable microzone formation, the likelihood of microzones
in nature, despite the presence of bioclogging, seems a real possibility. Sediment transport and
mobile bed conditions might naturally act to de-clog bioﬁlm that has accumulated at the inﬂuent of a
hyporheic zone, thereby facilitating the persistence of anoxic microzones and bulk oxic pore waters.
The simulation results for the microzones, suggest that while they may be small in extent relative
to bulk oxic or anoxic pore volumes, they can develop and persist. Thus, they are still a viable
hypothesis for the experimental and observation based studies that observe anaerobic products
occurring in oxic pore waters of hyporheic zones. Additional modeling and ﬁeld studies can use
the microzone formation criteria revealed in this study before attempting to upscale the overall
function of these microzones across diﬀerent types of hyporheic zones and across spatiotemporal
scales.

37

BIBLIOGRAPHY

38

BIBLIOGRAPHY

[1]

Joon Ho Ahn, Sungpyo Kim, Hongkeun Park, Brian Rahm, Krishna Pagilla, and Kartik
Chandran. N2O emissions from activated sludge processes, 2008-2009: Results of a national
monitoring survey in the united states. Environmental Science and Technology, 44(12):4505–
4511, 2010.

[2] T J Battin and D Sengschmitt. Linking Sediment Bioﬁlms, Hydrodynamics, and River Bed

Clogging: Evidence from a Large River. Microb Ecol, 37:185–196, 1999.

[3] Tom J Battin, Katharina Besemer, Mia M Bengtsson, Anna M Romani, and Aaron I Packmann.

The ecology and biogeochemistry of stream bioﬁlms. 2016.
Jacob. Bear. Dynamics of ﬂuids in porous media. American Elsevier Pub. Co, 1972.

[4]

[5] F. Boano, J. W. Harvey, A. Marion, A. I. Packman, R. Revelli, L. Ridolﬁ, and A. Wör-
man. Hyporheic ﬂow and transport processes: Mechanisms, models, and biogeochemical
implications. Reviews of Geophysics, 52(4):2012RG000417, dec 2014.

[6] H. (Hermann) Bothe, Stuart J. (Stuart John) Ferguson, and William E. (William Edward)

Newton. Biology of the nitrogen cycle. Elsevier, 2007.

[7] Martin A. Briggs, Frederick D. Day-Lewis, John B. T. Ong, Gary P. Curtis, and John W.
Lane. Simultaneous estimation of local-scale and ﬂow path-scale dual-domain mass transfer
parameters using geoelectrical monitoring. Water Resources Research, 49(9):5615–5630, sep
2013.

[8] Martin A. Briggs, Frederick D. Day-Lewis, Jay P. Zarnetske, and Judson W. Harvey. A physical
explanation for the development of redox microzones in hyporheic ﬂow. Geophysical Research
Letters, 42(11):4402–4410, jun 2015.

[9] Amy J. Burgin and Stephen K. Hamilton. Have we overemphasized the role of denitriﬁcation
in aquatic ecosystems? A review of nitrate removal pathways. Frontiers in Ecology and the
Environment, 5(2):89–96, mar 2007.

[10] Alice Caruso, Fulvio Boano, Luca Ridolﬁ, David L. Chopp, and Aaron Packman. Bioﬁlm-
induced bioclogging produces sharp interfaces in hyporheic ﬂow, redox conditions, and
microbial community structure. Geophysical Research Letters, 44(10):4917–4925, may 2017.
[11] Kartik Chandran, Lisa Y Stein, Martin G Klotz, and Mark C M van Loosdrecht. Nitrous
oxide production by lithotrophic ammonia-oxidizing bacteria and implications for engineered
nitrogen-removal systems. Biochemical Society transactions, 39(6):1832–7, dec 2011.

[12] C.c. Ezeuko, A. Sen, A. Grigoryan, and I.d. Gates. Pore-network modeling of bioﬁlm evolution

in porous media. Biotechnology and Bioengineering, 108(10):2413–2423, oct 2011.

39

[13] Chuanhui Gu, George M. Hornberger, Aaron L. Mills, Janet S. Herman, and Samuel A.
Flewelling. Nitrate reduction in streambed sediments: Eﬀects of ﬂow and biogeochemical
kinetics. Water Resources Research, 43(12), dec 2007.

[14] Charles F. Harvey, Roy Haggerty, and Steven M. Gorelick. Aquifer remediation: A method
for estimating mass transfer rate coeﬃcients and an evaluation of pulsed pumping. Water
Resources Research, 30(7):1979–1991, jul 1994.

[15] Judson W. Harvey, J. K. Böhlke, Mary A. Voytek, Durelle Scott, and Craig R. Tobias.
Hyporheic zone denitriﬁcation: Controls on eﬀective reaction depth and contribution to
whole-stream mass balance. Water Resources Research, 49(10):6298–6316, oct 2013.

[16] Judson W. Harvey, J. K. Böhlke, Mary A. Voytek, Durelle Scott, and Craig R. Tobias.
Hyporheic zone denitriﬁcation: Controls on eﬀective reaction depth and contribution to
whole-stream mass balance. Water Resources Research, 49(10):6298–6316, oct 2013.

[17] David. Knighton. Fluvial forms and processes : a new perspective. Arnold, 1998.
[18] K. Lansdown, C. M. Heppell, M. Trimmer, A. Binley, A. L. Heathwaite, P. Byrne, and
H. Zhang. The interplay between transport and reaction rates as controls on nitrate attenuation
in permeable, streambed sediments. Journal of Geophysical Research: Biogeosciences,
120(6):1093–1109, jun 2015.

[19] Yingyu Law, Bing-Jie Ni, Paul Lant, and Zhiguo Yuan. N 2 O production rate of an enriched
ammonia-oxidising bacteria culture exponentially correlates to its ammonia oxidation rate.
2012.

[20] Roxane Maranger, Stuart E. Jones, and James B. Cotner. Stoichiometry of carbon, nitrogen,
and phosphorus through the freshwater pipe. Limnology and Oceanography Letters, 3(3):89–
101, jun 2018.

[21] A. Marzadri, D. Tonina, A. Bellin, and J. L. Tank. A hydrologic model demonstrates ni-
trous oxide emissions depend on streambed morphology. Geophysical Research Letters,
41(15):5484–5491, aug 2014.

[22] Michael E. McClain, Elizabeth W Boyer, C Lisa Dent, Sarah E Gergel, Nancy B Grimm,
Peter M Groﬀman, Stephen C Hart, Judson W Harvey, Carol A Johnston, Emilio Mayorga,
William H. McDowell, and Gilles Pinay. Biogeochemical Hot Spots and Hot Moments at the
Interface of Terrestrial and Aquatic Ecosystems, 2003.

[23] F. J. Molz, M. A. Widdowson, and L. D. Beneﬁeld. Simulation of Microbial Growth Dynamics
Coupled to Nutrient and Oxygen Transport in Porous Media. Water Resources Research,
22(8):1207–1216, aug 1986.

[24] Patrick J. Mulholland, Erich R. Marzolf, Jackson R. Webster, Deborah R. Hart, and Susan P.
Hendricks. Evidence that hyporheic zones increase heterotrophic metabolism and phosphorus
uptake in forest streams. Limnology and Oceanography, 42(3):443–451, may 1997.

40

[25] Ignacio Peralta-Maraver, Jason Galloway, Malte Posselt, Shai Arnon, Julia Reiss, Jörg
Lewandowski, and Anne L. Robertson. Environmental ﬁltering and community delineation
in the streambed ecotone. Scientiﬁc Reports, 8(1):15871, dec 2018.

[26] Annika M Quick, W. Jeﬀery Reeder, Tiﬀany B Farrell, Daniele Tonina, Kevin P Feris, and
Shawn G Benner. Controls on Nitrous Oxide Emissions from the Hyporheic Zones of Streams.
Environmental Science and Technology, 50(21):11491–11500, 2016.

[27] A. Raoof, H.M. Nick, S.M. Hassanizadeh, and C.J. Spiers. PoreFlow: A complex pore-
network model for simulation of reactive transport in variably saturated porous media. Com-
puters & Geosciences, 61:160–174, dec 2013.

[28] Michael O. Rivett, Stephen R. Buss, Philip Morgan, Jonathan W N Smith, and Chrystina D.
Bemment. Nitrate attenuation in groundwater: A review of biogeochemical controlling
processes, 2008.

[29] A. H. Sawyer. Enhanced removal of groundwater-borne nitrate in heterogeneous aquatic

sediments. Geophysical Research Letters, 42(2):2014GL062234, jan 2015.

[30] R.R. Sharp, A.B. Cunningham, J. Komlos, and J. Billmayer. Observation of thick bioﬁlm
accumulation and structure in porous media and corresponding hydrodynamic and mass
transfer eﬀects. Water Science and Technology, 39(7):195–201, jan 1999.

[31] Will Steﬀen, Katherine Richardson, Johan Rockström, Sarah E Cornell, Ingo Fetzer, Elena M
Bennett, Reinette Biggs, Stephen R Carpenter, Wim de Vries, Cynthia A de Wit, Carl Folke,
Dieter Gerten, Jens Heinke, Georgina M Mace, Linn M Persson, Veerabhadran Ramanathan,
Belinda Reyers, and Sverker Sörlin. Sustainability. Planetary boundaries: guiding human
development on a changing planet. Science (New York, N.Y.), 347(6223):1259855, feb 2015.
[32] Deborah L Stoliker, Deborah A Repert, Richard L Smith, Bongkeun Song, Denis R. LeBlanc,
Timothy D. McCobb, Christopher H Conaway, Sung Pil Hyun, Dong Chan Koh, Hee Sun
Moon, and Douglas B Kent. Hydrologic Controls on Nitrogen Cycling Processes and Func-
tional Gene Abundance in Sediments of a Groundwater Flow-Through Lake. Environmental
Science and Technology, 50(7):3649–3657, 2016.

[33] Stewart W. Taylor and Peter R. Jaﬀé. Bioﬁlm growth and the related changes in the physical
properties of a porous medium: 3. Dispersivity and model veriﬁcation. Water Resources
Research, 26(9):2171–2180, sep 1990.

[34] Martin Thullner. Comparison of bioclogging eﬀects in saturated porous media within one-and

two-dimensional ﬂow systems. Ecological Engineering, 36:176–196, 2010.

[35] Martin Thullner and Philippe Baveye. Computational pore network modeling of the inﬂuence
of bioﬁlm permeability on bioclogging in porous media. Biotechnology and Bioengineering,
99(6):1337–1351, 2008.

[36] Martin Thullner, Josef Zeyer, and Wolfgang Kinzelbach. Inﬂuence of microbial growth on

hydraulic properties of pore networks. Transport in Porous Media, 49(1):99–122, 2002.

41

[37] Mark A. Widdowson, Fred J. Molz, and Larry D. Beneﬁeld. A numerical transport model for
oxygen- and nitrate-based respiration linked to substrate and nutrient availability in porous
media. Water Resources Research, 24(9):1553–1565, sep 1988.

[38] Jost Wingender, Thomas R. Neu, and Hans-C. (Hans-Curt) Flemming. Microbial extracellular

polymeric substances : characterization, structure, and function. Springer, 1999.

[39] Thomas C Winter, Judson W Harvey, O Lehn Franke, and William M Alley. Ground Water

and Surface Water A Single Resource. 1998.

[40] Jay P Zarnetske, Roy Haggerty, and Steven M Wondzell. Coupling multiscale observations
to evaluate hyporheic nitrate removal at the reach scale. Freshwater Science, 34(1):172–186,
2015.

[41] Jay P. Zarnetske, Roy Haggerty, Steven M. Wondzell, and Michelle A. Baker. Dynamics of
nitrate production and removal as a function of residence time in the hyporheic zone. Journal
of Geophysical Research, 116(G1):G01025, mar 2011.

[42] Jay P. Zarnetske, Roy Haggerty, Steven M. Wondzell, Vrushali A. Bokil, and Ricardo
González-Pinzón. Coupled transport and reaction kinetics control the nitrate source-sink
function of hyporheic zones. Water Resources Research, 48(11), 2012.

[43] Xia Zhu, Martin Burger, Timothy A Doane, and William R Horwath. Ammonia oxidation
pathways and nitriﬁer denitriﬁcation are signiﬁcant sources of N2O and NO under low oxygen
availability. Proceedings of the National Academy of Sciences of the United States of America,
110(16):6328–33, apr 2013.

42