U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2008–5066

National Water-Quality Assessment Program 

Hypothetical Modeling of Redox Conditions Within a 

Complex Ground-Water Flow Field in a Glacial Setting

4.000e-002

0.400

4.000

40.00

0 2,000 METERS

Concentration mg/L,

log scale

EXPLANATION

Hypothetical Modeling of Redox Conditions 

Within a Complex Ground-Water Flow Field 

in a Glacial Setting

By Daniel T. Feinstein and Mary Ann Thomas

Scientific Investigations Report 2008–5066 

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia: 2009

For more information on the USGS—the Federal source for science about the Earth, its natural and living resources, 

natural hazards, and the environment, visit http://www.usgs.gov or call 1-888-ASK-USGS

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Feinstein, D.T., and Thomas, M.A., 2009, Hypothetical modeling of redox conditions for a complex ground-water flow 

field in a glacial setting: U.S. Geological Survey Scientific Investigations Report 2008–5066, 27 p.

Cover. Map showing simulated concentration of dissolved Fe2+ for random organic carbon 

scenario, lower aquifer (layer 3).
iii

Foreword

The U.S. Geological Survey (USGS) is committed to providing the Nation with credible scientific 

information that helps to enhance and protect the overall quality of life and that facilitates 

effective management of water, biological, energy, and mineral resources (http://www.usgs.

gov/). Information on the Nation’s water resources is critical to ensuring long-term availability 

of water that is safe for drinking and recreation and is suitable for industry, irrigation, and fish 

and wildlife. Population growth and increasing demands for water make the availability of that 

water, now measured in terms of quantity and quality, even more essential to the long-term 

sustainability of our communities and ecosystems.

The USGS implemented the National Water-Quality Assessment (NAWQA) Program in 1991 

to support national, regional, State, and local information needs and decisions related to 

water-quality management and policy (http://water.usgs.gov/nawqa). The NAWQA Program is 

designed to answer: What is the condition of our Nation’s streams and ground water? How are 

conditions changing over time? How do natural features and human activities affect the quality 

of streams and ground water, and where are those effects most pronounced? By combining 

information on water chemistry, physical characteristics, stream habitat, and aquatic life, the 

NAWQA Program aims to provide science-based insights for current and emerging water issues 

and priorities. From 1991-2001, the NAWQA Program completed interdisciplinary assessments 

and established a baseline understanding of water-quality conditions in 51 of the Nation’s river 

basins and aquifers, referred to as Study Units (http://water.usgs.gov/nawqa/studyu.html). 

In the second decade of the Program (2001–2012), a major focus is on regional assessments 

of water-quality conditions and trends. These regional assessments are based on major river 

basins and principal aquifers, which encompass larger regions of the country than the Study 

Units. Regional assessments extend the findings in the Study Units by filling critical gaps in 

characterizing the quality of surface water and ground water, and by determining status and 

trends at sites that have been consistently monitored for more than a decade. In addition, the 

regional assessments continue to build an understanding of how natural features and human 

activities affect water quality. Many of the regional assessments employ modeling and other 

scientific tools, developed on the basis of data collected at individual sites, to help extend 

knowledge of water quality to unmonitored, yet comparable areas within the regions. The 

models thereby enhance the value of our existing data and our understanding of the hydrologic 

system. In addition, the models are useful in evaluating various resource-management scenarios 

and in predicting how our actions, such as reducing or managing nonpoint and point sources of 

contamination, land conversion, and altering flow and (or) pumping regimes, are likely to affect 

water conditions within a region.

Other activities planned during the second decade include continuing national syntheses of 

information on pesticides, volatile organic compounds (VOCs), nutrients, selected trace elements, and aquatic ecology; and continuing national topical studies on the fate of agricultural 

chemicals, effects of urbanization on stream ecosystems, bioaccumulation of mercury in stream 

ecosystems, effects of nutrient enrichment on stream ecosystems, and transport of contaminants to public-supply wells.
iv

The USGS aims to disseminate credible, timely, and relevant science information to address 

practical and effective water-resource management and strategies that protect and restore 

water quality. We hope this NAWQA publication will provide you with insights and information 

to meet your needs, and will foster increased citizen awareness and involvement in the protection and restoration of our Nation’s waters. 

The USGS recognizes that a national assessment by a single program cannot address all waterresource issues of interest. External coordination at all levels is critical for cost-effective management, regulation, and conservation of our Nation’s water resources. The NAWQA Program, 

therefore, depends on advice and information from other agencies—Federal, State, regional, 

interstate, Tribal, and local—as well as nongovernmental organizations, industry, academia, and 

other stakeholder groups. Your assistance and suggestions are greatly appreciated.

Matthew C. Larsen

Associate Director for Water
v

Contents

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................1

Methods and Assumptions...........................................................................................................................1

Flow Model.............................................................................................................................................2

Geochemical Model .............................................................................................................................3

Equations of the Geochemical Model ......................................................................................5

Initial Geochemical Conditions..................................................................................................6

Definition of Redox Zones ..........................................................................................................9

Hypothetical Scenarios—A Single Ground-Water Flow Field with Varied Distributions of 

Organic Carbon ..............................................................................................................................10

Ground-Water Flow System..............................................................................................................10

Organic Carbon Source Areas ........................................................................................................10

Random Source Area.........................................................................................................................12

Lacustrine Source Area ...........................................................................................................12

Wetland Source Area................................................................................................................12

Results of Hypothetical Scenarios............................................................................................................12

Equilibrium Organic Carbon Concentrations..................................................................................12

Equilibrium Redox Water-Chemistry Conditions............................................................................15

Sensitivity Analysis......................................................................................................................................19

Sensitivity to Solver Strategy............................................................................................................20

Sensitivity to Attenuation Parameters ............................................................................................20

Sensitivity to Concentrations of Electron Acceptors....................................................................20

Sensitivity to Electron-Donor Behavior...........................................................................................22

Influence of Pumping on Distribution of Reduced Waters ...................................................................23

Limitations of the Study ..............................................................................................................................24

Conclusions...................................................................................................................................................25

Acknowledgments.......................................................................................................................................25

References Cited..........................................................................................................................................25

Figures

1. Diagrammatic geologic section of the Susquehanna River Valley, 

Broome County, N.Y......................................................................................................................2

2. Plan view of model grid and surface-water network.............................................................4

3. Map showing simulated water-table elevation.......................................................................5

4. Graph showing relation between concentration and depletion rate of dissolved 

organic carbon for decay constant equal to 0.001 day-1 ........................................................8

5–13. Maps showing:

5. Model-computed flow velocities in upper and lower aquifers..................................11

6. Spatial distribution of ground-water age in lower aquifer..........................................12

7. Electron-donor source areas and simulated electron donor concentrations 

in source layer for three organic carbon source scenarios ......................................13
vi

8. Simulated dissolved organic carbon concentrations in the lower aquifer 

for three source scenarios...............................................................................................14

9. Simulated redox water chemistry in the lower aquifer (layer 3) for three 

organic carbon source scenarios...................................................................................16

10. Comparison of A, simulated redox reaction zones to B, simulated redox 

water-chemistry zones for random organic carbon scenario, lower aquifer 

(layer 3)................................................................................................................................18

11. Simulated concentration of dissolved Fe2+ for random organic carbon 

scenario, lower aquifer (layer 3).....................................................................................19

12. Random organic carbon scenario for sensitivity analysis (Csource = 100 mg/L); 

simulated redox water-chemistry conditions in lower aquifer (layer 3) ..................20

13. Simulated drawdown contours (meters) around pumping wells ..............................23

Tables

1. Concentrations of dissolved organic carbon in varied hydrogeologic settings ................6

2. Concentrations of electron acceptors and electron donors in ground-water 

samples from the glacial aquifer system in the northern United States, 1991–2002.........7

3. Assumed initial reactant concentrations and reaction yield coefficients..........................8

4. Definition of redox water-chemistry zones in terms of dissolved species .........................9

5. Simulated redox water-chemistry zones in the upper aquifer, confining unit, 

and lower aquifer by percent volume .....................................................................................17

6. Mass flux of redox constituents for organic carbon source scenarios for

 entire model................................................................................................................................17

7. Results of sensitivity simulations.............................................................................................21

8a. Simulated degradation rates of redox reactants for mobile versus fixed electron 

donor source ...............................................................................................................................22

8b. Simulated production rates of redox products for mobile versus fixed electron 

donor source ...............................................................................................................................22

9. Effect of pumping on concentrations of redox reactants for the three 

hypothetical scenarios ..............................................................................................................24
vii

Conversion Factors and Abbreviations

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

Flow rate

cubic meter per second (m3/s) 22.82 million gallons per day (Mgal/d)

Hydraulic conductivity and dispersivity

meter per day (m/d) 3.281 feet per day (ft/d)

Concentrations of chemical constituents are given in milligrams per liter (mg/L). Production and 

depletion rates of redox reactants are given in milligrams per liter per day (mg/L/d).

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29).

Elevation, as used in this report, refers to distance above the vertical datum.

Hypothetical Modeling of Redox Conditions Within a 

Complex Ground-Water Flow Field in a Glacial Setting

By Daniel T. Feinstein and Mary Ann Thomas

Abstract

This report describes a modeling approach for studying 

how redox conditions evolve under the influence of a complex 

ground-water flow field. The distribution of redox conditions 

within a flow system is of interest because of the intrinsic susceptibility of an aquifer to redox-sensitive, naturally occurring 

contaminants—such as arsenic—as well as anthropogenic contaminants—such as chlorinated solvents. The MODFLOWMT3D-RT3D suite of code was applied to a glacial valley-fill 

aquifer to demonstrate a method for testing the interaction of 

flow patterns, sources of reactive organic carbon, and availability of electron acceptors in controlling redox conditions. 

Modeling results show how three hypothetical distributions of 

organic carbon influence the development of redox conditions 

in a water-supply aquifer. The distribution of strongly reduced 

water depends on the balance between the rate of redox reactions and the capability of different parts of the flow system 

to transmit oxygenated water. The method can take account of 

changes in the flow system induced by pumping that result in a 

new distribution of reduced water.

Introduction

Redox conditions of ground water influence the mobilization and transport of chemical constituents such as arsenic 

(Welch and others, 2000), chlorinated solvents (Moran and 

others, 2007), and pesticides (Barbash, 2004). An increased 

understanding of the distribution of redox conditions within an 

aquifer would allow for a better assessment of its susceptibility to contamination by a wide range of natural and anthropogenic constituents. 

The current study is an outgrowth of a USGS National 

Water-Quality Assessment (NAWQA) regional study of 

arsenic in the glacial aquifer system of the northern United 

States (Thomas, 2007). One of the major findings was that the 

highest arsenic concentrations were associated with reducing 

ground water. The crucial importance of the relation between 

the mobility of contaminants such as arsenic and the degree 

to which ground water is reduced invites a research effort 

aimed at combining principles of subsurface chemistry with 

fine-scale simulations of ground-water flow fields in order to 

deepen what is known about the controls on redox conditions 

in aquifer systems.

The distribution of redox zones in the subsurface is a 

product of microbe-mediated reactions, and it is influenced 

by the interaction between the spatially varying ground-water 

flow field and the availability of electron donors and electron 

acceptors in the subsurface. Reducing conditions can be maintained only if the advection and diffusion of dissolved oxygen 

from the surface are less rapid than the consumption of oxygen 

by redox reactions (Smedley and Kinniburgh, 2002). 

In this report, the evolution of redox is evaluated in a 

glacial valley-fill aquifer characterized by complex flow paths, 

pronounced heterogeneity, and streams that are both gaining 

and losing. The primary objective of this study is to show the 

importance of the ground-water flow field in determining the 

spread and distribution of reduced waters in a glacial valleyfill aquifer. Although the flow system used to test the methodology is real, the geochemical conditions and geochemical 

simulations are hypothetical. A methodology for modeling 

redox has been applied to three hypothetical distributions of 

organic carbon and the robustness of the methodology has 

been tested through a sensitivity analysis of one scenario. 

Availability of a pumping history within the flow system also 

enabled investigation of the possible effect of stresses on the 

redistribution of reduced water emanating from the hypothetical source areas. 

Methods and Assumptions

The approach to the problem of simulating redox reactions in a field setting was to combine a heterogeneous, 

three-dimensional ground-water flow model—which incorporates a sophisticated geologic model of glacial deposition in 

bedrock valleys—with a relatively simple geochemical model 

for redox reactions. The MODFLOW-MT3D-RT3D suite 

of code (McDonald and Harbaugh, 1988; Zheng and Wang, 

1998; Clement, 2001) was used for this purpose. Although the 

flow model was calibrated to real conditions, the geochemical 

inputs to the transport model were not derived from local data 

but instead are representative of water-quality data from the 
2 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

glacial aquifer system of the northern United States, collected 

as part of the National Water-Quality Assessment (NAWQA) 

Program of the U.S. Geological Survey. The geometries of the 

organic-rich material horizons that function as electron donors 

for the modeled redox reactions are hypothetical.

Flow Model

The ground-water flow model adopted for this study is 

based on a three-facies interpretation of glacial deposition 

in the eastern United States (Randall, 2001). The conceptual 

model is applicable to bedrock valley deposits (also called 

valley-train deposits) that commonly contain sufficient coarse 

material from ice-contact and outwash deposits to support 

high-capacity public-supply wells. In Randall’s conceptual 

model, based on extensive field mapping across the northeastern United States, the ice-contact deposits form the bottom of 

the system and are separated from more recent outwash (and 

alluvial) deposits at the top of the system by generally finegrained layers of lacustrine origin. However, the geometry of 

the three facies can be complex. 

Randall incorporated the three-facies geologic model into 

a ground-water flow model of an important aquifer system in 

southern New York State associated with the Susquehanna 

River Valley (Randall, 1986); this original model was later 

upgraded by Wolcott and Coon (2001). A typical geologic section from the ground-water flow model shows the complicated 

geometry of the facies within the bedrock valley (fig. 1). 

The aquifer system consists of an upper aquifer containing outwash and alluvial deposits (mostly sand and gravel) and 

a lower aquifer dominated by ice-contact deposits (also generally coarse grained), separated by confining lenses of silt and 

clay of lacustrine origin. The thickness of the confining layer 

is varied and is locally absent, putting the two aquifers locally 

in direct hydraulic connection. The system is characterized by 

a pronounced heterogeneity due to wide ranges in horizontal 

and vertical hydraulic conductivity. In the model, the median 

horizontal hydraulic conductivity of the upper aquifer in layer 

1 is 40 m/d, with 95 percent of cells assigned values between 

2 and 70 m/d. The median value for the lacustrine material in 

layer 2 is 0.002 m/d, with 95 percent of cells between 0.0002 

and 0.5 m/d. The median value for the lower aquifer in layer 

Figure 1. Diagrammatic geologic section of the Susquehanna River Valley, Broome County, N.Y. [Modified from Randall (1986) 

and Wolcott and Coon (2001)].

West Cross section along row 24 East

Vertical exaggeration = 100X

EXPLANATION

Alluvium and glacial outwash (upper aquifer, model layer 1)

Glacial lacustrine deposits (model layer 2)

Glacial ice-contact deposits (lower aquifer, model layer 3)

Shale

0 2,000 METERS

Location of area shown below

L a k e E ri e

LakeOntario

OHIO

PENNSYLVANIA

Broome 

County

NEW YORK

VT

MAINE

NH

MA

CT

NJ

RI

0 100 MILES
Methods and Assumptions 3

3 is 66 m/d, with 95 percent of cells between 4 and 400 m/d. 

The hydraulic conductivity pattern for the upper aquifer is 

particularly heterogeneous. Hydraulic conductivities for the 

aquifer layers were based on about 260 geologic logs, specificcapacity tests from about 80 large-capacity or test wells, and 

pumping-test results from 20 wells (Randall, 1986). Hydraulic 

conductivities for the lacustrine material were estimated from 

observed texture and literature values (Randall, 1986). 

Ground water originates from three sources: recharge to 

the water table, lateral infiltration due to unchanneled runoff 

from the surrounding shale-covered uplands, and natural or 

induced infiltration along losing stretches of streams, including parts of the Susquehanna River. 

Figure 2 shows the model grid and surface-water network in plan view. The model grid contains 61 rows and 157 

columns for each layer, with nonuniform cell spacing. Average cell size is 100 m on a side. Recharge due to precipitation is added at the water table, whereas lateral infiltration is 

simulated by constant-flux boundaries. Rivers and streams are 

head-dependent boundaries. In some places, the surface water 

loses to the ground water. The contact between the unconsolidated glacial material and the underlying shale is a no-flow 

boundary.

The ground-water flow solution is solved by applying 

the flow-model code MODFLOW (McDonald and Harbaugh, 

1988). The water table generally slopes from the valley walls 

toward the surface-water network (fig. 3).

Geochemical Model

The geochemical model adopted for this study posits that 

reducing conditions occur in ground water because microbes 

gain energy from redox reactions involving electron donors 

and electron acceptors. Naturally occurring electron donors 

include reactive organic carbon derived from plant decay and 

sulfide minerals such as pyrite. For the hypothetical simulations, the electron donor is assumed to be organic carbon. The 

electron acceptors are dissolved oxygen, dissolved nitrate, 

solid-phase ferric iron, dissolved sulphate, and dissolved 

carbon dioxide. Typically, the microbes utilize the electron 

acceptors sequentially (Chapelle and others, 1995; Guoping 

and others, 1999; Schreiber and others, 2004). 

The geochemical model also requires the assumption that 

these reactions approximately follow a first-order law with 

respect to the electron donor and zero order with respect to the 

electron acceptors, implying that the rate of reaction between 

the electron donor and any electron acceptor is directly proportional to the concentration of the electron donor. Although 

more sophisticated Monod-type kinetic models could have 

been used (Bekins and others, 1998; Rolle and others, 2006), 

the difficulty in quantifying the parameters necessary for a 

Monod formulation led to adoption of the simpler model. 

This simple geochemical model does not directly simulate 

the reaction of the electron donor with solid-phase ferric iron, 

Fe3+. Instead, the assumption is that once a given amount of 

dissolved ferrous iron, Fe2+, is produced, then either the Fe3+

on sediments is no longer available to be reduced or (as is 

more likely) iron reduction continues but any additional Fe2+

derived from the solid phase precipitates out of solution or 

is sorbed to solid-phase minerals. At this point, it is assumed 

that sulfate reduction replaces iron reduction even though it is 

likely that iron reduction continues to occur (but at too slow 

a rate to offset ferrous iron precipitation and thereby increase 

the concentration of dissolved Fe2+). 

The sequential first-order geochemical model incorporates four key assumptions. First, the electron donor is 

represented by the simplest organic carbon molecule, CH2O, 

and the stoichiometry of the reactions corresponds to the degradation of this molecule. An alternative electron donor, such 

as pyrite, can replace CH2O in the redox reactions with little 

change in the reaction stoichiometry. Second, the availability 

of electron donors in local source areas (associated with individual model cells) is assumed to be unlimited, and the donors 

are assumed to move freely from the solid to liquid phase in 

those areas. Third, electron donors are assumed to be mobile; 

that is, dissolved species in ground water are able to move 

downgradient from the source areas and continue to react 

with dissolved electron acceptors according to the posited 

first-order rate law. (This assumption was relaxed for some 

sensitivity simulations.) Finally, sources and sinks for electron 

acceptors other than recharge and river water are neglected, 

even though transfers within the subsurface can at times be 

important. The transient effect of water-table fluctuations that 

sometimes influences the availability of electron acceptors is 

also neglected. Within the limits imposed by these simplifying assumptions, the methodological approach is designed to 

simulate the crucial interaction between the rate of influx of 

young oxygenated water and the rates of redox reaction in the 

presence of a variable velocity field.

The transport code MT3D (Chunmiao and Wang, 1998), 

which is linked to MODFLOW, incorporates advection, dispersion, and sink/sources for the electron donor and acceptors. 

The reaction code RT3D (Clement, 2001), which simulates 

reactive multispecies transport, is, in turn, linked to MT3D. 

We used RT3D to simulate strictly sequential first-order 

decay reactions between a single electron donor (naturally 

occurring organic carbon) and five electron acceptors via five 

pathways: aerobic respiration, denitrification, iron reduction, 

sulfate reduction, and methanogenesis. Reaction pathways 

switch from favored to less favored when the concentration 

of the favored electron acceptor falls to near zero. For example, when the dissolved oxygen concentration in a model 

cell falls below 0.01 mg/L, the redox reaction switches from 

aerobic respiration to denitrification. Similarly, when the 

nitrate concentration falls below 0.01 mg/L, the redox reaction 

switches from denitrification to iron reduction. However, the 

electron acceptor Fe3+ is not explicitly simulated in the routine 

because it is a solid-phase constituent. Therefore, the switching from iron reduction to sulfate reduction in the RT3D code 

is triggered by the buildup of the reaction product, Fe2+, to a 

maximum concentration that reflects iron reduction and the 
4 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting Figure 2. Plan view of model grid and surface-water network (green cells). Black cells are inactive. 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 102030405060

0 2,000 METERS EXPLANATION

Surface-water networkInactive cells

ROWS

COLUMNS
Methods and Assumptions 5

effect of sorption and precipitation reactions that remove Fe2+

from the dissolved phase. Finally, sulfate reduction conditions 

yield to methanogenesis when the sulfate concentration falls 

below 0.01 mg/L. 

The transient transport simulations begin with background conditions for electron acceptors that correspond to the 

concentrations in recharge and river water. The simulations are 

run until equilibrium conditions are attained with respect to all 

reactants in all model cells. 

Redox conditions become stable when the rate at which 

electron acceptors are supplied by incoming water balances 

the rate at which they are consumed by redox reactions.

Equations of the Geochemical Model

Module 3 of the RT3D code (Clement, 2001) was adapted 

to model the distribution of reduced ground water. The module 

simulates kinetic-limited degradation by way of multiple 

electron acceptors. It was intended originally to model BTEX 

degradation (Clement and others, 1997), but it was used 

in this study to model degradation of any electron donor, 

such as CH2O, via the five degradation pathways mentioned 

previously.

RT3D Module 3 contains one equation for each of the 

five reaction pathways. In this study, the equations were 

applied in a simplified form and effectively eliminated 

inhibition and Monod half-saturation variables. The resulting expressions are strictly sequential—with no overlapping 

reactions at a given location at a given time—and they are 

first-order with respect to the electron donor and zero-order 

with respect to the electron acceptors. The reformulated equations are as follows:

rDONOR,O2 =-kO2 × [DONOR] for [O2] ≥ 0.01 mg/L

rDONOR,NO3=-kNO3 × [DONOR] for [O2] < 0.01 mg/L and 

[NO3] ≥ 0.01 mg/L

rDONOR,Fe2+=-kFe3+ × [DONOR] for [O2] < 0.01 mg/L, [NO3] 

< 0.01 mg/L and [Fe3+] ≥ 0.01 mg/L

rDONOR,SO4=-kSO4 × [DONOR] for [O2] < 0.01 mg/L, [NO3] < 

0.01 mg/L, [Fe3+] < 0.01 mg/L and [SO4] ≥ 0.01 mg/L

rDONOR,CH4=-kCO2 × [DONOR] for [O2] < 0.01 mg/L, [NO3] < 

0.01 mg/L, [Fe3+] < 0.01 mg/L, [SO4] < 0.01 mg/L, and [CO2] 

≥ 0.01 mg/L

where

[ ] indicates concentration, in milligrams per liter, 

rDONOR is the destruction rate of the electron donor (for 

example, organic carbon), and 

k O2= kNO3= kFe3+=kSO4=kCO2 is a single first-order reaction 

constant.

Figure 3. Simulated water-table elevation.

River cells

Water-level elevation,

in meters above NGVD 29

240.0

245.0

250.0

255.0

260.0

Layer 1 

0 2,000 METERS
6 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Because the concentration of Fe3+ and CO2 are not dissolved species and cannot be measured under normal field 

conditions (Wiedermeier and others, 1995), these concentrations are replaced in module 3 of RT3D with “assimilative 

capacity terms” defined as

[Fe3+] = [Fe2+max] – [Fe2+]

[CO2] = [CH4max] - [CH4], 

where 

the concentration term for Fe3+ and CO2 should be considered as a hypothetical term that represents the iron 

reduction and methanogenic capacity of the system 

and where [Fe2+max] and [CH4max] are the maximum 

concentrations of Fe2+ and CH4 allowed to occur in the 

model. 

The RT3D code simulates the consumption of electron 

acceptor as well as electron donor for any given redox reaction. It applies the following equations:

rO2 =YO2/DONOR × rDONOR,O2

rNO3 =YNo3/DONOR × rDONOR,NO3

r

Fe2+) =-YFe2+/DONOR ×* rDONOR,Fe2+

rSO4 =YSO4/DONOR × rDONOR.,SO4

r

CH4 =-YCH4/DONOR ×* rDONOR,CH4

where 

Y is the appropriate yield coefficient in each chemical 

reaction; that is, the mass ratio of electron acceptor 

consumed to electron donor (for example, CH2O) 

consumed.

Because the reactions are strictly sequential, only one 

electron acceptor is consumed at any location (model cell). 

The active reaction at a location can change over time, but 

only one reaction is permanently active once equilibrium 

conditions are achieved. All reactions are assumed to occur in 

the aqueous phase.

Initial Geochemical Conditions

The model requires concentrations of electron donors 

and electron acceptors for recharge, river water, and initial 

ground water. The electron donor was assumed to be organic 

carbon. Input values for organic carbon and electron acceptors 

are hypothetical but are generally consistent with literature 

values (table 1) and with water-quality data collected from the 

northern U.S. glacial aquifer system as part of NAWQA studies during 1991–2002 (table 2). 

The literature values reflect a range of dissolved organic 

carbon concentrations in surficial settings and in carbon-rich 

Table 1. Concentrations of dissolved organic carbon in varied hydrogeologic settings. 

[All data in milligrams per liter; <, less than; --, not available; bold values indicate median value]

Site description Range Mean or median Source

Water in equilibrium with atmosphere -- 10 Drever, 2002

Soil solutions:

O/A horizon 21–32 -- Drever, 2002

B horizon 5–7 -- Drever, 2002

Ground water Most <2 0.7 Drever, 2002

Rivers 2–15 4–6 Drever, 2002

Waters draining wetlands 5–60 25 Drever, 2002

Lakes:

Low productivity 1–3 -- Drever, 2002

Eutrophic 2–5 -- Drever, 2002

Shallow ground water in glacial outwash, Wisconsin Up to 25 4 Keating and Bahr, 1998

Clay-rich confining unit, southern Alberta -- 16 (typical) Hendry and Wassenaar, 2005

Boom clay, Belguim 100–500 -- Hendry and Wassenaar, 2005

Mexico City clay plain confining unit 300–580 -- Hendry and Wassenaar, 2005
Methods and Assumptions 7

Table 2. Concentrations of electron acceptors and electron donors in ground-water samples from the glacial aquifer system in 

the northern United States, 1991–2002.

[mg/L, milligrams per liter; <, less than]

Acceptor or donor All samples1 Samples from the 

northeastern United States

Electron-acceptor concentrations

Dissolved oxygen

Median (mg/L) 1.5 3.0

Range (mg/L) <0.1–14 <0.1–13

Number of samples 799 170

Nitrate

Median (mg/L) 0.5 0.8

Range (mg/L) <0.06–77 <0.06–26

Number of samples 800 173

Iron

Median (mg/L) 0.015 <0.010

Range (mg/L) <0.010–38 <.0.010–38

Number of samples 810 172

Sulfate

Median (mg/L) 30 13

Range (mg/L) <0.3–2,600 <0.3–340

Number of samples 811 172

Electron-donor concentrations

Dissolved organic carbon

Median (mg/L) 1.0 0.6

Range (mg/L) 0.2–260 0.2–14

Number of samples 769 172

1

 Samples collected by the National Water-Quality Assessment (NAWQA) Program during 1991–2003 and used for analysis in Thomas (2007).
8 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

subsurface areas (for example, waters draining wetlands and 

some confining units) but show concentrations consistently 

close to 1.0 mg/L in most shallow ground water (table 1). The 

NAWQA glacial-aquifer database indicates some variability 

in the central tendency of electron-acceptor and dissolved 

organic carbon concentrations, especially for dissolved nitrate 

and sulfate (table 2). The values selected for the hypothetical 

simulations (table 3) are generally in line with the median values from the NAWQA glacial-aquifer dataset for the northeastern United States and are considered representative of shallow, 

recently recharged ground water.

Reactive organic carbon is inserted into the model in 

two ways: as diffuse background loadings, and as elevated 

concentrations in source areas, such as buried wetlands. The 

diffuse loading occurs mostly through recharge but also from 

losing river reaches. The organic carbon concentration in 

recharge, river water, and initial background ground water was 

set constant at 1 mg/L, which is intended to represent typical 

conditions for shallow ground water. As the transport simulation proceeds, the organic carbon concentrations outside the 

source areas are free to change with the arrival of upgradient 

water. Redox reactions occur until an equilibrium is reached 

at any non-source location between the decay of the organic 

carbon and the amount replenished by incoming ground-water 

flow in the same period of time. The elevated organic carbon 

concentrations assigned to the source areas are discussed in a 

later section of the report. 

Other inputs to the reaction model in MT3D/RT3D 

include the decay-rate constants for each sequential reaction, 

the reaction mass ratios, the maximum allowed Fe2+ and CH4

concentration, and the solver scheme. 

For the primary simulations, a first-order rate constant of 

0.001 day-1 was applied to all redox reactions. The resulting 

depletion rate of CH2O as a function of CH2O concentration is 

shown in fig. 4. 

The selected first-order decay rate falls between a 

reported high of 0.1 day-1 for deep-sea sediments and a 

reported low of 0.000001 day-1, for regional aquifer systems in 

the United States with long flow paths (derived from Chapelle, 

1993, fig. 7.6, assuming that the electron-donor compound is 

CH20). It seems reasonable that the decay rate in a shallow 

glacial system would fall between these end-members. The 

use of a single decay rate for all electron acceptors is, in part, 

a simplification in the absence of hard data on individual rates. 

However, this adapted convention is also consistent with findings by Bohlke and others (2002) that first-order rate constants 

in glacial sand aquifers were similar for oxygen reduction, 

nitrate reduction, and sulfate reduction. 

Figure 4. Relation between 

concentration and depletion 

rate of dissolved organic 

carbon for decay constant 

equal to 0.001 day-1.

Table 3. Assumed initial reactant concentrations 

and reaction yield coefficients.

[mg/L, milligrams per liter]

Reactant Concentration 

(mg/L)

Yield 

coefficient

Organic carbon, as CH2O 1 1

O2 6 1.07

NO3 0.3 1.65

Fe2+ 0.005 7.44

SO4 10 1.60

CH4 0 0.267

0.001 0.01 0.1 1 10 100

CONCENTRATION OF DISSOLVED ORGANIC CARBON,

IN MILLIGRAMS PER LITER

DEPLETION RATE OF DISSOLVED ORGANIC CARBON,

IN MILLIGRAMS PER LITER PER DAY 

0.000001

0.00001

0.0001

0.001

0.01

0.1
Methods and Assumptions 9

The yield coefficients for redox reactions between CH2O 

and electron acceptors are listed in table 3. They account for 

the consumption of electron acceptors during redox reactions. For example, if 10 mg/L of dissolved organic carbon 

is present in a model cell and the concentration of oxygen is 

greater than 0.01 mg/L, then the depletion rate of the organic 

carbon is 0.001 x 10 or 0.01 mg/L/d and the depletion rate of 

dissolved oxygen, using its yield coefficient, is 1.07 x 0.01 or 

0.0107 mg/L/d. If dissolved oxygen and nitrate concentrations 

are less than 0.01 mg/L and the concentration of Fe2+ is less 

than 40 mg/L, then the depletion rate of the organic carbon is 

still 0.01 mg/L/d but the production rate of Fe2+ is 7.44 x 0.01 

or 0.0744 mg/L/d.

The maximum allowable Fe2+ was set at 40 mg/L. This 

is the concentration at which modeled redox reactions switch 

from iron reducing to sulfate reducing. One way to think about 

this 40 mg/L threshold is that it represents the reduced aqueous concentration achieved when all available Fe3+ on the sediments has been completely transferred to the aqueous phase. 

However, this assumption is wrong because there is no direct 

relation between the Fe3+ available in the glacial material and 

the resulting Fe2+ that could be produced by redox reactions. 

In other words, it is incorrect to assume that the dissolved 

Fe2+ measured in the aquifer pore water has a valid quantitative relation to the assimilative iron-reducing capacity of the 

sediments. For example, at the Bemidji, Minnesota site, if all 

the bioavailable sediment Fe3+ were converted to Fe2+, then the 

resulting Fe2+ concentration would reach about 5,700 mg/L. 

However, the observed concentration at Bemidji is only 

50 mg/L (Essaid and others, 1995; Bekins and others, 2005). 

The discrepancy is probably due to removal of aqueous Fe2+

reduced from Fe3+ through precipitation as sulfides, carbonates, or mixed-iron-phase oxides (Bekins and others, 2001). 

The simple geochemical model adopted in this study cannot 

reproduce the transfer of iron from Fe3+ to Fe2+ and the subsequent removal of some of the Fe2+ by other chemical reactions. 

Instead, it simply assigns a maximum Fe2+ concentration of 

40 mg/L that is assumed to be the maximum concentration that 

would arise in the absence of any sorption and precipitation 

of the Fe2+ derived from solid-phase Fe3+. Work at contaminated sites (Bekins and others, 2001; Chapelle and others, 

2002) suggests that the 40 mg/L value could underestimate 

the concentration of dissolved Fe2+ concentration that would 

arise from the reservoir of bioavailable iron on sediments 

and would then be subject to sorption and precipitation. The 

implications of underestimating the bioavailable iron from the 

standpoint of the modeling are to shorten the iron-reducing 

phase in the redox reaction chain and to minimize the importance of iron in oxidizing organic carbon relative to the other 

strongly reduced phases. For the purposes of this study, it is 

effectively assumed that naturally occurring bioavailable iron 

is quite low in the glacial sediments from the start, thereby 

allowing the redox reactions to continue beyond iron reduction 

to sulfate-reducing and methanogenic phases without exhausting available organic carbon. 

The maximum allowable CH4 concentration was set to 

10 mg/L. This value is a likely maximum because it is roughly 

equal to the aqueous solubility of methane in shallow ground 

water. This concentration of CH4 would be simulated only 

under extremely reducing conditions.

The transport model was solved by means of implicit 

finite-difference methods. Routine application of an alternative solver scheme that would minimize numerical dispersion 

was impractical because of the duration of the simulations, but 

alternative schemes were applied in the sensitivity phase.

Definition of Redox Zones

In the context of the model simulations, redox conditions 

can be presented in two ways: in terms of the active redox 

reaction in each model cell or in terms of the concentrations of 

redox reactants or products in each cell. The two approaches 

differ somewhat because the concentration of redox reactants 

or products in each cell depends not only on the active redox 

reaction but also on the chemistry of the water that is transported into the cell by advective and dispersive flow. As a 

result, concentrations of redox reactants or products can overlap in ways that depart from the strict sequencing of reactions, 

a condition commonly detected in field ground-water sampling 

(for example, as described by Schreiber and others, 2004). 

For this reason, we present the model results in terms of the 

concentration of redox reactants or products, using a method 

based on Paschke (2007). Six redox water-chemistry zones are 

listed in table 4. The first and second zones are dominated by 

Table 4. Definition of redox water-chemistry zones in terms of dissolved species.

[mg/L, milligrams per liter; ≥, greater than or equal to; <, less than; strongly reducing zones are shaded]

Water-chemistry zone Speciation

Oxic O2 ≥ 0.5 mg/L; Fe2+ < 0.1 mg/L

Nitrate reducing O2 < 0.5 mg/L; NO3 ≥ 0.1 mg/L; Fe2+ < 0.1 mg/L

Mixed: Oxic with Fe2+ O2 ≥ 0.5 mg/L; NO3 < 0.1 mg/L; Fe2+ ≥ 0.1 mg/L

Mixed: Nitrate reducing with Fe2+ O2 < 0.5 mg/L; NO3 ≥ 0.1 mg/L; Fe2+ ≥ 0.1 mg/L

Iron reducing with high sulfate O2 < 0.5 mg/L; NO3 < 0.1 mg/L; Fe2+ ≥ 0.1 mg/L; SO4 ≥ 4 mg/L

Iron reducing with low sulfate O2 < 0.5 mg/L; NO3 < 0.1 mg/L; Fe2+ ≥ 0.1 mg/L; SO4 < 4 mg/L
10 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

the dissolved electron acceptors oxygen and nitrate, respectively. The third and fourth represent the overlap of some dissolved Fe2+ with oxygen or nitrate. The last two zones (negligible dissolved oxygen and nitrate, dissolved Fe2+ greater than 

0.1 mg/L) are collectively referred to as “strongly reducing.” 

Hypothetical Scenarios—A Single 

Ground-Water Flow Field with Varied 

Distributions of Organic Carbon 

The MODFLOW-MT3D-RT3D suite of codes combines 

a flow field with controls on transport that allow the simulation of sequential redox reactions for particular configurations 

of electron donors and electron acceptors. Numerical experiments were conducted for a single flow system using three 

different hypothetical distributions of elevated organic carbon 

concentrations (source areas). The flow system and source 

areas are described below. 

Ground-Water Flow System

In general, ground-water flow fields are complex in the 

sense that waters of different ages and origins mix under the 

influence of different scales of subsurface heterogeneity. This 

mixing typically yields waters of varying chemistry within a 

single aquifer or confining unit.

In the Susquehanna River model, pronounced heterogeneity in hydraulic conductivities that is associated with smallscale grain-size differences in glacial deposits gives rise to a 

complicated velocity and age-distribution pattern. The velocity 

is also influenced by the effective porosity of the sediments, 

assumed to be 0.2 in aquifer sediments and 0.02 in lacustrine 

sediments. Velocities in the aquifers range over several orders 

of magnitude (fig. 5A, 5B), and the rate of vertical flow also 

is highly variable in magnitude and direction (fig. 5C). Flow 

between the aquifers is mostly downward except along the 

course of the Susquehanna River. 

The transport model considers not only movement due to 

the advective velocity in each cell but also the movement due 

to assumed heterogeneity at a scale smaller than the grid size, 

the so-called macrodispersion term. Because the cell-by-cell 

heterogeneity is so large in this model, the assumed dispersion 

terms are relatively small, equal to 0.3 m/d in the longitudinal 

direction, 0.03 m/d in the transverse horizontal direction, and 

0.003 m/d in the transverse vertical direction. 

The ages of the waters in the model were calculated by 

use of the MT3D transport model to take account of movement both by advection and dispersion (following the method 

of Goode, 1996). The average age of ground water in the 

upper aquifer is 3.2 years, with ground water in 95 percent of 

the cells between 0.2 and 10 years. The average age of ground 

water in the lacustrine confining unit is 33.4 years, with 

ground water in 95 percent of the cells between 2 and 132 

years. The average age of ground water in the lower aquifer 

is 10.4 years, with ground water in 95 percent of the cells 

between 1 and 31 years. The spatial distribution of groundwater age for the lower aquifer is shown in figure 6. 

If ground water less than 10 years old is assumed to be 

young, then these results indicate that the water in the upper 

aquifer is young, the water in the lower aquifer is young or not 

very old, and the water in the confining unit is mostly old.

Organic Carbon Source Areas 

In this part of the Susquehanna Valley, the most likely 

sources of elevated organic carbon are landfills (Allan D. 

Randall, U.S. Geological Survey (retired), written commun., 

2006). However, for this study, the distribution of elevated 

organic carbon concentrations (source areas) was hypothetical. For each hypothetical source area, an unlimited reservoir 

of organic carbon was assumed to exist on the solid phase in 

contact with the ground water. The dissolved organic carbon 

concentration in contact with this solid phase over a model cell 

was assumed to be a constant value between 10 and 100 mg/L, 

which is at the high end of dissolved organic carbon concentrations detected in the glacial aquifer system (table 2). The 

source areas were configured to represent naturally occurring 

deposits of elevated organic carbon that might be present in a 

glacial aquifer system. 

The three scenarios are not mutually exclusive and do 

not represent competing hypotheses to explain the redox 

distribution observed at an actual study site. Rather, they are 

presented to test a method for simulating redox under a range 

of conditions. By considering multiple scenarios, it is possible 

to demonstrate how different organic carbon distributions 

interact with the ground-water flow field to produce different 

patterns of reduced water. The three hypothetical source areas 

are described below. 
Hypothetical Scenarios—A Single Ground-Water Flow Field with Varied Distributions of Organic Carbon 11

Figure 5. Model-computed flow velocities in upper and lower aquifers. A, Magnitude of horizontal 

seepage velocity in upper aquifer (model layer 1). B, Magnitude of horizontal seepage velocity in lower 

aquifer (model layer 3). C, Downward flow into lower aquifer (model layer 3). Note: white areas are zones of 

upward flow from lower aquifer.

Horizontal velocity,

in meters per day

1.000e-006

1.000e-004

1.000e-002

1.000

100.0

0 2,000 METERS

Magnitude of horizontal 

seepage velocity in upper 

aquifer (model layer 1).

Magnitude of horizontal 

seepage velocity in lower 

aquifer (model layer 3).

Downward flow into lower 

aquifer (model layer 3). Note: 

white areas are zones of 

upward flow from lower 

aquifer.

B

C

A

Horizontal velocity,

in meters per day

1.000e-006

1.000e-004

1.000e-002

1.000

100.0

0 2,000 METERS

1.000e-008

1.000e-006

1.000e-004

1.000e-002

Downward Darcy velocity, 1.000

in meters per day

0 2,000 METERS
12 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Random Source Area

For this scenario, dissolved organic carbon concentrations 

of 20 mg/L were added at random locations in the upper aquifer (model layer 1), occupying 3.8 percent of the total saturated volume of the upper aquifer (fig. 7A). This source area is 

intended to represent randomly distributed pockets of organic 

carbon associated with buried wetlands or kettle-hole infilling. 

Lacustrine Source Area 

For this scenario, fine-grained material in the confining 

unit (model layer 2) was assumed to have a dissolved organic 

carbon concentration of 10 mg/L. The fine-grained material 

amounted to 78 percent of the volume of the layer (fig. 7B). It 

is likely that the glacial lacustrine facies identified by Randall (1986) was deposited under permafrost conditions and is 

not a host of organic material (Jeffrey Starn, U.S. Geological 

Survey, written commun., 2006). This scenario is included 

nonetheless to investigate how a clay-rich confining unit with 

elevated concentrations of organic carbon would affect redox 

conditions in the adjacent aquifer. 

Wetland Source Area

For this scenario, dissolved organic carbon concentrations of 100 mg/L were added to the upper aquifer (model 

layer 1) at four locations, amounting to 4.6 percent of the 

saturated volume (fig. 7C). This scenario was intended to represent active wetlands or marshland with very high concentrations of organic carbon at land surface. Of the four wetland 

source areas, two are along the rivers in areas of upward 

ground-water flow, and two are in pitted depressions in areas 

of downward ground-water flow. 

The flow system transports dissolved organic carbon 

from the source areas; redox reactions occur at any model-cell 

location as long as the concentrations of the organic carbon 

and at least one electron acceptor are greater than 0.01 mg/L. 

The redox simulations were allowed to proceed in time until 

stable conditions were achieved everywhere in the three-layer 

system for both the electron donor and the five electron acceptors. In all cases, stable conditions were achieved in less than 

500 years. This result suggests that the flow system and the 

geochemical reactions interact at a timescale that is short relative to the age of the glacial deposits.

Results of Hypothetical Scenarios

Equilibrium Organic Carbon Concentrations

The interaction of the ground-water flow field, the presence of localized sources of organic carbon, and the sequence 

of redox reactions yield an equilibrium distribution of dissolved organic carbon in the ground water that is very different from the initial uniform concentration of 1 mg/L. The 

equilibrium concentration of dissolved organic carbon at any 

model cell is between 0 and the concentration at source areas, 

which is fixed to 20 mg/L in the random scenario, 10 mg/L in 

the lacustrine scenario, and 100 mg/L in the wetland scenario. 

Figure 8 shows the spatial distribution of the dissolved organic 

carbon that results at equilibrium in the lower aquifer for the 

three simulated scenarios.

Figure 6. Spatial distribution of ground-water age in the lower aquifer.

Age, in years

0.1

1.0

10.0

100.0

0 2,000 METERS
Results of Hypothetical Scenarios 13

Figure 7. Electron-donor source areas and simulated electron donor concentrations in source layer for three organic carbon 

source scenarios. A, Random organic carbon scenario (Csource = 20 mg/L). Electron donor source in upper aquifer. B, Lacustrine 

organic carbon scenario (Csource = 10 mg/L). Electron donor source in aquitard. C, Wetland organic carbon scenario (Csource = 

100 mg/L). Electron donor source in upper aquifer. [Csource, source concentration; mg/L, milligrams per liter]

0 2,000 METERS

Concentration mg/L,

log scale

2.000e-003

2.000e-002

0.200

2.000

20.00

Upward

Upward

Downward

Downward

flow flow

flow

flow

1.000e-003

1.000e-002

0.100

1.000

10.00

Concentration mg/L,

log scale

1.000e-002

0.100

1.000

10.00

100.0

Concentration mg/L,

log scale

C Wetland organic carbon scenario (Csource = 100 mg/L)

Purple areas show electron donor source in upper aquifer. Upper aquifer absent where white. 

B Lacustrine organic carbon scenario (Csource = 10 mg/L)

Purple areas show electron donor source in aquitard. Aquitard absent where white.

A Random organic carbon scenario (Csource = 20 mg/L)

Purple areas show electron donor source in upper aquifer. Upper aquifer absent where white.
14 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Figure 8. Simulated dissolved organic carbon concentrations in the lower aquifer for three source scenarios. A, Random 

organic carbon scenario (Csource = 20 mg/L). B, Lacustrine organic carbon scenario (Csource = 10 mg/L). C, Wetland organic carbon 

scenario (Csource = 100 mg/L). [Csource, source concentration; mg/L, milligrams per liter]

0 2,000 METERS

A

B

C

Random organic carbon scenario (Csource = 20 mg/L)

Lacustrine organic carbon scenario (Csource = 10 mg/L)

Wetland organic carbon scenario (Csource = 100 mg/L)

2.000e-003

2.000e-002

0.200

2.000

20.00

1.000e-003

1.000e-002

0.100

1.000

10.00

1.000e-002

0.100

1.000

10.00

100.0

Concentration mg/L,

log scale

Concentration mg/L,

log scale

Concentration mg/L,

log scale
Results of Hypothetical Scenarios 15

Equilibrium Redox Water-Chemistry Conditions

The output of the three hypothetical scenarios is presented in terms of redox water-chemistry zones in the lower 

aquifer (fig. 9; table 5).

For the random scenario, elevated organic carbon 

(20 mg/L) is randomly distributed within the upper aquifer to 

simulate buried former wetlands. At equilibrium, the redox 

water chemistry in the lower aquifer is predominantly oxic 

(55 percent) and mixed (39 percent), with isolated areas that 

are strongly reduced (6 percent) (fig. 9A; table 5). The ratio of 

strongly reduced water in the upper aquifer to that in the lower 

aquifer is close to 1, indicating similar volumes of reduced 

water in the upper and lower aquifer—although with very different spatial patterns. 

For the lacustrine scenario, elevated organic carbon 

(10 mg/L) is associated with fine-grained material within the 

confining unit to simulate an old, buried source of dispersed 

organic carbon that is less reactive than an active or buried 

wetland. At equilibrium, the redox water chemistry in the 

lower aquifer is predominantly mixed (69 percent), with lesser 

amounts of oxic water (20 percent) and strongly reduced water 

(11 percent) (fig. 9B; table 5). For this scenario, about 7 times 

more strongly reduced water is in the confining unit than in the 

lower aquifer. 

For the wetland scenario, four source areas were located 

at the water table to simulate active wetlands with very high 

concentrations of organic carbon (100 mg/L). Two of the 

wetland source areas are in areas of upward ground-water flow 

(river wetlands), and the other two are in areas of downward 

flow (pitted wetlands). At equilibrium, the redox water chemistry in the lower aquifer is predominantly oxic (93 percent), 

with a small, strongly reducing zone beneath the pitted wetlands (fig. 9C; table 5). For this scenario, about twice as much 

strongly reduced water is in the upper aquifer than in the lower 

aquifer. 

The loading rates of oxygen and sulfate are constant 

across the three scenarios, but the loading rates of the organic 

carbon differ markedly, as do the rates at which Fe2+ and methane are produced (table 6). 

Comparison of mass-flux data for the three scenarios 

shows that rates of organic carbon loading do not correspond 

to rates of Fe2+ and methane production (table 6). In particular, 

although the random scenario sustains the highest overall rate 

of organic carbon loading to the ground water, the wetland 

scenario yields by far the largest mass of reduced products. 

The random and wetland source areas cover about the same 

fraction of the upper aquifer, but half of the wetland source 

areas are in ground-water discharge areas and therefore do not 

contribute to organic carbon loading of the lower aquifer. In 

contrast, for the random-source scenario, the great majority 

of elevated organic carbon is in recharge areas and therefore 

increases organic carbon loading to the lower aquifer. Despite 

the generally more favorable location of the random sources 

relative to the wetland sources with respect to the spread of 

reduced water, the wetland scenario yields a larger overall 

mass of reduced products as a result of higher source intensity. 

The large contribution of reduced products in the wetland 

scenario is associated with the very strong reducing conditions 

that are generated around and under the wetlands, which in 

turn is associated with the elevated concentrations (100 mg/L) 

of dissolved organic carbon hypothesized to occur. These 

elevated concentrations promote sequential redox reactions 

that go beyond the aerobic and nitrate-reducing stage to ironreducing and even methanogenic conditions, as indicated by 

the iron-reducing, low-sulfate water-chemistry zone in fig. 9C.

It is also noteworthy that the lacustrine scenario yields 

the smallest overall mass of reduced products (table 6), even 

though it has a very large volumetric extent and yields a 

large volume of strongly reduced ground water (77 percent 

of the confining unit; table 5). This discrepancy is related to 

the location of the lacustrine confining unit with respect to 

the ground-water flow field. Although the lacustrine material 

separates the upper from the lower aquifer over more than 

three-quarters of the area, about 95 percent of the exchange 

between the two aquifers occurs where they are in direct connection. Consequently, the amount of water passing through 

the source areas is small, and the rate of transfer of organic 

carbon from the solid to dissolved phase is limited. Moreover, 

the reduced water leaking out of the lacustrine source areas 

to the lower aquifer is strongly diluted by oxic water moving 

directly downward from the upper aquifer. However, this dilution effect is not sufficient to prevent reducing waters in one 

large area of the lower aquifer (fig. 9B).

A comparison of a map of simulated redox reaction zones 

(defined by the RT3D equations) to the simulated water chemistry (based on the classification in table 4) demonstrates one 

difficulty in relating observed reactant concentrations to the 

reducing conditions that produce them. For the purposes of the 

comparison, a special simulation based on the random source 

scenario was done with the cutoff for switching from aerobic 

reactions to more reducing reactions set to 0.5 mg/L instead 

of 0.01 mg/L. In this way, the cutoff for switching reactions 

is the same as the cutoff for distinguishing oxic redox water 

chemistry from other redox water chemistry. The modelgenerated redox reaction zones, which cannot be observed 

directly in the field, show that almost all the lower aquifer 

has sufficient oxygen to support aerobic reactions (fig. 10A). 

The accompanying map of model-generated water-chemistry 

zones, which can, in principle, be observed in the field, show 

that most of the lower aquifer falls into a mixed zone where 

dissolved oxygen is greater than 0.5 mg/L and Fe2+ is greater 

than 0.1 mg/L (fig. 10B). The overlap of dissolved iron with 

dissolved oxygen can be attributed to the spreading by advection and (or) dispersion of iron-enriched water (from reduced 

zones influenced by organic carbon source areas in the upper 

aquifer) with oxic water (which did not pass through those 

source areas and which contains a background Fe2+ concentration of 0.005 mg/L). 
16 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Figure 9. Simulated redox water chemistry in the lower aquifer (layer 3) for three organic carbon source scenarios. Cells 

shown in white represent areas where the lower aquifer is pinched out. A, Random organic carbon scenario (Csource = 20 mg/L). 

B, Lacustrine organic carbon scenario (Csource = 10 mg/L). C, Wetland organic carbon scenario (Csource = 100 mg/L). [Csource, source 

concentration; mg/L, milligrams per liter]

EXPLANATION

Simulated water-chemistry zones

Oxic

Nitrate reducing

Iron reducing with low sulfate

Iron reducing with high sulfate

Mixed: oxic with Fe2+

Mixed: nitrate reducing with Fe2+

0 2,000 METERS

A

B

C

Random organic carbon scenario (Csource = 20 mg/L)

Lacustrine organic carbon scenario (Csource = 10 mg/L)

Wetland organic carbon scenario (Csource = 100 mg/L)
Results of Hypothetical Scenarios 17

Table 5. Simulated redox water-chemistry zones in the upper aquifer, confining unit, and lower aquifer by percent volume.

[Csource, source concentration; mg/L, milligrams per liter; strongly reducing categories are shaded]

Water-chemistry zones

Percent volume

Upper aquifer Confining unit Lower aquifer

Random organic carbon scenario (Csource= 20 mg/L)

Source formation is upper aquifer.

Ratio of volume of strongly reduced water, upper aquifer to lower aquifer = 1.14

Oxic 61.7 59.2 54.9

Nitrate reducing 0.0 0.0 0.2

Mixed: Oxic with Fe2+ 31.9 25.3 39.2

Mixed: Nitrate reducing with Fe2+ 0.0 0.1 0.1

Iron reducing with high sulfate 6.0 12.9 5.6

Iron reducing with low sulfate 0.4 2.6 0.0

TOTAL 100.0 100.1 100.0

Lacustrine organic carbon scenario (Csource= 10 mg/L)

Source formation is confining unit.

Ratio of volume of strongly reduced water, upper aquifer to lower aquifer = 7.25

Oxic 38.9 20.9 20.4

Nitrate reducing 0.0 0.0 0.0

Mixed: Oxic with Fe2+ 61.0 2.2 69.0

Mixed: Nitrate reducing with Fe2+ 0.0 0.0 0.0

Iron reducing with high sulfate 0.1 6.6 9.7

Iron reducing with low sulfate 0.0 70.3 0.9

TOTAL 100.0 100.0 100.0

Wetland organic carbon scenario (Csource = 100 mg/L)

Source formation is upper aquifer.

Ratio of volume of strongly reduced water, upper aquifer to lower aquifer = 1.95

Oxic 88.7 94.0 92.7

Nitrate reducing 0.0 0.0 0.0

Mixed: Oxic with Fe2+ 3.2 1.5 3.2

Mixed: Nitrate reducing with Fe2+ 0.0 0.0 0.0

Iron reducing with high sulfate 1.3 0.8 0.8

Iron reducing with low sulfate 6.7 3.7 3.3

TOTAL 99.9 100.0 100.0

Table 6. Mass flux of redox constituents for organic carbon source scenarios for entire model.

[All data are in kilograms per day; O2, SO4, and organic carbon loaded to system are either consumed by redox reactions or discharged to surface water.]

Redox constituent Recharge and river, 

all scenarios

Organic carbon source scenario

Random Lacustrine Wetland

O2 loading 729 0 0 0

SO4 loading 1,193 0 0 0

Organic carbon loading 117 267 48 193

Fe2+ produced by reduction of Fe3+ 0 157 133 422

CH4 produced by reduction of CO2 0 0.2 8 40
18 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Figure 10. Comparison of A, simulated redox reaction zones to B, simulated redox water-chemistry zones for random organic 

carbon scenario, lower aquifer (layer 3). Cells shown in white represent areas where the lower aquifer is pinched out. The lower 

aquifer is represented by layer 3 in the model.

0 2,000 METERS

0 2,000 METERS

EXPLANATION

Inactive

Simulated redox reaction zones

Fe2+

CH4

DO2

SO4

NO3

EXPLANATION

Simulated water-chemistry zones

A

B

Oxic

Nitrate reducing

Iron reducing with low sulfate

Iron reducing with high sulfate

Mixed: oxic with Fe2+

Mixed: nitrate reducing with Fe2+
Sensitivity Analysis 19

The variability of the dissolved Fe2+ concentration in the 

lower aquifer for this simulation (fig. 11) shows the extent of 

the mixing that occurs as a result of the interaction of the flow 

system with the electron-donor source areas. Although almost 

all of the lower aquifer in this hypothetical simulation is oxic, 

there are broad areas of elevated dissolved Fe2+.

Sensitivity Analysis

The methodology described here for simulating redox 

at the field scale depends on inputs that control the solver 

strategy, rates of reaction, the availability of terminal electron 

acceptors, and the reactivity of electron donors. An analysis 

was performed to evaluate the sensitivity of the results to 

each category of input using a version of the random-source 

scenario in which the source concentration cells are set to 

constant values of 100 mg/L instead of 20 mg/L. The high 

concentration of electron donors results in large areas across 

the lower aquifer with strongly reduced water for the base run 

(fig. 12). Each sensitivity run is a variation of the base run 

with one input changed. The change in outcome relative to the 

base run was quantified in terms of the redox water-chemistry 

zones (table 4)—specifically, the percentage of the loweraquifer volume classified as strongly reducing and the ratio of 

the strongly reduced volume in the upper aquifer to that in the 

lower aquifer.

In the base sensitivity simulation, the volume of reduced 

water in the lower aquifer is 29.6 percent, and the ratio of the 

reduced volume in the upper to lower aquifer is 0.88, implying 

slightly more reduced water in the lower aquifer than in the 

upper aquifer.

Figure 11. Simulated concentration of dissolved Fe2+ for random organic carbon scenario, lower aquifer (layer 3). Cells shown in 

white represent areas where the concentration is below 4.000e-002.

4.000e-002

0.400

4.000

40.00

0 2,000 METERS

Concentration mg/L,

log scale

EXPLANATION
20 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Sensitivity to Solver Strategy

The following was performed to test the robustness of 

the numerical solution: (1) the vertical discretization was 

increased within the confining unit, (2) the finite-difference 

solution, which is subject to numerical dispersion, was 

replaced with the TVD (total-variation-diminishing) solver, 

which requires longer runtimes but is less subject to numerical 

dispersion (Chunmiao and Wang, 1998), and (3) the transport 

time step was decreased from 5 days to 0.5 day. These input 

changes yield very small changes with respect to the reduced 

volume in the lower aquifer or the ratio of reduced volume 

in the upper to the lower aquifer (table 7). They imply that 

numerical dispersion does not strongly influence the solution.

Sensitivity to Attenuation Parameters

Attenuation parameters quantify processes that cause the 

concentration of a redox reactant to decrease. They include the 

assumed decay rate for the first-order reactions, the consideration of retardation, and the degree of assumed dispersion 

mostly to account for heterogeneities at a scale smaller than 

the grid size (table 7). The sensitivity runs show that the solution is insensitive to the first-order decay constant assumed in 

the organic carbon source areas, but it is moderately sensitive to an overall increase of the decay rate from 0.001 to 

0.01 day-1 and strongly sensitive to a decrease from 0.001 to 

0.0001 day-1. In the latter case, the slow rate of reaction means 

that although the equilibrium volume of reduced water in the 

lower aquifer is small (8.5 percent), it is 4 times the volume in 

the upper aquifer. 

Retardation of the electron donor has negligible effect 

on the equilibrium solution, whereas the choice of dispersivity parameters has moderate influence within the range tested 

(that is, an increase of 10 times causes a small increase of 

reduced water in the lower aquifer). 

Sensitivity to Concentrations of Electron 

Acceptors

The background water quality was altered to test the sensitivity of the model to changes in concentrations of electron 

acceptors (table 7). For the alternate background water quality, 

concentrations were decreased for dissolved O2 (from 6 to 

4 mg/L), increased for NO3 (from 0.2 to 9 mg/L), increased for 

Fe2+ (from 0.005 to 0.01 mg/L), and increased for SO4 (from 

10 to 100 mg/L). The alternate background water quality was 

judged to be representative of high-nitrate and high-sulfate 

waters from the glacial aquifer system in parts of Nebraska 

(Thomas, 2007). Results show that increasing the concentrations of electron acceptors in the recharge and river water 

yielded a lower volume of strongly reduced water in the lower 

Figure 12. Random organic carbon scenario for sensitivity analysis (Csource = 100 mg/L); simulated redox water-chemistry 

conditions in lower aquifer (layer 3). Cells shown in white represent areas where the lower aquifer is pinched out.

EXPLANATION

Observed water-chemistry zones

0 2,000 METERS

Oxic

Nitrate reducing

Iron reducing with low sulfate

Iron reducing with high sulfate

Mixed: oxic with Fe2+

Mixed: nitrate reducing with Fe2+
Sensitivity Analysis 21

Table 7. Results of sensitivity simulations.

[mg/L, milligrams per liter; --, not applicable]

Simulation

Volume of strongly reduced water Sensitivity

Lower aquifer

(percent)

Ratio of upper aquifer

to lower aquifer

Sensitivity to solver strategy

Scenario: Random donor source concentration in upper aquifer = 20 mg/L 5.6 1.14 --

Base run: Source concentration changed to 100 mg/L 29.6 0.88 --

1. Lacustrine confining unit represented as 3 layers instead of 1 29.0 0.89 Weak

2. TVD solver instead of finite difference solver 29.2 0.83 Weak

3. Transport time step = 0.5 day instead of 5 days 29.5 0.88 Weak

Sensitivity to attenuation parameters

Scenario: Random donor source concentration in upper aquifer = 20 mg/L 5.6 1.14 --

Base run: Source concentration changed to 100 mg/L 29.6 0.88 --

1. First-order rate constant in source areas = 0.1 instead of 0.001 day-1 29.9 0.89 Weak

2. First-order rate constant everywhere = 0.01 instead of 0.001 day-1 25.9 1.01 Moderate

3. First-order rate constant everywhere = 0.0001 instead of 0.001 day-1 8.5 0.27 Strong

4. Donor retardation = 10 instead of 1 outside of source areas 29.6 0.88 Negligible

5. Increased dispersivity by 10x: longitudinal = 10 feet instead of 1 foot 32.0 0.84 Moderate

Sensitivity to concentrations of electron acceptors

Scenario: Random donor source concentration in upper aquifer = 20 mg/L 5.6 1.14 --

Base run: Source concentration changed to 100 mg/L 29.6 0.88 --

1. Alternate background water quality1 17.5 0.99 Strong

2. Nitrate loading: concentration in recharge set to 10 mg/L for 20 years 15.7 1.00 Strong

3. Maximum allowable dissolved Fe2+ concentration 10 mg/L instead of 

40 mg/L

29.7 0.89 Weak

Sensitivity to electron-donor behavior

Scenario: Random donor source concentration in upper aquifer = 20 mg/L 5.6 1.14 --

Base run: Source concentration changed to 100 mg/L 29.6 0.88 --

1. No electron-donor source areas, only in recharge and river water 0.0 -- Strong

2. Redox reactions only occur in electron donor source areas 

(fixed instead of mobile source)

0.1 37 Strong

1

 Recharge and river concentrations changed: O2 from 6 to 4 mg/L, NO3 from 0.2 to 9 mg/L, Fe2+ from 0.005 to 0.01 mg/L, and SO4 from 10 to 

100 mg/L.
22 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

aquifer at equilibrium conditions. Specifically, the relatively 

high concentration of dissolved nitrate causes more of the 

available organic carbon to be exhausted when going from aerobic to nitrate conditions and, therefore, redox reactions cease 

in many areas before attaining strongly reduced conditions. 

The influence of nitrate loading (due to fertilizer applications, for example) was evaluated by extending the equilibrium base solution for 20 years with 10 mg/L of nitrate 

added to recharge. The effect is similar to the case above: the 

sequential reactions in much of the lower aquifer never pass 

beyond the nitrate-reducing phase, so there is less strongly 

reduced water (table 7). 

To test the sensitivity of the solution to the assumption 

that Fe2+ produced by the reduction of solid-phase Fe3+ cannot exceed an assumed amount, the maximum dissolved Fe2+

threshold was reduced from 40 mg/L to 10 mg/L. The effect of 

this change on the overall amount of strongly reduced water 

in the lower aquifer is weak (table 7), but the distribution of 

reduced water shifted from plentiful iron-reduced water in the 

base case to more sulfate-reduced water in the sensitivity run.

Sensitivity to Electron-Donor Behavior

A sensitivity run was done to test the effect of limiting 

the availability of electron donors to the 1 mg/L dissolved 

organic carbon concentration assumed in the recharge water 

and river water (that is, the source areas of organic carbon 

were eliminated). Under this assumption, the oxygen supplied 

by the inflow to the upper aquifer overwhelms the low availability of electron donors to the point that there is no strongly 

reduced water in the lower aquifer (table 7).

The proposed methodology includes the key assumption 

that redox reactions between electron donors and electron 

acceptors continue to occur in the dissolved phase downgradient from the fixed electron donor source areas. In other words, 

it is assumed that a mobile dissolved source of electron donor 

carried downgradient by the ground water causes redox reactions to spread, although they tend to occur at a deceasing rate 

as the concentration of the dissolved electron donor declines. 

This conceptual framework is not necessarily consistent with 

typical behavior in the presence of localized anthropogenic 

electron donor sources, where almost all the reducing activity 

appears to occur in close proximity to the fixed, solid-phase 

source areas. It is not clear whether this model associated 

with contaminated sites is applicable to natural systems. But 

if, for testing purposes, the assumption is made that the locus 

of reaction in natural systems is limited to the solid-phase 

organic carbon, thereby eliminating downgradient reactions in 

the dissolved phase, then the volume of strongly reduced water 

in the lower aquifer drops to almost zero (table 7). 

By calculating the difference in the rates of reaction 

between the base run with both fixed and mobile sources of 

reaction and the sensitivity run with only fixed sources, it is 

possible to determine what fraction of the reducing activity 

occurs in model cells defined as electron-donor source areas 

and what fraction occurs in downgradient waters. The rates 

of consumption of dissolved oxygen and sulfate (table 8) 

and the rates of production of Fe2+ and methane (table 8) are 

each about 4 times those in the base case than in the case of 

no downgradient reactions. Consequently, for the hypothetical geochemical conditions of the random-source case, about 

three-quarters of the activity occurs downgradient from the 

Table 8a. Simulated degradation rates of redox reactants for mobile versus fixed electron donor source.

[All data are in kilograms per year; Csource, source concentration; mg/L, milligram per liter]

Simulated degradation rates of redox reactants (base sensitivity run: random organic carbon source, Csource= 100 mg/L)

Location of redox reaction Organic carbon 

degradation

Dissolved oxygen 

degradation

Sulfate 

degradation

Redox reactions occur everywhere (fixed and mobile donor source) 715 385 252

Redox reactions occur in source areas only (fixed donor source) 0.0 86 68

Table 8b. Simulated production rates of redox products for mobile versus fixed electron donor source.

[All data are in kilograms per year; Csource, source concentration; mg/L, milligram per liter]

Simulated production rates of redox products (base sensitivity run: random donor source, Csource= 100 mg/L)

Location of redox reaction Ferrous iron 

production

Methane

production

Redox reactions occur everywhere (fixed and mobile donor source) 1,411 64

Redox reactions occur in source areas only (fixed donor source) 414 16
Influence of Pumping on Distribution of Reduced Waters 23

Figure 13. Simulated drawdown contours (meters) around pumping wells.

Layer 3 

1

1

1

2

1

2

2

1

1

2

2

1

4

4

EXPLANATION

Well location in lower aquifer

0 2,000 METERS

source areas. Further research would be needed to determine 

whether this is a realistic ratio or whether a greater proportion 

of the redox activity occurs in the fixed source areas. 

Influence of Pumping on Distribution of 

Reduced Waters

The lower aquifer (model layer 3) in the bedrock valley 

along the Susquehanna River is important for water supply. 

Records show that in the 1980s, 21 municipal water-supply 

wells collectively discharged about 21 Mgal/d from the 

lower aquifer (Walcott and Coon, 2001). Model simulations 

of pumping produce a steady-state flow field with cones of 

depression around the wells (fig. 13).

According to the flow model, pumping discharges about 

55 percent of inflow to the flow system from recharge and 

river loss. Pumping causes increased downward leakage from 

the upper to lower aquifer and an increased proportion of 

leakage (from 95 percent to 98 percent) in areas where the two 

aquifers are in contact.

The transport code was used to simulate the hypothetical geochemical response of the lower aquifer to the new flow 

field after 20 years of pumping for each of the three source 

scenarios. A major effect is to induce more oxygenated water 

from rivers into the subsurface and lower aquifer, thereby 

moderately increasing the extent of the oxic zones and moderately shrinking the size of the strongly reducing zones in the 

lower aquifer.

The effect of pumping on concentrations of redox 

reactants for the three hypothetical scenarios (table 9) indicates that, in general, the increased circulation of water due 

to pumping has a more pronounced effect on the chemistry of 

water entering wells than it does on the average concentration 

of reactants across the entire lower aquifer. For example, 

the decrease in average Fe2+ concentrations is appreciable in 

well discharge for all three scenarios. This reflects the close 

connection of the deep aquifer to the shallow water bodies in 

the Susquehanna River Valley intervening low-permeability 

lacustrine deposits. The aquifers under study constitute an 

open system from the standpoint of redox, and pumping wells 

capture mostly young water.

The differences among scenarios are also noteworthy. For 

the random scenarios, the average dissolved oxygen concentration in the lower aquifer is minimally affected by pumping, but the dissolved oxygen in well water increases by 50 

percent. In contrast, for the lacustrine and wetland scenarios, 

the dissolved oxygen concentrations in the aquifer and well 

water do not change greatly. These differences are at least 

partly attributable to the exact location of wells in relation to 

the organic carbon source areas and to areas of good hydraulic connection between the two aquifers. In the case of the 

random scenario, the wells are able to draw mostly young 

water from well-oxygenated shallow areas that are not in close 

contact with electron-donor source areas. In the case of the 

lacustrine scenario, mixing of older water emanating from the 

large volume of organic-rich fine-grained material likely limits 

the increase of dissolved oxygen at the wells. In the case of 

the wetland scenarios, which locally produce the most reduced 

water, some wells show increased dissolved oxygen; however, 

the average value in all well water decreases slightly because 

some wells underlying wetlands capture very reduced water. 

The wetland scenario also shows, on average, less reduction in 

dissolved Fe2+ in well water than do the other two scenarios. 

These subtle differences indicate the difficulty in predicting 

the effect of pumping on redox conditions in well water without a good knowledge of both the flow field and the distribution of electron donors. This difficulty carries over to predicting the effect of pumping on redox-sensitive constituents.
24 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Limitations of the Study 

The primary objective of this study was to show the 

importance of the ground-water flow field in determining the 

spread and distribution of reduced waters in a glacial valleyfill aquifer. Very different patterns of redox conditions are 

produced by the interaction of distinct organic carbon source 

areas with aspects of the ground-water flow field, such as 

(1) the degree of vertical connectedness between aquifers, and 

(2) local velocity variations induced by contrasts in hydraulic 

conductivity values of the unconsolidated glacial sediments. 

Pumping alters the ground-water flow field and, therefore, 

changes the distribution of strongly reduced water, most notably in the vicinity of a pumping well. 

A chief limitation of this study is that the assumed distribution and concentration of organic carbon is hypothetical, 

although the flow model is calibrated to field observations 

and the transport inputs reflect real-world conditions. Consequently, the model results are useful primarily for advancing 

the general understanding of subsurface redox processes. 

Extension of this hypothetical methodology to other glacial 

settings (for example till and outwash plains) would add to the 

general understanding of the processes, but calibration to field 

studies would be necessary to test the predictive capability of 

the method. 

A further limitation arises from the simplicity of the 

geochemical model adopted. The geochemical framework 

could be made more sophisticated in a variety of ways—for 

example, by substituting Monod for first-order kinetics (Waddill and Widdowson, 2000), by allowing for overlapping redox 

reactions (Schreiber and others, 2004), by allowing for sources 

and sinks of constituents such as sulfate within the subsurface 

(Chapelle and others, 1995) or—perhaps most importantly—

by explicitly accounting for the transfer of iron from the solid 

to dissolved phase (Rolle and others, 2006). However, if the 

methodology was to be applied to a field site where understanding of redox is important, such as an area of natural 

arsenic contamination, it would be advantageous to combine a 

sophisticated flow model with a simple geochemical framework requiring few parameters for calibration. Complexity 

in the geochemical model could be limited by calibrating to 

distinct first-order decay constants for each electron acceptor. A field application of the modeling methodology would 

work best in an area where the rate of flow is fast enough that 

the time to equilibrium of the reactants is relatively short (on 

the order of hundreds of years), resulting in relatively short 

computer runtimes. 

An inherent assumption in the adopted geochemical 

model is that redox reactions occur away from the source areas 

of naturally occurring electron donors, although at a reduced 

rate according to a first-order law proportional to electron 

donor concentration. It is possible that reducing reactions are 

limited chiefly to the proximity of source areas so that the 

ground-water flow system serves only to spread reduced water 

rather than to spread both reduced water and reducing reactions. Application of the modeling methodology to an actual 

ground-water system might contribute to a better understanding of the amount of redox activity that actually occurs downgradient from natural electron-donor sources.

An important assumption particular to this modeling 

effort is that the iron reduction state can be exhausted in the 

ground-water system because the reservoir of bioavailable 

iron in the glacial sediments is quite limited. It is assumed 

that its reduction in the absence of any sorption or precipitation would yield a maximum of only 40 mg/L of dissolved 

Fe2+. The distribution of bioavailable iron in distinct sediment 

types is still poorly understood, although some studied sites 

show levels of bioavailable iron that would support aqueous 

Fe2+ concentrations, net of sorption and precipitation, that are 

much higher than 40 mg/L (Bekins, 2001; Chapelle, 2002). 

Table 9. Effect of pumping on concentrations of redox reactants for the three hypothetical scenarios.

[Concentrations are in milligrams per liter]

Scenario Reactant

Average concentration in well 

boreholes in lower aquifer

Average concentration in 

entire lower aquifer

No pumping 20 years 

of pumping No pumping 20 years of pumping

Random donor source O2 2.8 4.2 4.2 4.1

Fe2+ 4.5 0.7 1.7 1.6

CH4 0.0 0.0 0.0 0.0

Lacustrine donor source O2 4.5 5.0 4.6 4.4

Fe2+ 4.9 1.5 4.5 3.1

CH4 0.2 0.1 0.1 0.1

Wetland donor source O2 3.9 3.7 5.1 4.8

Fe2+ 10.2 7.3 1.7 2.1

CH4 1.2 0.6 0.3 0.2
References Cited 25

In natural settings, it would be particularly useful to obtain a 

suite of bioavailable iron measurements on sediments to better judge the likelihood that, in environments with plentiful 

electron-donor source areas, reducing conditions would stop 

at the iron-reducing phase (indicative of relatively high levels 

of bioavailable iron) or, alternatively, would move beyond it 

to sulfate-reducing and methanogenic conditions (indicative of 

low or possibly very low levels of bioavailable iron). It would 

also be useful to correlate concentrations of bioavailable iron 

in sediments with the age of the associated ground water and 

with naturally occurring contaminants that tend to accumulate 

in particular redox phases—arsenic, for example, which tends 

to dissolve in iron-reducing and methanogenic phases and to 

preciptate in sulfate-reducing phases (Smedley and Kinniburgh, 2002). 

Conclusions

Recent studies have emphasized the connection between 

the dynamics of a ground-water flow field and the mobilization of redox-sensitive contaminants in terms of the relation between flow rates and the kinetics of redox reactions 

(Smedley and Kinniburgh, 2002; Harvey and others, 2006). 

To explore this link, the MODFLOW-MT3D-RT3D suite 

of code was applied to combine a ground-water flow model 

of a complex depositional system with a relatively simple 

geochemical model. The geochemical model posits that (1) 

reducing conditions are present in ground water when, due 

to microbial respiration, redox reactions involving electron 

donors and electron acceptors occur and (2) reactions are 

approximately first-order with respect to the electron donor 

and sequential with respect to electron acceptors. Hypothetical 

distributions of electron-donor sources such as organic carbon 

were simulated in conjunction with a real flow field typical of 

a heterogeneous glacial valley-fill aquifer to demonstrate how 

strongly reducing conditions can arise in areas not immediately adjacent to electron-donor source areas, such as an 

underlying aquifer, because of flow and mixing of water along 

three-dimensional flow paths. The modeling simulates the 

parts of the aquifer where strongly reduced water is present 

because the rate at which oxygenated water is introduced from 

the water table is insufficient to offset the rate of redox reactions fueled by available electron acceptors. By comparing the 

simulated spatial distribution of active redox reactions with 

the simulated spatial distribution of water chemistry, it is possible to identify the role of advection and dispersion in mixing 

waters and in transporting, for example, dissolved iron into 

oxic zones. Discharge from wells increases flow and mixing 

which, where the pumped aquifer is in close connection with 

shallow oxygenated waters, can reverse reduced conditions in 

ground water. 

Acknowledgments

This work was supported by the U.S. Geological Survey 

(USGS) National Water-Quality Assessment Program as part 

of the regional assessment of the glacial aquifer system in 

the northern United States. Review comments from Barbara 

Bekins and J. Jeffrey Starn, both of the USGS, dramatically 

improved the organization, clarity, and technical substance 

of the report. For the availability of the flow model used in 

this study we depended on the generous cooperation of Allan 

Randall (USGS-retired) and Richard Yager of the USGS. 

Substantial assistance was provided by USGS researchers and 

editors including Peter McMahon, Michelle Greenwood, and 

Michael Eberle. We extend special thanks to Roberta Bellini 

who assisted with the modeling on leave from her Ph.D. work 

in Milan.

References Cited

Barbash, J.E., 2004, The geochemistry of pesticides, in Lollar, 

B.S., ed., Environmental Geochemisty: Oxford, United 

Kingdom, Elsevier, v. 9, p. 541–577. 

Bekins, B.A., Warren, Ean, and Godsy, E.M., 1998, A comparison of zero-order, first-order, and Monod biotransformation 

models: Ground Water, v. 36, no. 2, p. 261–268.

Bekins, B.A, Cozzarelli, I.M., Godsy, E.M., Warren, E., 

Essaid, H.I., and Tuccillo, M.E., 2001, Progression of natural attenuation processes at a crude oil spill site; II. Controls 

on spatial distribution of microbial populations: Journal of 

Contaminant Hydrology, v. 53, p. 387–406.

Bekins, B.A., Cozzarelli, I.M., and Curtis, G.P., 2005, A simple method for calculating growth rates of petroleum hydrocarbon plumes: Ground Water, v. 43, no. 6, p. 817–826.

Böhlke, J.K., Wanty, R., Tuttle, M., Delin, G., and Landon 

M., 2002, Denitrification in the recharge area of a transient agricultural nitrate plume in a glacial outwash 

sand aquifer: Water Resources Research, v. 38, no. 7, 

doi:10.1029/2001WR000663, 26 p.

Chapelle, F.H., 1993, Ground-water microbiology and geochemistry: New York, John Wiley and Sons, 424 p.

Chapelle, F.H., McMahon, P.B., Dubrovsky, N.M., Fujii, 

R.F., Oaksford, E.T., and Vroblesky. D.A., 1995, Deducing 

the distribution of terminal electron-accepting processes 

in hydrologically diverse groundwater systems: Water 

Resources Research, v. 31, no. 2, p. 359–371.
26 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Chapelle, F.H., Bradley, P.M., Lovley, D.R., O’Neill, K, and 

Landmeyer, J.E., 2002, Rapid evolution of redox processes 

in a petroleum hydrocarbon-contaminated aquifer: Ground 

Water, v. 40, no. 4, p. 353–360.

Clement, T.P., 2001, A modular computer model for simulating 

reactive multi-species in three-dimensional ground water 

systems: Richland, Wash., Pacific Northwest National Laboratory, PNNL–SA–11720, 59 p. 

Clement, T.P., Sun, Y. Hooker, B.S. and Petersen, J.N., 1997, 

Modeling natural attenuation of contaminants in saturated 

groundwater, in Leeson, A., and Alleman, B.C., eds., In situ 

and on site bioremediation: Columbus, Ohio, Battelle Press, 

v. 1 (4–1), p. 37–42.

Drever, J.I., 2002, The geochemistry of natural waters: Surface 

and groundwater environments (3d ed): Upper Saddle 

River, N.J., Prentice Hall, 436 p.

Essaid, H.I., Bekins, B.A., Godsy, E.M., Warren, Ean, 

Baedecker, M.J., and Cozzarelli, I.M.,1995, Simulation of 

aerobic and anaerobic biodegradation processes at a crude 

oil spill site: Water Resources Research, v. 31, no. 12, 

p. 3309–3327.

Goode, D.J., 1996, Direct simulation of groundwater age: 

Water Resources Research, v. 32, no. 2, p. 289–296.

Guoping Lu, Clement, T.P., Chumiao Zheng, and Wiedemeier, 

T.H., 1999, Natural attenuation of BTEX compounds—

Model development and field-scale application: Ground 

Water, v. 37, no. 5, p. 707–717.

Harvey, C.F., Khandaker, N.A., Yu, Winston, Badruzzaman, 

A.B.M., Ali, M.A., Oates, P.M., Michael, H.A., Neumann, R.B., Beckie, R., Islam, S., and Ahmed, M.F., 2006, 

Groundwater dynamics and arsenic contamination in Bangladesh: Chemical Geology, v. 228, p. 112–136.

Hendry, M.J., and Wassenaar, L.I., 2005, Origin and migration 

of dissolved organic carbon fractions in a clay-rich aquitard—14C and d13C evidence: Water Resources Research, 

v. 41, no. 2, W2021, doi:10.1029/2004WR003157, 10 p.

Keating, E.H., and Bahr, J.M., 1998, Reactive transport 

modeling of redox geochemistry—Approaches to chemical 

disequilibrium and reaction rate estimation at a site in northern Wisconsin: Water Resources Research, v. 34, no. 12, 

p. 3573–3584. 

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference flow model: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A1, 586 p.

McNab, W.W., Jr., and Narasimhan, T.N., 1994, Modeling 

reactive transport of organic compounds in groundwater using a partial redox disequilibrium approach: Water 

Resources Research, v. 30, no. 9, p. 2619–2635. 

Moran, M.J., Zogorski, J.S., and Squillace. P.J., 2007, Chlorinated solvents in groundwater of the United States: Environmental Science & Technology, v. 41, no. 1, p. 74–81.

Paschke, S.S., ed., 2007, Hydrogeologic settings and groundwater flow simulations for regional studies of the transport 

of anthropogenic and natural contaminants to public-supply 

wells—Studies begun in 2001: U.S. Geological Survey 

Professional Paper 1737–A, 288 p.

Randall, A.D., 1986, Aquifer model of the Susquehanna River 

Valley in southwestern Broome County, New York: U.S. 

Geological Survey Water-Resources Investigations Report 

85–4099, 38 p.

Randall, A.D., 2001, Hydrogeologic framework of stratifieddrift aquifers in the glaciated northeastern United States: 

U.S. Geological Survey Professional Paper 1415–B, 179 p.

Rolle, M., Zolla, V., Sethi, R., Di Molfetta, A., and Clement. 

T.P., 2006, Modeling TEAPs and computing redox zonation in contaminated aquifer, in Poeter, Eileen; Chunmiao 

Zheng; Hill, Mary; and Seo, Sophia, MODFLOW and 

More 2006—Managing ground-water systems, conference 

proceedings: Golden, Colo., International Ground Water 

Modeling Center, p. 215–219.

Schreiber, M.E., Carey, G.R., Feinstein, D.T., and Bahr, J.M., 

2004, Mechanisms of electron acceptor utilization—Implications of simulating anaerobic biodegradation: Journal of 

Contaminant Hydrology, v. 73, p. 99–127.

Smedley, P.L., and Kinniburgh, D.G., 2002, A review of the 

source, behaviour and distribution of arsenic in natural 

waters: Applied Geochemistry, v. 17, p. 517–568.

Thomas, M.A, 2007, The association of arsenic with redox 

conditions, depth, and ground-water age in the glacial aquifer system of the Northern United States: U.S. Geological 

Survey Scientific Investigations Report 2007–5036, 26 p.

Waddill, D.W., and Widdowson, M.A., 2000, SEAM3D—A 

numerical model for three-dimensional solute transport 

and sequential electron acceptor-based bioremediation in 

groundwater: Vicksburg, Miss., U.S. Army Corps of Engineers, Engineer Research and Development Center, [110] p.

Warner, K.L., and Arnold, T.L., 2005, Framework for regional 

synthesis of water-quality data for the glacial aquifer system 

in the United States: U.S. Geological Survey Scientific 

Investigations Report 2005–5223, 6 p.

Wiedmeier, T.H., Wilson, J.T., Kampbell, D.H., Miller, R.N., 

and Hansen, J., 1995, Technical protocol for implementing 

intrinsic remediation with long-term monitoring for natural 

attenuation of fuel contamination dissolved in groundwater: San Antonio, Tex., Air Force Center for Environmental 

Excellence, Technology Transfer Division, Brooks AFB, 

v. 1 and 2.
References Cited 27

Welch, A.H., Westjohn, D.B., Helsel, D.R., and Wanty, R.B., 

2000, Arsenic in ground water of the United States—

Occurrence and geochemistry: Ground Water, v. 38, no. 4, 

p. 589–604.

Wolcott, S.W., and Coon, W.F., 2001, Simulation of a val-

ley fill aquifer system to delineate flow paths, contributing 

areas, and traveltime to wellfields in southwestern Broome 

County, New York: U.S. Geological Survey WaterResources Investigations Report 01–4171, 25 p.

Zheng, Chunmiao, and Wang, P.P., 1998, MT3DMS—Docu-

mentation and user’s guide: Prepared for U.S. Army Corps 

of Engineers [variously paged].
28 Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting

Prepared by the Wisconsin Water Science Center:

U.S. Geological Survey

Wisconsin Water Science Center

8505 Research Way

Middleton, WI 53562

tel: (608) 828-9901

fax: (608) 821-3817

This publication is available online at URL:

http://pubs.water.usgs.gov/sir2008-5066

Information regarding the Wisconsin Water Science Center

is available at URL:

http://wi.water.usgs.gov

Feinstein and Thomas—Hypothetical Modeling of Redox Conditions Within a Complex Ground-Water Flow Field in a Glacial Setting—Scientific Investigations Report 2008-5066

Science Center Director

USGS Wisconsin Water Science Center

8505 Research Way

Middleton, Wisconsin 53562

http://wi.water.usgs.gov