Model Report for Christina River Basin,
Pennsylvania, Delaware, and Maryland,
High-Flow Nutrient and DO TMDL Development
Draft
April 22, 2005
U.S. Environmental Protection Agency
Region 3
1650 Arch Street
Philadelphia, Pennsylvania
-------
CONTENTS
Page
1 INTRODUCTION 1-1
2 WATERSHED LOADING MODEL 2-1
2.1 HSPF Model Overview 2-1
2.2 XP-SWMM Model Overview 2-2
2.3 Modeling Assumptions 2-5
2.4 HSPF Model Configuration 2-5
2.4.1 HSPF Subbasins 2-5
2.4.2 Land Use Classifications 2-5
2.4.3 Nutrient Sources 2-11
2.4.4 Time Step and Simulation Duration 2-11
2.5 Model Testing and Calibration 2-11
3 EFDC HYDRODYNAMIC MODEL 3-1
3.1 General 3-1
3.2 Hydrodynamics and Salinity and Temperature Transport 3-4
3.3 Sediment Transport 3-5
3.4 Water Quality and Eutrophication Simulation 3-5
3.5 Toxic Contaminant Transport and Fate 3-6
3.6 Finfish and Shellfish Transport 3-6
3.7 Near-field Discharge Dilution and Mixing Zone Analysis 3-7
3.8 Spill Trajectory and Search and Rescue Simulation 3-7
3.9 Wetland, Marsh, and Tidal Flat Simulation 3-7
3.10 Nearshore Wave-induced Currents and Sediment Transport 3-8
3.11 User Interface 3-8
3.12 Preprocessing Software 3-9
3.13 Program Configuration 3-9
3.14 Run-Time Diagnostics 3-9
3.15 Model Output Options 3-9
3.16 Postprocessing, Graphics and Visualization 3-10
3.17 Documentation 3-10
3.18 Computer Requirements 3-10
4 EFDC WATER QUALITY MODEL 4-1
4.1 Introduction 4-1
4.2 Conservation of Mass Equation 4-4
4.3 Algae 4-6
4.4 Organic Carbon 4-11
4.5 Phosphorus 4-16
4.6 Nitrogen 4-21
4.7 Silica 4-26
4.8 Chemical Oxygen Demand 4-28
4.9 Dissolved Oxygen 4-29
4.10 Total Active Metal 4-31
4.11 Fecal Coliform Bacteria 4-32
4.12 Method of Solution 4-32
4.13 Macroalgae (Periphyton) State Variable 4-33
-------
5 EFDC SEDIMENT PROCESS MODEL 5-1
5.1 Depositional Flux 5-3
5.2 Diagenesis Flux 5-5
5.3 Sediment Flux 5-6
5.4 Silica 5-18
5.5 Sediment Temperature 5-19
5.6 Method of Solution 5-20
6 EFDC WATER QUALITY MODEL CALIBRATION 6-1
6.1 Modeling Assumptions 6-1
6.2 Model Configuration 6-1
6.2.1 Segmentation 6-3
6.2.2 Streamflow Estimation 6-3
6.2.3 Atmospheric and Tidal Boundary Conditions 6-3
6.2.4 Initial Conditions 6-5
6.2.5 Point and Nonpoint Source Representation 6-5
6.2.6 Time Step and Simulation Duration 6-7
6.3 Model Calibration Results 6-8
6.3.1 Tide Elevation and Phase 6-8
6.3.2 Water Depth and Stream Velocity 6-9
6.3.3 Sediment Oxygen Demand andBenthic Nutrient Fluxes 6-10
6.3.4 Water Quality Calibration Results 6-11
6.3.4.1 Mean Error Statistic 6-13
6.3.4.2 Absolute Mean Error Statistic 6-13
6.3.4.3 Root-Mean Error Statistic 6-13
6.3.4.4 Relative Error Statistic 6-14
6.3.4.5 Statistic Results 6-14
7 REFERENCES 7-1
Appendix A Model-Data Time-Series Graphics
Appendix B Summary of Wilmington CSO Volumes and Loads
Appendix C Miscellaneous Information
Appendix D Listing of Input Data Files for EFDC 1994-1998 Calibration Run
-------
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
1 -INTRODUCTION
A scientifically justifiable Total Maximum Daily Load (TMDL) for a waterbody can only be developed
based on a quantitative understanding of the system. In practice, water quality modeling offers a feasible
tool to establish this quantitative understanding. A water quality model that is customized for a specific
waterbody can simulate the major physical, chemical, and biological processes that occur in the system,
and thus provide quantitative relationships between the water quality response and external forcing
functions. A customized modeling framework was developed to support determination of nutrient and
dissolved oxygen TMDLs for the Christina River Basin. The TMDLs are presented in the report titled
Total Maximum Daily Loads for Nutrients and Dissolved Oxygen in the Christina River Basin,
Pennsylvania-Delaware-Maryland (USEPA, 2005). This report is intended to accompany the TMDL
report and provide a more detailed discussion on the models used for the nutrient TMDL analysis,
including assumptions, parameters, and references.
The modeling framework used in this study consisted of three major components: (1) a series of
watershed loading models (HSPF) developed for each of the four primary subwatersheds in the Christina
River Basin (Senior and Koerkle, 2003a, 2003b, 2003c, 2003d), (2) a CSO flow model (XP-SWMM)
developed by the City of Wilmington, and (3) a hydrodynamic model developed using the computational
framework of the Environmental Fluid Dynamics Code (EFDC) (Hamrick, 1992). A linkage interface
was also developed to allow for the transfer of model data results from the HSPF and XP-SWMM model
components to the EFDC water quality model.
Under the HSPF model framework, the Christina River Basin was configured into 70 subbasins (see
Figure 1-1 and Table 1-1) with each subbasin having 12 land use categories. The XP-SWMM model
calculated hourly CSO flow rates from rainfall events. Storm monitoring data were used to determine
event mean concentrations to estimate CSO loads for nutrients. The EFDC model framework includes
the main channels of Brandywine Creek, East Branch Brandywine Creek, West Branch Brandywine
Creek, Buck Run, Red Clay Creek, White Clay Creek, Christina River, Delaware River, and several other
smaller tributaries. The EFDC receiving water model was linked to the HSPF and XP-SWMM models to
incorporate watershed and CSO loads. The EFDC hydrodynamic and water quality model was used to
predict the dissolved oxygen and nutrient concentrations in the main channels of the Christina River,
Brandywine Creek, White Clay Creek, and Red Clay Creek watersheds. The water quality constituents
were calibrated using monitoring data for the period October 1, 1994 to October 1, 1998 (a period of 4
years). This period included two dry summers (1995 and 1997) as well as a number of high-flow
periods, both of which are important to satisfy the TMDL seasonality requirements.
1-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Table 1-1. Subbasins in the HSPF models of Christina River Basin
Subbasin
Stream Name
Area (mi2)
Brandywine Creek Watershed
601
B02
B03
B04
605
B06
B07
608
609
610
B11
612
B13
614
B15
B16
B17
618
619
620
621
622
623
624
625
626
627
628
629
B30
B31
B32
B33
B34
635
Upper 8randywine Creek West 6r.
8randywine Creek West 6r.
8randywine Creek West 6r.
8randywine Creek West 6r.
8randywine Creek West 6r.
8randywine Creek West 6r.
8randywine Creek West 6r.
8randywine Creek West 6r.
Upper 8randywine Creek East 6r.
8randywine Creek East 6r.
Brandywine Creek East Br.
Brandywine Creek East 6r.
Brandywine Creek East Br.
Brandywine Creek East 6r.
Brandywine Creek
Brandywine Creek
Brandywine Creek
Brandywine Creek
8randywine Creek
Upper 6uck Run
Upper Doe Run
Lower Doe Run
Lower 6uck Run
Trib. To 6road Run
6road Run
Marsh Creek
Marsh Creek
Trib. To Valley Creek
Valley Creek
Beaver Creek
Pocopson Creek
Birch Run
Rock Run
Lower Brandywine Creek
Upper Marsh Creek
18.39
7.38
6.76
0.80
8.82
8.06
13.46
3.62
14.68
18.31
6.31
3.70
7.94
12.92
10.36
14.06
7.51
10.37
8.64
25.54
11.05
10.96
1.95
0.60
5.83
2.61
11.54
2.40
18.21
18.08
9.19
4.66
8.03
6.05
5.80
Subbasin
Stream Name
Area (mi2)
White Clay Creek Watershed
W01
W02
W03
W04
W05
W06
W07
W08
W09
W10
W11
W12
W13
W14
W15
W16
W17
White Clay Creek West Br.
Upper White Clay Creek Middle 6r.
White Clay Creek Middle Br.
Trib. To White Clay Creek East 6r.
Trib. To White Clay Creek East Br.
Upper White Clay Creek East 6r.
Trout Run
White Clay Creek East 6r.
White Clay Creek East Br.
White Clay Creek
White Clay Creek
White Clay Creek
White Clay Creek
White Clay Creek
Muddy Run
Pike Creek
Mill Creek
10.23
9.51
6.35
6.20
2.65
8.57
1.37
7.47
6.85
3.58
6.53
8.76
2.08
3.41
3.89
6.65
13.00
Red Clay Creek Watershed
R01
R02
R03
R04
R05
R06
R07
R08
R09
Upper Red Clay Creek West Br.
Red Clay Creek West Br.
Red Clay Creek East Br.
Red Clay Creek
Red Clay Creek
6urroughs Run
Hoopes Reservoir
Red Clay Creek
Red Clay Creek
10.08
7.39
9.90
5.11
5.24
7.10
2.10
5.38
1.72
Christina River Watershed
C01
C02
COS
C04
COS
C06
C07
COS
C09
Christina River West 6r.
Upper Christina River
Christine River
Upper Little Mill Creek
Little Mill Creek
Muddy Run
8elltown Run
Christina River
Lower Christina River
6.70
9.73
4.47
5.37
3.84
8.64
6.37
10.70
21.90
1-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
N
KILOMETERS
• •
024
10
Q
Figure 1-1. Christina River Basin showing HSPF model subbasins and EFDC model grid
1-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
1-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
2 - WATERSHED LOADING MODELS
A watershed runoff and loading model (HSPF) was developed for the Christina River Basin to estimate
the amount of nutrients and oxygen demanding substances introduced to the receiving streams during
rainfall-runoff events. In addition, an urban storm water runoff model (XP-SWMM) was developed by
the City of Wilmington and was used to estimate combined sewer overflow (CSO) flows and loads to
local receiving waters.
2.1 HSPF Model Overview
The Hydrologic Simulation Program—Fortran (HSPF), is a U.S. EPA supported model for simulation of
watershed hydrology and water quality for both conventional and toxic organic pollutants. The HSPF
model uses information such as the time history of rainfall, temperature and solar radiation; land surface
characteristics such as land-use patterns; and land management practices to simulate the processes that
occur in a watershed. The result of this simulation is a time history of the quantity and quality of runoff
from an urban or agricultural watershed. Flow rate, sediment load, and nutrient and pesticide
concentrations are predicted. HSPF includes an internal database management system to process the large
amounts of simulation input and output. HSPF includes the source code, executable version, user's guide,
and technical support. The HSPF model incorporates the watershed-scale Agricultural Runoff Model
(ARM) and Non-Point Source (NPS) models into a basin-scale analysis framework that includes
pollutant transport and transformation in stream channels.
The Christina River Basin drains 565 square miles in Pennsylvania, Delaware, and Maryland. Water
from the basin is used for recreation, drinking-water supply, and to support aquatic life. The Christina
River Basin includes four main watersheds: Brandywine Creek, Red Clay Creek, White Clay Creek, and
Christina River. Brandywine Creek is the largest of the watersheds and drains an area of 327 square
miles. Water quality in some parts of the Christina River Basin is impaired and does not support
designated uses of the streams.
A multi-agency water-quality management strategy included a modeling component to evaluate the
effects of point and nonpoint-source contributions of nutrients and suspended sediment on stream water
quality. To assist in nonpoint-source evaluation, four independent models, one for each of the four main
watersheds of the Christina River Basin, were developed and calibrated using the HSPF modeling
framework.
The HSPF models simulate streamflow, suspended sediment, nitrogen, phosphorus, BOD, water
temperature, and dissolved oxygen. For the models, the Christina River Basin was subdivided into 70
2-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
reaches. Twelve different pervious land uses and two impervious land uses were selected for simulation.
Land-use areas were determined from 1995 land-use data. The predominant land uses in the basin are
forested, agricultural, residential, and urban.
The hydrologic component of the model was run at an hourly time step and calibrated using streamflow
data for eight U.S. Geological Survey (USGS) streamflow measurement stations for a period covering
four water years from October 1, 1994 to October 1, 1998. Daily precipitation data for three National
Oceanic and Atmospheric Administration (NOAA) gages and hourly data for one NOAA gage were used
for model input. More detailed descriptions of the HSPF models developed for the Christina River Basin
can be found in Senior and Koerkle (2003a, 2003b, 2003c, and 2003d).
2.2 XP-SWMM Model Overview
The City of Wilmington has developed a model (XP-SWMM) to simulate stormwater flows and CSO
events in the city's sewer collection system. XP-SWMM is a link-node model that performs hydrology,
hydraulics, and water quality analysis of stormwater and wastewater drainage systems including sewage
treatment plants, water quality control devices, and best management practices (BMPs). XP-SWMM can
be used to model the full hydrologic cycle from stormwater and wastewater flow and pollutant generation
to simulation of the hydraulics in any combined system of open and/or closed conduits with any
boundary conditions.
Typical XP-SWMM applications include predicting combined sewer overflows (CSOs) and sanitary
sewer overflows (SSOs), interconnected pond analysis, open and closed conduit flow analysis,
major/minor flow analysis, design of new developments, and analysis of existing stormwater and sanitary
sewer systems.
XP-SWMM uses a self-modifying dynamic wave solution algorithm. Like all implicit solutions, which
solve for the unknown values at a given time simultaneously, XP-SWMM is not Courant-limited.
However, XP-SWMM uses the Courant number as a guide, to prevent numerical attenuation that can
occur if excessively large time steps are used. This is important in models where pumps are involved or
in urban systems where steeply rising hydrographs, requiring responses in seconds or fractions of a
second will predominate, or where checks are being made against empirical procedures like the FHWA
inlet control scheme for culverts. XP-SWMM will use small time steps when required and larger time
steps when appropriate.
XP-SWMM has three computational modules. There is a stormwater module for hydrology and water
quality generation, a wastewater module for generation of wastewater flows including Storage/Treatment
2-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
for BMP and water quality routing, and a hydrodynamic hydraulics module for the hydraulic simulation
of open and closed conduit wastewater or stormwater systems.
Hourly flow rates at each of the city's 38 CSO outfalls were calculated by XP-SWMM for the 1994-1998
calibration period based on hourly rainfall measured at New Castle County Airport and Porter Reservoir.
Water quality was monitored at three CSO locations (CSO 25, CSO 4b, and the 11th Street Pump
Station) for two storm events on 10/27/2003 and 12/17/2003. Event mean concentrations (EMCs) were
estimated for nutrients and oxygen demanding substances (see Tables 2-la, b, and c). The monitoring
included 20-day CBOD (CBOD20), 5-day CBOD (CBOD5), dissolved organic carbon (DOC), total
organic carbon (TOC), ammonia nitrogen (NH3-N), nitrite+nitrate nitrogen (NOxN), total Kjeldahl
nitrogen (TKN), total nitrogen (TN), dissolved orthophosphate (DOrthP), total phosphorus (TP), and
total suspended solids (TSS). The EMCs were used in conjunction with the CSO flow rates to estimate
daily loads for each CSO outfall. The CSO flows and loads were then input to the EFDC receiving water
model to simulate their impact on nutrient and dissolved oxygen concentrations in the tidal Christina
River and tidal Brandywine Creek. The annual loads from each of the CSO outfalls for the calibration
period are tabulated in Appendix B. The locations of the CSOs are shown in Appendix B, Figure B-l.
Table 2-1 a. Storm monitoring at Wilmington CSO 4b.
Date
1 0/27/2003
1 0/27/2003
1 0/27/2003
1 0/27/2003
1 0/27/2003
1 0/27/2003
Time
11:40
12:10
12:40
13:10
13:40
14:10
CBOD20
mg/L
14.62
13.60
10.20
14.48
13.98
13.50
CBOD5
mg/L
11.70
5.82
5.64
7.85
7.65
10.60
DOC
mg/L
6.6
2.9
6.1
5.9
6.8
7.3
TOC
mg/L
9.1
3.7
6.2
7.1
8.3
8.9
NH3-N
mg/L
0.362
0.137
0.189
0.238
0.244
0.238
NOxN
mg/L
0.969
0.248
0.502
0.831
1.070
1.290
TKN
mg/L
1.400
0.275
0.644
1.080
1.210
1.370
TN
mg/L
2.369
0.248
0.502
1.911
2.280
2.660
DOrthP
mg/L
0.004
0.020
0.100
0.126
0.141
0.159
TP
mg/L
0.238
0.320
0.219
0.270
0.219
0.216
TSS
mg/L
298
278
195
177
75
32
12/17/2003
12/17/2003
12/17/2003
12/17/2003
12/17/2003
12/17/2003
09:00
09:30
10:00
10:30
11:00
11:30
16.20
16.10
23.80
16.20
12.10
10.60
9.20
8.65
12.80
10.60
8.18
6.86
4.9
4.7
6.8
5.9
5.5
5.0
6.8
6.2
8.4
6.1
6.0
6.2
0.403
0.480
4.520
0.504
0.486
0.357
0.627
0.855
1.210
1.360
1.710
1.970
2.650
2.790
4.830
3.060
2.610
1.950
3.277
3.645
6.040
4.420
4.320
3.920
0.203
0.180
0.222
0.192
0.138
0.112
0.388
0.382
0.546
0.416
0.306
0.194
35
34
25
17
19
19
EMC 10/27/2003
EMC 12/17/2003
EMC both storms
13.40
15.83
14.62
8.21
9.38
8.80
5.94
5.47
5.70
7.22
6.62
6.92
0.235
1.125
0.680
0.818
1.289
1.054
0.997
2.982
1.989
1.662
4.270
2.966
0.092
0.175
0.133
0.247
0.372
0.310
176
25
100
2-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Table 2-1 b. Storm monitoring at Wilmington CSO 25.
Date
1 0/27/2003
1 0/27/2003
1 0/27/2003
1 0/27/2003
1 0/27/2003
1 0/27/2003
Time
11:00
11:30
12:00
12:30
13:30
14:00
CBOD20
mg/L
13.88
14.76
7.83
12.14
14.10
14.26
CBOD5
mg/L
13.88
14.76
5.36
12.14
14.10
14.26
DOC
mg/L
11.8
10.3
3.8
70.5
10.6
10.8
TOC
mg/L
14.4
11.6
4.3
80.0
11.6
12.0
NH3-N
mg/L
0.325
0.294
0.136
0.421
0.352
0.455
NOxN
mg/L
0.516
0.503
0.215
0.634
0.820
1.160
TKN
mg/L
1.270
1.050
0.392
3.070
1.900
2.480
TN
mg/L
1.786
1.553
0.215
3.704
2.720
3.640
DOrthP
mg/L
0.234
0.286
0.113
1.870
0.249
0.354
TP
mg/L
0.296
0.397
0.178
1.620
0.450
0.642
TSS
mg/L
32
33
51
39
26
15
12/17/2003
12/17/2003
12/17/2003
08:45
09:15
09:45
15.00
28.30
28.76
9.48
19.60
28.76
6.3
9.1
40.8
6.6
10.2
44.6
0.350
0.500
3.720
0.547
0.839
1.030
1.850
3.140
5.500
2.397
3.979
6.530
0.202
0.317
1.560
0.102
0.296
1.580
27
22
14
EMC 10/27/2003
EMC 12/17/2003
EMC both storms
12.83
24.02
16.56
12.42
19.28
14.70
19.65
18.73
19.34
22.32
20.47
21.70
0.331
1.523
0.728
0.641
0.805
0.696
1.694
3.497
2.295
2.270
4.302
2.947
0.518
0.693
0.576
0.597
0.659
0.618
33
21
29
Table 2-1 c. CSO Storm monitoring at Wilmington 11th Street Pumping Station
Date
10/27/2003
10/27/2003
10/27/2003
10/27/2003
10/27/2003
10/27/2003
Time
11:20
11:50
12:10
12:50
13:20
13:50
CBOD20
mg/L
11.76
10.88
10.88
12.98
11.82
11.66
CBOD5
mg/L
11.76
10.88
10.88
9.02
11.82
11.66
DOC
mg/L
23.5
9.5
7.7
4.6
13.9
6.8
TOC
mg/L
29.6
11.9
9.6
5.8
15.3
8.5
NH3-N
mg/L
4.040
3.070
1.520
2.200
1.720
2.340
NOxN
mg/L
0.467
1.100
0.545
0.517
0.646
0.753
TKN
mg/L
7.250
3.820
1.450
1.400
0.964
1.880
TN
mg/L
7.717
4.920
1.995
1.917
1.610
0.753
DOrthP
mg/L
0.262
0.433
0.202
0.003
0.167
0.311
TP
mg/L
1.470
0.520
0.357
0.366
0.289
0.420
TSS
mg/L
454
71
166
144
104
106
12/17/2003
12/17/2003
12/17/2003
12/17/2003
12/17/2003
12/17/2003
08:50
09:20
09:50
10:20
10:50
11:20
82.32
26.50
29.60
20.80
42.40
82.05
29.30
13.80
15.40
14.30
23.70
82.05
8.5
5.3
6.0
6.7
7.3
21.4
10.4
6.3
8.2
9.1
11.3
25.5
3.040
4.520
1.650
3.530
2.940
1.150
0.682
0.732
0.820
0.842
1.200
1.140
6.790
4.880
4.900
4.670
5.910
6.810
7.472
0.732
5.720
5.512
7.110
7.950
0.157
0.129
0.004
0.019
0.004
0.341
1.160
0.630
0.632
0.645
0.883
0.909
143
86
91
73
106
64
EMC 10/27/2003
EMC 12/17/2003
EMC both storms
11.66
47.28
29.47
11.00
29.76
20.38
11.01
9.20
10.10
13.45
11.80
12.63
2.482
2.805
2.643
0.671
0.903
0.787
2.794
5.660
4.227
3.152
5.749
4.451
0.230
0.109
0.169
0.570
0.810
0.690
174
94
134
2-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
2.3 Modeling Assumptions
The simulation of streamflow in the Christina River Basin HSPF models considered the following
assumptions: (1) inputs of hourly precipitation would be estimated reasonably well by disaggregated 24-
hour precipitation data; (2) the average precipitation over a given land segment would be represented
adequately by weighted data from a single precipitation gage; and (3) a simplified set of impervious land
uses (PERLND) and impervious land uses (IMPLND) would not limit a satisfactory hydrologic
calibration (Senior and Koerkle, 2003a).
The simulation of water quality in the HSPF models considered the following assumptions: (1) land-
based contributions of sediment and nutrients could be simulated by a simplified set of land-use
categories; (2) water quality could be represented by the condition where chemical transformation of
nutrients are simulated explicitly in the stream channel but not in land processes; and (3) the contribution
of sediment from bank erosion in the stream channel can be estimated by sediment from pervious land
areas (Senior and Koerkle, 2003a).
The simulation of CSO nutrient loads assumes that the event mean concentrations are the same no matter
what the intensity or duration of the storm event. Also, since CSO concentrations were monitored only at
two outfalls (CSO 4b and CSO 25), the concentrations at the remaining 36 CSO outfalls were assumed to
be equivalent to those measured at the 11th Street Pump Station.
2.4 HSPF Model Configuration
2.4.1 HSPF Subbasins
Four separate HSPF models were developed to simulate watershed runoff and nutrient loading in the
Christina River Basin. One model was developed for each of the four main watersheds: Brandywine
Creek watershed, White Clay Creek watershed, Red Clay Creek watershed, and Christina River
watershed. The Christina River Basin was delineated into 70 subbasins (or reaches) for the modeling
effort (see Figure 1-1). The size of the subbasins ranged from 0.6 to 25.5 mi2. The subbasins were
delimited based on major tributary inflows, calibration locations (stream gages and water quality
monitoring stations), and time-of-travel considerations.
2.4.2 Land Use Classifications
Spatial data input to the HSPF model are used to define the structure fixed characteristics of the model.
The principal structural unit of the HSPF model is the hydrologic response units PERLND (pervious
land) and IMPLND (impervious land). Fifteen original land-use categories (circa 1995) from several
2-5
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
sources were simplified and reclassified into 10 pervious and 2 impervious land-use categories that were
expected to have distinct nonpoint-source water-quality characteristics (Table 2-2).
Agricultural land use was divided into three characteristic subtypes for the model. Agricultural-livestock
land use identifies relatively small acreage farms with high animals-per-acre densities, limited pasture
areas, and rowcrops. Small acreage dairy operations typify this land-use type. Agricultural-rowcrop land
use identifies farms with lower animals-per-acre densities (typically beef cattle and horses) and
substantial pasture and crop acreage. Agricultural-mushroom land use is the third type of agriculture land
use delimited, but mushroom production operations are much more prevalent in the Red Clay Creek and
White Clay Creek Basins than in the Brandywine Creek Basin. Residential land use is distributed
throughout the basin and is divided into two types: sewered and non-sewered. Sewered residential areas
tend to have higher housing densities and are nearer to urban/suburban areas than non-sewered area.
Non-sewered residential areas tend to have lower densities and are more rural. Other urban land use is in
small boroughs and along major roadways. Forested land is distributed throughout the basin and tends to
be along stream channels. The land use delineations for each of the four main watersheds in the
Christina River Basin are presented in Tables 2-3 through 2-6.
Table 2-2. Land-use categories used in HSPF models for Christina River Basin
Land-use category for HSPF model
Pervious
Impervious
Residential-septic
Residential-sewer
Urban
Agricultural-livestock
Agricultural-rowcrop
Agricultural-mushroom
Open
Forested
Wetlands/water
Undesignated
Residential
Urban
Description
Residential land not within a sewer service area
Residential land within a sewer service area
Commercial, industrial, institutional, and transportation uses
Predominantly mixed agricultural activities of dairy cows, pasture,
and other livestock operations
Predominantly row crop cultivation (corn, soybean, alfalfa), may
include some hay or pasture land
Mushroom-growing activities including compost preparation,
mushroom-house operations, spent compost processing
Recreational and other open land not used for agricultural
Predominantly forested land
Wetlands and open water
Land use not defined
Impervious residential land
Impervious commercial, industrial, and other urban land
2-6
-------
Table 2-3. Land use characteristics (in percent) for Brandywine Creek Watershed.
Reach
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
Total
Length
(mi)
6.60
7.60
2.94
1.85
2.91
2.93
7.80
2.19
7.10
12.10
1.79
2.02
3.86
4.86
2.49
2.88
4.15
3.39
2.71
8.66
6.73
3.18
0.87
3.14
3.14
1.60
4.80
2.00
7.20
4.09
4.09
2.00
2.75
4.46
4.00
1 44.88
Area
(sq.mi.)
18.39
7.38
6.76
0.80
8.82
8.06
13.46
3.62
14.68
18.31
6.31
3.70
7.94
12.92
10.36
14.06
7.51
10.37
8.64
25.54
11.05
10.96
1.95
0.60
5.83
2.61
11.54
2.40
18.21
18.08
9.19
4.66
8.03
6.05
5.80
324.59
Resident.
Septic
4.1
17.3
22.3
0.0
1.5
17.1
5.9
9.2
6.1
17.0
4.7
8.4
6.5
8.8
17.6
25.0
12.2
9.2
10.6
7.7
3.5
0.7
0.0
73.2
15.2
8.1
21.5
0.1
4.3
12.2
22.7
11.3
12.2
1.9
6.3
10.5
Resident.
Sewer
1.4
0.6
0.3
7.1
11.0
0.5
0.0
0.0
0.5
0.2
11.6
18.7
10.1
10.8
7.2
0.0
0.0
3.5
10.3
1.8
0.0
0.0
0.0
4.9
3.7
0.0
0.1
37.6
12.9
6.6
0.0
0.0
3.5
2.5
0.0
3.9
Urban
0.6
1.8
1.2
2.1
10.5
1.5
1.5
0.6
0.4
1.2
1.9
4.4
4.3
3.5
1.9
2.4
0.1
1.6
3.4
1.1
0.4
0.9
0.01
0.0
2.3
2.2
0.9
6.5
3.5
4.7
0.8
0.8
1.3
28.0
1.1
2.7
Agriculture
livestock
45.6
9.4
7.5
0.0
0.0
4.0
0.0
0.0
27.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
2.1
0.0
5.9
7.6
7.9
4.9
0.0
0.0
6.5
8.9
0.0
0.0
0.0
0.0
15.8
4.2
0.0
12.1
6.3
Agriculture
row crop
22.5
19.0
22.6
14.9
19.1
35.6
49.0
62.6
27.0
36.0
33.1
11.4
14.3
31.9
40.7
25.7
27.0
19.1
4.1
52.9
68.6
71.3
44.4
3.5
40.7
19.6
20.6
3.0
20.9
32.4
48.8
15.8
38.0
1.6
36.3
32.7
Agriculture
mushroom
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
Forested
20.1
46.4
39.8
68.8
34.8
35.4
38.2
24.9
32.6
40.3
35.6
38.9
47.9
30.2
16.8
38.7
48.6
38.2
16.5
25.5
17.3
17.7
49.4
8.2
30.4
59.5
33.9
20.5
35.1
30.0
22.1
52.9
29.8
13.9
34.1
31.8
Open
2.7
0.5
2.0
0.1
3.6
1.8
1.9
0.0
2.2
1.2
4.7
2.2
3.1
3.2
6.9
1.6
6.1
14.6
40.3
1.3
1.1
0.0
0.0
0.0
1.7
0.3
2.4
5.7
5.0
2.2
1.8
0.9
4.5
12.9
0.3
3.8
Wetland-
water
0.5
0.8
0.5
1.7
1.5
0.5
1.2
1.2
2.8
0.6
0.5
1.3
1.4
1.0
1.0
0.9
1.3
1.1
1.0
0.4
0.1
0.3
1.3
0.0
0.1
0.5
7.4
0.0
2.5
0.2
0.3
0.1
2.1
2.6
7.8
1.2
Undesig-
nated
0.9
0.2
0.0
0.2
2.4
0.0
0.1
0.1
0.2
0.2
0.6
1.3
2.7
1.4
1.0
0.5
0.4
5.8
4.6
0.8
0.5
0.2
0.0
0.0
0.3
0.1
1.1
3.6
3.2
2.7
0.3
0.3
0.4
7.3
0.2
1.3
Impervious-
residential
1.1
2.2
2.6
3.0
4.9
2.1
0.7
1.0
0.9
2.0
5.5
8.9
5.1
5.6
5.0
2.8
4.1
2.5
5.6
1.6
0.4
0.1
0.0
10.3
3.3
0.9
2.4
16.1
6.0
4.2
2.5
1.3
2.9
1.3
0.7
2.9
Impervious-
urban
0.7
1.8
1.2
2.1
10.7
1.5
1.5
0.6
0.4
1.2
1.9
4.5
4.6
3.6
1.9
2.4
0.3
2.2
3.6
1.1
0.4
0.9
0.01
0.0
2.3
2.2
0.9
6.7
6.7
5.0
0.8
0.8
1.3
28.2
1.1
2.8
a
a
b
o
9
-------
fable 2-4. Land use characteristics (in percent
Reach
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Total
Length
(mi)
7.33
6.57
7.18
6.02
2.49
6.16
1.75
4.09
4.46
1.67
4.02
5.28
2.21
2.97
4.08
5.85
9.76
81.89
Area
(sq.mi.)
10.23
9.51
6.35
6.20
2.65
8.57
1.37
7.47
6.85
3.58
6.53
8.76
2.08
3.41
3.89
6.65
13.00
107.10
Resident.
Septic
15.6
11.3
16.4
6.8
1.5
1.5
5.8
11.6
17.5
11.2
1.2
0.0
0.0
0.0
0.0
0.0
0.5
6.5
Resident.
Sewer
0.0
1.8
0.0
2.6
0.0
0.8
5.1
0.5
7.2
4.5
8.1
24.1
6.7
10.6
14.4
38.9
33.9
11.1
Urban
1.0
0.8
0.0
1.3
0.0
1.3
1.5
0.5
0.7
0.0
4.3
10.2
14.4
11.4
1.8
6.2
6.6
3.4
for White Clay Creek watershed.
Agriculture
livestock
10.4
15.9
9.0
11.5
14.7
13.4
0.0
0.0
6.6
5.3
0.0
0.0
0.0
0.0
0.0
0.0
1.1
5.8
Agriculture
row crop
36.2
47.5
33.5
40.2
52.1
47.3
5.8
20.1
22.9
21.8
15.5
9.4
11.1
0.0
29.6
8.4
8.7
25.3
Agriculture
mushroom
5.2
0.0
2.2
5.8
7.5
6.8
56.2
30.3
3.2
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.1
4.9
Forested
26.1
17.7
35.9
23.5
23.0
21.8
17.5
32.4
31.7
53.1
54.8
10.7
11.5
14.1
42.2
12.9
11.7
24.8
Open
2.1
1.1
0.6
2.4
0.8
2.9
0.0
1.2
2.3
0.3
7.0
9.8
7.2
21.7
3.3
9.2
10.2
4.9
Wetland-
water
0.1
0.2
0.5
0.5
0.0
0.2
1.5
0.5
1.6
0.6
0.8
0.9
1.4
14.4
0.0
0.0
0.0
0.9
Undesig-
nated
0.7
0.9
0.0
2.1
0.4
2.1
2.9
0.8
0.6
0.3
0.3
10.2
27.4
10.6
0.5
1.4
4.1
2.9
Impervious-
residential
1.8
2.0
1.9
1.9
0.0
0.5
2.9
1.5
5.0
3.1
3.7
10.4
2.9
4.7
6.2
16.7
15.1
5.5
Impervious-
urban
1.0
0.8
0.0
1.5
0.0
1.3
1.5
0.5
0.7
0.0
4.3
14.4
17.3
12.6
2.1
6.3
7.1
4.0
a
a
b
o
9
-------
Table 2-5. Land use characteristics (in percent) for Red Clay Creek watershed
Reach
1
2
3
4
5
6
7
8
9
Total
Length
(mi)
5.00
4.90
7.20
3.40
5.10
5.00
1.70
4.30
0.84
37.44
Area
(sq.mi.)
10.08
7.39
9.90
5.11
5.24
7.10
2.10
5.38
1.72
54.02
Resident.
Septic
10.4
11.8
14.9
35.5
32.2
23.9
26.1
1.5
0.0
17.1
Resident.
Sewer
1.6
0.5
1.9
2.3
2.2
0.0
0.0
35.3
45.5
6.1
Urban
2.1
0.7
1.0
1.0
0.2
0.2
0.5
5.3
4.8
1.6
Agriculture
livestock
5.9
0.0
0.0
1.8
0.0
0.0
0.0
0.0
0.0
1.3
Agriculture
row crop
46.8
40.9
33.1
14.2
14.6
42.4
7.2
1.6
0.0
29.2
Agriculture
mushroom
5.9
17.5
14.2
1.8
0.0
0.0
0.0
0.0
0.0
6.3
Forested
18.4
25.2
22.5
28.4
37.7
25.1
44.4
13.4
10.2
24.0
Open
2.8
0.3
5.8
7.9
6.0
5.1
3.0
13.8
3.8
5.2
Wetland-
water
0.4
0.2
0.6
0.8
1.2
0.1
14.3
1.4
0.4
1.2
Undesig-
nated
1.6
0.7
2.3
0.6
0.9
0.3
1.1
6.3
9.9
2.0
Impervious-
residential
1.9
1.5
2.5
4.9
4.6
2.7
2.9
15.3
19.5
4.5
Impervious-
urban
2.2
0.7
1.1
1.0
0.2
0.2
0.5
6.0
5.9
1.7
a
a
b
o
9
-------
Table 2-6. Land use characteristics (in percent) for Christina River watershed
Reach
1
1
2
2
3
4
5
6
7
8
9
Total
Length
(mi)
1.46
3.33
4.96
3.20
3.24
4.65
2.30
6.73
2.08
6.21
15.09
53.25
Area
(sq.mi.)
1.87
4.83
7.08
2.65
4.47
5.37
3.84
8.64
6.37
10.70
21.90
77.70
Resident.
Septic
9.1
7.9
25.5
7.7
4.8
1.5
0.1
6.0
0.6
0.2
0.5
4.6
Resident.
Sewer
0.0
4.0
4.7
25.9
17.5
28.1
20.9
8.8
20.1
19.3
10.3
13.7
Urban
0.4
7.9
0.3
9.6
12.7
12.8
17.6
6.6
6.9
8.4
17.8
10.8
Agriculture
livestock
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
Agriculture
row crop
58.1
27.1
41.1
4.5
9.2
1.4
0.0
15.5
8.5
10.4
1.8
11.9
Agriculture
mushroom
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
Forested
29.7
23.4
21.7
10.7
18.4
18.3
8.7
38.4
34.3
27.0
12.5
21.6
Open
1.1
12.5
0.6
13.8
8.0
9.2
15.2
10.6
8.1
7.0
20.0
11.6
Wetland-
water
0.0
0.0
0.1
0.0
0.2
0.2
1.1
1.1
0.8
0.5
4.3
1.6
Undesig-
nated
0.2
6.2
1.0
13.3
9.0
2.6
9.3
1.8
4.7
9.3
9.9
6.5
Impervious-
residential
1.0
2.6
4.8
12.0
8.0
12.2
8.9
4.5
8.7
8.3
4.5
6.4
Impervious-
urban
0.4
8.3
0.3
10.2
13.5
13.4
25.6
6.8
7.4
9.2
18.5
11.4
a
a
b
o
9
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
2.4.3 Nutrient Sources
The HSPF models required large amounts of data to characterize the hydrologic and water quality
response of the watershed to precipitation and other inputs. Data used in creating the model structure
and parameters were derived primarily from spatial analysis of basin characteristics and other published
information. Spatial data analyzed for model construction included land use, land-surface slope, and soil
associations. Time-series inputs for streamflow and water-quality simulation included meteorologic,
precipitation quality, water-use, and point source quantity and quality data. Nonpoint sources of
nutrients were calculated by the model based on build up, storage, and wash off processes inherent in the
HSPF model.
2.4.4 Time Step and Simulation Duration
The HSPF models were executed on a 1-hour time step. The duration of the calibration runs was from
October 1, 1994 to October 1, 1998, a period that covered 4 consecutive water years.
2.5 Model Testing and Calibration
Complete descriptions of the calibration of each of the four HSPF models for the Christina River Basin
can be found in the USGS Water Resources Investigation Reports (Senior and Koerkle, 2003a, 2003b,
2003c, and 2003d), which are available in Portable Document File (PDF) format at the following
website: http://pa.water.usgs.gov/pa pubs.html
2-11
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
2-12
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
3 - EFDC HYDRODYNAMIC MODEL
Modeling the physics, chemistry, and biology of the receiving waters of streams, lakes, estuaries, or
coastal regions requires a model that incorporates all the major processes. Transport processes for this
study were simulated using the three-dimensional EFDC hydrodynamic model that includes temperature
transport. The EFDC hydrodynamic model was developed by Hamrick (1992a). The model formulation
was based on the principles expressed by the equations of motion, conservation of volume, and
conservation of mass. Quantities computed by the model included three-dimensional velocities, surface
elevation, vertical viscosity and diffusivity, temperature, salinity, and density.
3.1 General
The Environmental Fluid Dynamics Code is a general purpose modeling package for simulating
three-dimensional flow, transport, and biogeochemical processes in surface water systems including
rivers, lakes, estuaries, reservoirs, wetlands, and coastal regions. The EFDC model was originally
developed at the Virginia Institute of Marine Science for estuarine and coastal applications and is
considered public domain software. In addition to hydrodynamic and salinity and temperature transport
simulation capabilities, EFDC is capable of simulating cohesive and noncohesive sediment transport,
near field and far field discharge dilution from multiple sources, eutrophication processes, the transport
and fate of toxic contaminants in the water and sediment phases, and the transport and fate of various life
stages of finfish and shellfish. Special enhancements to the hydrodynamic portion of the code, including
vegetation resistance, drying and wetting, hydraulic structure representation, wave-current boundary
layer interaction, and wave-induced currents, allow refined modeling of wetland marsh systems,
controlled flow systems, and near-shore wave induced currents and sediment transport. The EFDC model
has been extensively tested and documented for more than 20 modeling studies. The model is presently
being used by a number of organizations including universities, governmental agencies, and
environmental consulting firms.
The structure of the EFDC model includes four major modules: (1) a hydrodynamic model, (2) a water
quality model, (3) a sediment transport model, and (4) a toxics model (see Figure 3-1). The EFDC
hydrodynamic model itself, which was used for this study, is composed of six transport modules
including dynamics, dye, temperature, salinity, near field plume, and drifter (see Figure 3-2). Various
products of the dynamics module (i.e., water depth, velocity, and mixing) are directly coupled to the
water quality, sediment transport, and toxics models as shown in the following figures. Schematic
diagrams for the water quality model, the sediment transport model, and the toxics model are shown in
Figures 3-3, 3-4, and 3-5, respectively.
3-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
EFDC Model
Hydrodynamics
Water
Quality
Sediment
Transport
Toxics
Figure 3-1. Primary modules of the EFDC model
Dynamics ;; D
(E, u, v, w, mixing) ^ y
\
Hydrodynamics ;
• •"..,;.
•:: Temperature
,"i. .-;..•;..•;..•;..•;. .-i-V
„ .. ., Near Field _ .,,
Sa mity _. Drifter
* Pume
Figure 3-2. Structure of the EFDC hydrodynamic model
Hydrodynamic V
Model
Figure 3-3. Structure of the EFDC water quality model
3-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Hydrodynamic
Model
Sediment Transport
Model
Water Column
Sediment Bed
Cohesive
Noncohesive
Cohesive
Noncohesive
Figure 3-4. Structure of the EFDC sediment transport model
Hydrodynamic
Model
Sediment Transport
Model
£>..
Cohesive
Sediment
Toxic Model
Figure 3-5. Structure of the EFDC toxics model
3-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
3.2 Hydrodynamics and Salinity and Temperature Transport
The physics of the EFDC model and many aspects of the computational scheme are equivalent to the
widely used Blumberg-Mellor model (Blumberg and Mellor 1987). The EFDC model solves the
three-dimensional, vertically hydrostatic, free surface, turbulent averaged equations of motions for a
variable density fluid. Dynamically coupled transport equations for turbulent kinetic energy, turbulent
length scale, salinity, and temperature are also solved. The two turbulence parameter transport equations
implement the Mellor-Yamada level 2.5 turbulence closure scheme (Mellor and Yamada 1982; Galperin
et al. 1988). The EFDC model uses a stretched or sigma vertical coordinate and Cartesian, or curvilinear,
orthogonal horizontal coordinates.
The numerical scheme employed in EFDC to solve the equations of motion uses second order accurate
spatial finite differencing on a staggered or C grid. The model's time integration employs a second order
accurate three-time level, finite difference scheme with an internal-external mode splitting procedure to
separate the internal shear or baroclinic mode from the external free surface gravity wave or barotropic
mode. The external mode solution is semi-implicit and simultaneously computes the two-dimensional (2-
D) surface elevation field by a preconditioned conjugate gradient procedure. The external solution is
completed by the calculation of the depth average barotropic velocities using the new surface elevation
field. The model's semi-implicit external solution allows large time steps that are constrained only by the
stability criteria of the explicit central difference or high order upwind advection scheme (Smolarkiewicz
and Margolin 1993) used for the nonlinear accelerations. Horizontal boundary conditions for the external
mode solution include options for simultaneously specifying the surface elevation only, the characteristic
of an incoming wave (Bennett and Mclntosh 1982), free radiation of an outgoing wave (Bennett 1976;
Blumberg and Kantha 1985), or the normal volumetric flux on arbitrary portions of the boundary. The
EFDC model's internal momentum equation solution, at the same time step as the external solution, is
implicit with respect to vertical diffusion. The internal solution of the momentum equations is in terms
of the vertical profile of shear stress and velocity shear, which results in the simplest and most accurate
form of the baroclinic pressure gradients and eliminates the over determined character of alternate
internal mode formulations. Time splitting inherent in the three-time-level scheme is controlled by
periodic insertion of a second order accurate two-time-level trapezoidal step. EFDC is also readily
configured as a 2-D model in either the horizontal or vertical planes.
The EFDC model implements a second order accurate in space and time, mass conservation fractional
step solution scheme for the Eulerian transport equations for salinity, temperature, suspended sediment,
water quality constituents, and toxic contaminants. The transport equations are temporally integrated at
the same time step or twice the time step of the momentum equation solution (Smolarkiewicz and
Margolin 1993). The advective step of the transport solution uses either the central difference scheme
3-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
used in the Blumberg-Mellor model or a hierarchy of positive definite upwind difference schemes. The
highest accuracy upwind scheme, second order accurate in space and time, is based on a flux-corrected
transport version Smolarkiewicz's multidimensional positive-definite advection transport algorithm
(Smolarkiewicz and Clark, 1986; Smolarkiewicz and Grabowski 1990), which is monotonic and
minimizes numerical diffusion. The horizontal diffusion step, if required, is explicit in time, whereas the
vertical diffusion step is implicit. Horizontal boundary conditions include time variable material inflow
concentrations, upwind outflow, and a damping relaxation specification of climatological boundary
concentration. The NOAA Geophysical Fluid Dynamics Laboratory's atmospheric heat exchange model
(Rosati and Miyakoda 1988) is implemented for the temperature transport equation.
3.3 Sediment Transport
The EFDC code is capable of simulating the transport and fate of multiple size classes of cohesive and
noncohesive suspended sediment including bed deposition and resuspension. Water column transport is
based on the same high order advection-diffusion scheme used for salinity and temperature. A number of
options are included for the specification of settling velocities. For the transport of multiple size classes
of cohesive sediment, an optional flocculation model (Durban et al. 1989, 1990) can be activated.
Sediment mass conservative deposited bed formulations are included for both cohesive and noncohesive
sediment. The deposited bed may be represented by a single layer or multiple layers. The multiple bed
layer option provides a time since deposition versus vertical position in the bed relationship to be
established. Water column/sediment bed interface elevation changes can be optionally incorporated into
the hydrodynamic continuity equation. An optional one-dimensional (1-D) in the vertical, bed
consolidation calculation can be performed for cohesive beds.
3.4 Water Quality and Eutrophication Simulation
The EFDC code includes two internal eutrophication submodels for water quality simulation (Park et al.
1995). The simple or reduced eutrophication model is functionally equivalent to the WASPS EUTRO
model (Ambrose et al. 1993). The complex or full eutrophication model is functionally equivalent to the
CE-QUAL-ICM or Chesapeake Bay Water Quality model (Cerco and Cole 1993). Both water column
eutrophication models are coupled to a functionally equivalent implementation of the CE-QUAL-ICM
sediment diagenesis or biogeochemical processes model (DiToro and Fitzpatrick 1993). The
eutrophication models can be executed simultaneously with the hydrodynamic component of EFDC, or
EFDC simulated hydrodynamic transport fields can be saved, allowing the EFDC code to be executed in
a water quality only simulation mode.
The computational scheme used in the internal eutrophication models employs a fractional step extension
of the same advective and diffusive algorithms used for salinity and temperature, which guarantee
3-5
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
positive constituent concentrations. A novel ordering of the reaction sequence in the reactive source and
sink fractional step allows the linearized reactions to be solved implicitly, further guaranteeing positive
concentrations. The eutrophication models accept an arbitrary number of point and nonpoint source
loadings as well as atmospheric and ground water loadings.
In addition to the internal eutrophication models, the EFDC model can be externally linked to the
WASPS model. In the external linking mode, the EFDC model generates WASPS input files describing
cell geometry and connectivity as well as advective and diffusive transport fields. For estuary simulation,
the transport fields may be intratidally time averaged or intertidally time averaged using the averaging
procedure described by Hamrick (1994).
3.5 Toxic Contaminant Transport and Fate
The EFDC code includes two internal submodels for simulating the transport and fate of toxic
contaminants. A simple, single contaminant submodel can be activated from the master input file. The
simple model accounts for water and suspended sediment phase transport with equilibrium partitioning
and a lumped first order reaction. Contaminant mass per unit area in the sediment bed is also simulated.
The second, more complex, submodel simulates the transport and fate of an arbitrary number of reacting
contaminants in the water and sediment phases of both the water column and sediment bed. In this mode,
the contaminant transport and fate simulation is functionally similar to the WASPS TOXIC model
(Ambrose et al. 1993), with the added flexibility of simulating an arbitrary number of contaminants, and
the improved accuracy of utilizing more complex three-dimensional physical transport fields in a highly
accurate numerical transport scheme. Water-sediment phases interaction may be represented by
equilibrium or nonlinear sorption processes. In this mode, the multilayer sediment bed formulation is
active, with sediment bed water volume and dissolved contaminant mass balances activated to allow
contaminants to reenter the water column by sediment resuspension, pore water expulsion due to
consolidation, and diffusion from the pore water into the water column. The complex contaminant model
activates a subroutine describing reaction processes with appropriate reaction parameters provided by the
toxic reaction processes input file.
3.6 Finfish and Shellfish Transport
The EFDC code includes the capability of simulating the transport and fate of various life stages of
finfish and shellfish. In addition to advection and diffusion by the ambient flow, mortality, predation,
toxicity, and swimming behavior are simulated. Organism age and ambient environment queued vertical
and horizontal swimming and settling is simulated. Environmental queues include light intensity,
temperature, salinity, and tidal phases.
3-6
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
3.7 Near-Field Discharge Dilution and Mixing Zone Analysis
In addition to the far-field transport and fate simulation capability incorporated into the EFDC code's
water quality and toxic contaminant modules, the code includes a near-field discharge dilution and
mixing zone module. The near field model is based on a Lagrangian buoyant jet and plume model (Frick
1984; Lee and Cheung 1990) and allows representation of submerged single and multiple port diffusers
and buoyant surface jets. The near field model provides analysis capabilities similar to CORMIX (Jirka
and Doneker 1991; Jirka and Akar 1991) while offering two distinct advantages. The first advantage is
that a more realistic representation of ambient current and stratification conditions, provided directly by
the EFDC hydrodynamic module, is incorporated into the analysis. The second advantage is that multiple
discharges and multiple near field analysis times may be specified to account for varying ambient current
and stratification conditions. For example, the analysis of 10 discharges under six ambient conditions
each would require 60 executions of CORMIX, while the entire analysis of the 60 situations would be
produced in a single EFDC simulation. The near-field simulation may be executed in two modes. The
first provides virtual source information for representing the discharges in a standard EFDC far-field
transport and fate simulation. In the second mode the near-field and far-field transport are directly
coupled, using a virtual source formulation, to provide simultaneous near and far field transport and fate
simulation.
3.8 Spill Trajectory and Search and Rescue Simulation
In addition to the Eulerian transport equation formulation used for far field analysis and the Lagrangian
jet and plume module used for near field analysis, the EFDC code incorporates a number of Lagrangian
particle transport formulations based on an implicit trilinear interpolation scheme (Bennett and Clites
1987). The first formulation allows release of neutrally buoyant or buoyant drifters at user specified
locations and times. This formulation is useful in simulating spill trajectories, search and rescue
operations, and oceanographic instrument drifters. The second formulation releases drifters in each
three-dimensional model cell at a specified sequence of times and calculates the generalized Lagrangian
mean velocity field (Andrews and Mclntyre 1978) relative to a user-specified averaging interval.
3.9 Wetland, Marsh, and Tidal Flat Simulation Extension
The EFDC model provides a number of enhancements for the simulation of flow and transport in
wetlands, marshes, and tidal flats. The code allows for drying and wetting in shallow areas by a mass
conservative scheme. The drying and wetting formulation is coupled to the mass transport equations in a
manner that prevents negative concentrations of dissolved and suspended materials. A number of
alternatives are in place in the model to simulate general discharge control structures such as weirs,
spillways, culverts, and water surface elevation activated pumps. The effect of submerged and emergent
plants is incorporated into the turbulence closure model and flow resistance formulation. Plant density
3-7
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
and geometric characteristics of individual and composite plants are required as input for the vegetation
resistance formulation. A simple soil moisture model, allowing rainfall infiltration and soil water loss
due to evapotranspiration under dry conditions, is implemented. To represent narrow channels and
canals in wetland, marsh and tidal flat systems, a subgrid scale channel model is implemented. The
subgrid channel model allows a 1-D network in the horizontal channels to be dynamically coupled to the
two-dimensional horizontal grid representing the wetland, marsh, or tidal flat system. Volume and mass
exchanges between 2-D wetland cells and the 1-D channels are accounted for. The channels may
continue to flow when the 2-D wetland cells become dry.
3.10 Nearshore Wave-Induced Currents and Sediment Transport Extensions
The EFDC code includes a number of extensions for simulation of nearshore wave-induced currents and
noncohesive sediment transport. The extensions include a wave-current boundary layer formulation
similar to that of Grant and Madsen (1986); modifications of the hydrodynamic model's momentum
equations to represent wave period averaged Eulerian mean quantities; the inclusion of the
three-dimensional wave-induced radiation or Reynold's stresses in the momentum equations; and
modifications of the velocity fields in the transport equations to include advective transport by the wave-
induced Stake's drift. High frequency surface wave fields are provide by an external wave
refraction-diffraction model or by an internal mild slope equation submodel similar to that of Madsen and
Larsen (1987). The internal refraction-diffraction computation is executed on a refined horizontal grid
coincident with the main model's horizontal grid.
3.11 User Interface
The EFDC modeling package's user interface is based on text input file templates. This choice was
selected in the interest of maintaining model portability across a range of computing platforms and
readily allows the model user to modify input files using most text editing software. The text interface
also allows modification of model files on remote computing systems and in hetrogeneous network
environments. All input files have standard templates available with the EFDC code and in the digital
version of the user's manual. The file templates include extensive built-in documentation and an
explanation of numerical input data quantities. Actual numerical input data are inserted into the text
template in a flexible free format as internally specified in the file templates. Extensive checking of
input files is implemented in the code and diagnostic on screen messages indicate the location and nature
of input file errors. All input files involving dimensional data have unit conversion specifications for the
Meters-Kilograms-Seconds (MKS) international system of units used internally in the model.
3-8
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
3.12 Preprocessing Software
The EFDC modeling package includes a grid generating preprocessor code, GEFDC, which is used to
construct the horizontal model grid, and interpolate bathymetry and initial fields such as water surface
elevation, salinity, to the grid cells. EFDC inputs files specifying the grid geometry and initial fields are
generated by the preprocessor. The preprocessor is capable of generating Cartesian and
curvilinear-orthogonal grids using a number of grid generation schemes (Mobley and Stewart 1980;
Ryskin and Leal 1983; Kang and Leal 1992).
3.13 Program Configuration
The EFDC code exists in only one generic version. A model application is specified entirely by
information in the input files. To minimize memory requirements for specific applications, an executable
file is created by adjusting the appropriate variable array size in the model's parameter file and compiling
the source code. The EFDC model can be configured to execute all or a portion of a model application in
reduced spatial dimension mode including 2-D depth or width averaged and 1-D cross section averaged.
The number of layers used in the 3-D mode or 2-D width averaged mode is readily changed by one line
of model input. Model grid sections specified as 2-D width averaged are allowed to have depth varying
widths to provide representations equivalent to those of 2-D width averaged estuarine and reservoir
models such as CE-QUAL-W2 (Cole and Buchak 1994).
3.14 Run-Time Diagnostics
The EFDC modeling package includes extensive built-in run-time diagnostics that may be activated in
the master input file by the model user. Representative diagnostics include records of maximum CFL
numbers, times and locations of negative depths, a variety of volume and mass balance checks, and
global mass and energy balances. An on screen print of model variables in a specified cell can be
activated during modeling execution. The model generates a number of log files that allow additional
diagnostics of any run-time problems encountered during the set-up of a new application.
3.15 Model Output Options
A wide variety of output options are available for the EFDC model, including (1) specification of output
files for horizontal plane and vertical plane transect plotting of vector and scalar field at a specified time;
(2) the generation of time series of model variables at selected locations and time intervals; (3) grab
sample simulation at specified times and locations; and (4) the specification of least squares analysis of
selected model variables at a defined location over a specified interval. A general three-dimensional
output option allows saving of all major model variables in a compressed-file format at specified times.
A restart file is generated at user-specified intervals during model execution.
3-9
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
3.16 Postprocessing, Graphics, and Visualization
The generic model output files can be readily processed by a number of third party graphics and
visualization software packages, often without the need for intermediate processing (Rennie and Hamrick
1992). The availability of the source code to the user allows the code to be modified for specific output
options. Graphics and visualization software successfully used with EFDC output include: APE, AVS,
IDL, Mathematica, MatLab, NCAR Graphics, PV-Wave, Techplot, Site View, Spyglass Transform and
Slicer, Voxelview, and GrADS. The model developer currently uses Spyglass and Voxelview and a
number of postprocessor applications are available for special image enhancement for these products.
3.17 Documentation
Extensive documentation of the EFDC model is available. Theoretical and computational aspects of the
model are described by Hamrick (1992a). The model user's manual (Hamrick 1996) provides details on
use of the GEFDC preprocessor and set-up of the EFDC input files. Input file templates are also
included. A number of papers describe model applications and capabilities (Hamrick 1992b; Hamrick
1994; Moustafa and Hamrick 1994; Hamrick and Wu 1996; and Wu et al. 1996).
3.18 Computer Requirements
The EFDC modeling system is written in FORTRAN 77. The few nonstandard VAX FORTRAN
language extensions in the code are supported by a wide variety of ANSI standard FORTRAN 77
compilers. The generic or universal source code has been compiled and executed on most UNIX
workstations (DEC Alpha, Hewlett-Packard, IBM RISC6000, Silicon Graphics, Sun and Spare
compatibles) Cray and Convex supercomputers, and PC compatible and Macintosh personal computers.
Absoft, Lahey, and Microsoft compilers are supported on PC compatibles, while Absoft, Language
Systems, and Motorola compilers are supported on Macintosh and compatible systems.
3-10
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4 - EFDC WATER QUALITY MODEL
4.1 Introduction
The central issues in the water quality model are primary production of carbon by algae and
concentration of dissolved oxygen. Primary production provides the energy required by the ecosystem to
function. However, excessive primary production is detrimental since its decomposition in the water and
sediments consumes oxygen. Dissolved oxygen is necessary to support the life functions of higher
organisms and is considered an indicator of the health of estuarine systems. To predict primary
production and dissolved oxygen, a large suite of model state variables is necessary (Table 4-1). The
nitrate state variable in the model represents the sum of nitrate and nitrite nitrogen. The three variables
(salinity, water temperature, and total suspended solids) needed for computation of the above 21 state
variables are provided by the EFDC hydrodynamic model. The interactions among the state variables is
illustrated in Figure 4-1. The kinetic processes included in the EFDC water quality model are mostly
from the Chesapeake Bay three-dimensional water quality model, CE-QUAL-ICM (Cerco and Cole
1994). The kinetic sources and sinks, as well as the external loads for each state variable, are described
in Sections 4.3 to 4.11. The kinetic processes include the exchange effluxes at the sediment-water
interface, including sediment oxygen demand, which are explained in Section 5 (EFDC Sediment Process
Model) of this report. The description of the EFDC water column water quality model in this section is
from Park etal. (1995).
Table 4-1. EFDC model water quality state variables
(1) cyanobacteria (12) labile particulate organic nitrogen
(2) diatom algae (13) dissolved organic nitrogen
(3) green algae (14) ammonia nitrogen
(4) refractory particulate organic carbon (15) nitrate nitrogen
(5) labile particulate organic carbon (16) particulate biogenic silica
(6) dissolved organic carbon (17) dissolved available silica
(7) refractory particulate organic phosphorus (18) chemical oxygen demand
(8) labile particulate organic phosphorus (19) dissolved oxygen
(9) dissolved organic phosphorus (20) total active metal
(10) total phosphate (21) fecal coliform bacteria
(11) refractory particulate organic nitrogen (22) macroalgae
4-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
* TSS from hydrodynamic
model
Figure 4-1. Schematic diagram for the EFDC water column water quality model
4.1.1 Algae
Algae are grouped into four model classes: cyanobacteria, diatoms, greens, and macroalgae. The
grouping is based upon the distinctive characteristics of each class and upon the significant role the
characteristics play in the ecosystem. Cyanobacteria, commonly called blue-green algae, are
characterized by their abundance (as picoplankton) in saline water and by their bloom-forming char-
acteristics in fresh water. Cyanobacteria are unique in that some species fix atmospheric nitrogen,
although nitrogen fixers are not believed to be predominant in many river systems. Diatoms are distin-
guished by their requirement of silica as a nutrient to form cell walls. Diatoms are large algae
characterized by high settling velocities. Settling of spring diatom blooms to the sediments may be a
significant source of carbon for sediment oxygen demand. Algae that do not fall into the preceding two
groups are lumped into the heading of green algae. Green algae settle at a rate intermediate between
4-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
cyanobacteria and diatoms and are subject to greater grazing pressure than cyanobacteria. Macroalgae
are almost always attached to a stable substrate and are therefore most abundant in the areas of harbors
and near shore. The waters in many stream systems are characterized by various rooted macrophytes and
periphyton. All species of macroalgae in this study have been lumped into a single class of macroalgae.
Because of their attachment to the substrate, they are limited to growing in the bottom water-column
layer and are not subject to physical transport.
4.1.2 Organic Carbon
Three organic carbon state variables are considered: dissolved, labile particulate, and refractory
particulate. Labile and refractory distinctions are based upon the time scale of decomposition. Labile
organic carbon decomposes on a time scale of days to weeks whereas refractory organic carbon requires
more time. Labile organic carbon decomposes rapidly in the water column or the sediments. Refractory
organic carbon decomposes slowly, primarily in the sediments, and may contribute to sediment oxygen
demand years after deposition.
4.1.3 Nitrogen
Nitrogen is first divided into organic and mineral fractions. Organic nitrogen state variables are dissolved
organic nitrogen, labile particulate organic nitrogen, and refractory particulate organic nitrogen. Two
mineral nitrogen forms are considered: ammonium and nitrate. Both are utilized to satisfy algal nutrient
requirements, although ammonium is preferred from thermodynamic considerations. The primary reason
for distinguishing the two is that ammonium is oxidized by nitrifying bacteria into nitrate. This oxidation
can be a significant sink of oxygen in the water column and sediments. An intermediate in the complete
oxidation of ammonium, nitrite, also exists. Nitrite concentrations are usually much less than nitrate, and
for modeling purposes, nitrite is combined with nitrate. Hence the nitrate state variable actually
represents the sum of nitrate plus nitrite.
4.1.4 Phosphorus
As with carbon and nitrogen, organic phosphorus is considered in three states: dissolved, labile
particulate, and refractory particulate. Only a single mineral form, total phosphate, is considered. Total
phosphate exists as several states within the model ecosystem: dissolved phosphate, phosphate sorbed to
inorganic solids, and phosphate incorporated in algal cells. Equilibrium partition coefficients are used to
distribute the total among the three states.
4.1.5 Silica
Silica is divided into two state variables: available silica and particulate biogenic silica. Available silica
is primarily dissolved and can be utilized by diatoms. Particulate biogenic silica cannot be utilized. In the
4-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
model, participate biogenic silica is produced through diatom mortality. Particulate biogenic silica
undergoes dissolution to available silica or else settles to the bottom sediments.
4.1.6 Chemical Oxygen Demand
In the context of this study, chemical oxygen demand is the concentration of reduced substances that are
oxidizable by inorganic means. The primary component of chemical oxygen demand is sulfide released
from sediments. Oxidation of sulfide to sulfate may remove substantial quantities of dissolved oxygen
from the water column.
4.1.7 Dissolved Oxygen
Dissolved oxygen is required for the existence of higher life forms. Oxygen availability determines the
distribution of organisms and the flows of energy and nutrients in an ecosystem. Dissolved oxygen is a
central component of the water quality model.
4.1.8 Total Active Metal
Both phosphate and dissolved silica sorb to inorganic solids, primarily iron and manganese. Sorption and
subsequent settling is one pathway for removal of phosphate and silica from the water column.
Consequently, the concentration and transport of iron and manganese are represented in the model.
Limited data do not allow a complete treatment of iron and manganese chemistry, however. Rather, a
single-state variable, total active metal, is defined as the total concentration of metals that are active in
phosphate and silica transport. Total active metal is partitioned between particulate and dissolved phases
by an oxygen-dependent partition coefficient.
4.1.9 Salinity
Salinity is a conservative tracer that provides verification of the transport component of the model and
facilitates examination of conservation of mass. Salinity also influences the dissolved oxygen saturation
concentration and is used in the determination of kinetics constants that differ in saline and fresh water.
4.1.10 Temperature
Temperature is a primary determinant of the rate of biochemical reactions. Reaction rates increase as a
function of temperature, although extreme temperatures result in the mortality of organisms.
4.2 Conservation of Mass Equation
The governing mass-balance equation for each of the water quality state variables may be expressed as:
4-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
+ d(vQ + 8(wC) =
dt dx dy dz
K .KtK
dx( * dx) dy( y dy) dz( * dz) c '
C = concentration of a water quality state variable
u, v, w = velocity components in the x-, y-, and z-directions, respectively
Kx, Ky, Kz = turbulent diffusivities in the x-, y-, and z-directions, respectively
Sc = internal and external sources and sinks per unit volume.
The last three terms on the left-hand side (LHS) of Eq. 4-1 account for the advective transport, and the
first three terms on the right-hand side (RHS) of Eq. 4-1 account for the diffusive transport. These six
terms for physical transport are analogous to, and thus the numerical method of solution is the same as,
those in the mass-balance equation for salinity in the hydrodynamic model (Hamrick 1992a). The last
term in Eq. 4-1 represents the kinetic processes and external loads for each of the state variables. The
present model solves Eq. 4-1 after decoupling the kinetic terms from the physical transport terms. The
solution scheme for both the physical transport (Hamrick 1 992a) and the kinetic equations is second-
order accurate.
The governing mass-balance equation for water quality state variables (Eq. 4-1) consists of
physical transport, advective and diffusive, and kinetic processes. When solving Eq. 4-1, the kinetic
terms are decoupled from the physical transport terms. The mass-balance equation for physical transport
only, which takes the same form as the salt-balance equation, is:
dC + d(uQ + 8(vQ + 8(wQ = B K 8C + _d K dC + d_ K
dt dx dy dz dx( xdx) dy( y dy ) dz( z dz ) ^ '
The equation for kinetic processes only, which will be referred to as the kinetic equation, is:
dC „
(4'3)
which may be expressed as:
dc
dt
= K-C + R (4-4)
where K is kinetic rate (time"1) and R is source/sink term (mass volume"1 time"1). Equation 4-4 is
obtained by linearizing some terms in the kinetic equations, mostly Monod type expressions. Hence, K
4-5
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
and R are known values in Eq. 4-4. Equation 4-2 is identical to, and thus its numerical method of
solution is the same as, the mass-balance equation for salinity (Hamrick 1992a).
The remainder of this chapter details the kinetics portion of the mass-conservation equation for each state
variable. Parameters are defined where they first appear. All parameters are listed, in alphabetical order,
in an appendix. For consistency with reported rate coefficients, kinetics are detailed using a temporal
dimension of days. Within the CE-QUAL-ICM computer code, kinetics sources and sinks are converted
to a dimension of seconds before employment in the mass-conservation equation.
4.3 Algae
Algae, which occupies a central role in the model (Figure 4-1), are grouped into three model state
variables: cyanobacteria (blue-green algae), diatoms, and green algae. The subscript, x, is used to denote
four algal groups: c for cyanobacteria, d for diatoms, g for green algae, and m for macroalgae. Sources
and sinks included in the model are
• growth (production)
• basal metabolism
• predation
• settling
• external loads
Equations describing these processes are largely the same for the four algal groups with differences in the
values of parameters in the equations. The kinetic equation describing these processes is:
dB a WB
~^ = (PX- BMx - PRx)Bx + ±(WS,-B} + -jl (4-5)
Bx = algal biomass of algal group x (g C m"3)
t = time (day)
Px = production rate of algal group x (day"1)
BMX = basal metabolism rate of algal group x (day"1)
PRx = predation rate of algal group x (day"1)
WSX = settling velocity of algal group x (m day"1)
WBX = external loads of algal group x (g C day"1)
V = cell volume (m3).
The model simulates the total biomass of the macroalgae rather than the size of the macroalgae;
therefore, they can be treated as other groups of algae. Since macroalgae attach to the bottom, they are
limited to growing in the bottom layer only and are not be transported through water movement.
4-6
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.3.1 Production (Algal Growth)
Algal growth depends on nutrient availability, ambient light, and temperature. The effects of these
processes are considered to be multiplicative:
Px = PM,/1(JV)-/2(/)-/3(7) (4-6)
PMX = maximum growth rate under optimal conditions for algal group x (day"1)
fj(N) = effect of suboptimal nutrient concentration (0 < f, < 1)
f2(I) = effect of suboptimal light intensity (0 < f2 < 1)
f3(T) = effect of suboptimal temperature (0 < f, < 1).
The freshwater cyanobacteria may undergo rapid mortality in salt water, e.g., freshwater organisms in the
Potomac River (Thomann et al. 1985). For the freshwater organisms, the increased mortality may be
included in the model by retaining the salinity toxicity term in the growth equation for cyanobacteria:
Pc = PMc -f,(N) -/2(7) -/3(7) -/4(5) (4-7)
f4(S) = effect of salinity on cyanobacteria growth (0 < f4 < 1).
Activation of the salinity toxicity term, f4 (S), is an option in the source code.
4.3.2 Effect of Nutrients on Algal Growth
Using Liebig's "law of the minimum" (Odum 1971) that growth is determined by the nutrient in least
supply, the nutrient limitation for growth of cyanobacteria and green algae is expressed as:
,,,„ . . NH4 + NO3 PO4d
f,(N) - minimum
1
,
KHNr + NH4 + NOB KHPr + PO4d
X X
NH4 = ammonium nitrogen concentration (g N m"3)
NO 3 = nitrate nitrogen concentration (g N m"3)
KHNX = half-saturation constant for nitrogen uptake for algal group x (g N m"3)
PO4d = dissolved phosphate phosphorus concentration (g P m"3)
KHPX = half-saturation constant for phosphorus uptake for algal group x (g P m"3).
Some cyanobacteria, e.g., Anabaena, can fix nitrogen from atmosphere and thus are not limited by
nitrogen. Hence, Eq. 4-8 is not applicable to the growth of nitrogen fixers.
Since diatoms require silica as well as nitrogen and phosphorus for growth, the nutrient limitation for
diatoms is expressed as:
4-7
-------
= minimum
Model Report for Christina River Basin, Nutrient and DO TMDL
NH4 + NO3 PO4d SAd
KHNd + NH4 + NO3 KHPd + PO4d KHS + SAd\
(4-9)
SAd = concentration of dissolved available silica (g Si m"3)
KHS = half-saturation constant for silica uptake for diatoms (g Si m"3).
4.3.3 Effect of Light on Algal Growth
The daily and vertically integrated form of Steele's equation is:
2.718-FD
-£ - 7-
KCSS-&Z
.
«T\
I
- e I (4-10)
(4-11)
(4-12)
FD = fractional daylength (0 < FD < 1)
Kess = total light extinction coefficient (m"1)
Az = layer thickness (m)
I0 = daily total light intensity at water surface (langleys day"1)
(Is)x = optimal light intensity for algal group x (langleys day"1)
HT = depth from the free surface to the top of the layer (m).
Light extinction in the water column consists of three fractions in the model: a background value
dependent on water color, extinction due to suspended particles, and extinction due to light absorption by
ambient chlorophyll:
V- Bx
Kess = Keb + KeTSS-TSS + Kecu- ^ (——-) (4-13)
x=c,d,g CCm..
x
Keb = background light extinction (m"1)
KeTSS = light extinction coefficient for total suspended solid (m"1 per g m"3)
TSS = total suspended solid concentration (g m"3) provided from the hydrodynamic model
Kechl = light extinction coefficient for chlorophyll 'a' (m"1 per mg Chi m"3)
CChlx = carbon-to-chlorophyll ratio in algal group x (g C per mg Chi).
4-8
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Since macroalgae only attach to the bottom, they are not included in computation of the light extinction
Self shading is not considered for macroalgae for the present model. For a model application that does
not simulate TSS, the KeTSS term may be set to zero and Keb may be estimated to include light extinction
due to suspended solid.
Optimal light intensity (Is) for photosynthesis depends on algal taxonomy, duration of exposure,
temperature, nutritional status, and previous acclimation. Variations in Is are largely due to adaptations
by algae intended to maximize production in a variable environment. Steel (1962) noted the result of
adaptations is that optimal intensity is a consistent fraction (approximately 50%) of daily intensity.
Kremer and Nixon (1981) reported an analogous finding that maximum algal growth occurs at a constant
depth (approximately 1 m) in the water column. Their approach is adopted so that optimal intensity is
expressed as:
(I)x - maximum {(/0)^ e ~^'(ZV* , (/,)fflin} (4-14)
(Dopt)x = depth of maximum algal growth for algal group x (m)
(Uavg = adjusted surface light intensity (langleys day"1).
A minimum, (Is)mm, in Eq. 4-14 is specified so that algae do not thrive at extremely low light levels. The
time required for algae to adapt to changes in light intensity is recognized by estimating (Is)x based on a
time-weighted average of daily light intensity:
(Uv. = CI-I0 + db-I, + CI-I2 (4-15)
Ij = daily light intensity 1 day preceding model day (langleys day"1)
I2 = daily light intensity 2 days preceding model day (langleys day"1)
CIa, CIb, CIC = weighting factors for I0, Ij and I2, respectively: CIa + CIb + CIC = 1.
4.3.4 Effect of Temperature on Algal Growth
A Gaussian probability curve is used to represent temperature dependency of algal growth:
/3(7) = exp(-KTGlx[T - TMJ2) if T < TMx
= exp(-KTG2x[TMx - T]2) if T > TMx (4-16)
T = temperature (°C) provided from the hydrodynamic model
TMX = optimal temperature for algal growth for algal group x (°C)
KTG1X = effect of temperature below TMX on growth for algal group x (°C~2)
KTG2X = effect of temperature above TMX on growth for algal group x (°C~2).
4-9
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.3.5 Effect of Salinity on Growth of Freshwater Cyanobacteria
The growth of freshwater cyanobacteria in salt water is limited by:
STOX2
STOX = salinity at which Microcystis growth is halved (ppt)
S = salinity in water column (ppt) provided from the hydrodynamic model.
4.3.6 Algal Basal Metabolism
Algal biomass in the present model decreases through basal metabolism (respiration and excretion) and
predation. Basal metabolism in the present model is the sum of all internal processes that decrease algal
biomass and consists of two parts; respiration and excretion. In basal metabolism, algal matter (carbon,
nitrogen, phosphorus, and silica) is returned to organic and inorganic pools in the environment, mainly to
dissolved organic and inorganic matter. Respiration, which may be viewed as a reversal of production,
consumes dissolved oxygen. Basal metabolism is considered to be an exponentially increasing function
of temperature:
BMx = BMRx-exp(KTBx[T - TRJ) (4-18)
BMR,, = basal metabolism rate at TRX for algal group x (day"1)
KTBX = effect of temperature on metabolism for algal group x (''C"1)
TR,, = reference temperature for basal metabolism for algal group x (°C).
4.3.7 Algal Predation
The present model does not include zooplankton. Instead, a constant rate is specified for algal predation,
which implicitly assumes zooplankton biomass is a constant fraction of algal biomass. An equation
similar to that for basal metabolism (Eq. 4-18) is used for predation:
PRX - PRRx-exp(KTBx[T - TRx]) (4-19)
PRRx = predation rate at TR,, for algal group x (day"1).
The difference between predation and basal metabolism lies in the distribution of the end products of the
two processes. In predation, algal matter (carbon, nitrogen, phosphorus, and silica) is returned to organic
and inorganic pools in the environment, mainly to particulate organic matter. The predation for
macroalgae is a lumped parameter that includes losses due to grazing, frond breakage, and other losses.
This implicitly assumes that the losses are a fraction of the biomass.
4-10
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.3.8 Algal Settling
Settling velocities for four algal groups, WSC, WSd , WSg, and WSm, are specified as an input. Seasonal
variations in settling velocity of diatoms can be accounted for by specifying time-varying WSd.
4.4 Organic Carbon
The present model has three state variables for organic carbon: refractory particulate, labile particulate,
and dissolved.
4.4.1 Particulate Organic Carbon
Labile and refractory distinctions are based on the time scale of decomposition. Labile particulate
organic carbon with a decomposition time scale of days to weeks decomposes rapidly in the water
column or in the sediments. Refractory particulate organic carbon with a longer-than-weeks
decomposition time scale decomposes slowly, primarily in the sediments, and may contribute to sediment
oxygen demand years after decomposition. For labile and refractory particulate organic carbon, sources
and sinks included in the model are (Fig. 4-1):
• algal predation
• dissolution to dissolved organic carbon
• settling
• external loads.
The governing equations for refractory and labile particulate organic carbons are:
FCRP-PR-Bx - KSPOC-SPOC + (WS^-RPOQ + - (4-20)
x=c,d,g,m OZ V
FCLP-PR-Bx - KLPOC-LPOC +
O x=c,d,g,m OZ
RPOC = concentration of refractory particulate organic carbon (g C m"3)
LPOC = concentration of labile particulate organic carbon (g C m"3)
FCRP = fraction of predated carbon produced as refractory particulate organic carbon
FCLP = fraction of predated carbon produced as labile particulate organic carbon
= dissolution rate of refractory particulate organic carbon (day"1)
= dissolution rate of labile particulate organic carbon (day"1)
WSRP = settling velocity of refractory particulate organic matter (m day"1)
WSLP = settling velocity of labile particulate organic matter (m day"1)
WRPOC = external loads of refractory particulate organic carbon (g C day"1)
WLPOC = external loads of labile particulate organic carbon (g C day"1).
4-11
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.4.2 Dissolved Organic Carbon
Sources and sinks for dissolved organic carbon included in the model are (Fig. 4-1):
• algal excretion (exudation) and predation
• dissolution from refractory and labile particulate organic carbon
• heterotrophic respiration of dissolved organic carbon (decomposition)
• denitrification
• external loads
The kinetic equation describing these processes is:
dpoc
dt
- E
x=c,d,g,m
FCDx + (1 - FCDx)
KHR
KHR + DO
x
BM + FCDP-PR
•B
KLPOC-LPOC - KHR-DOC - Denit-DOC + WDOC (4.22)
DOC = concentration of dissolved organic carbon (g C m"3)
FCDX = fraction of basal metabolism exuded as dissolved organic carbon at infinite dissolved oxygen
concentration for algal group x
KHR,, = half-saturation constant of dissolved oxygen for algal dissolved organic carbon excretion for
group x (g O2 m"3)
DO = dissolved oxygen concentration (g O2 m"3)
FCDP = fraction of predated carbon produced as dissolved organic carbon
KHR = heterotrophic respiration rate of dissolved organic carbon (day"1)
Denit = denitrification rate (day"1) given in Eq. 4-34
WDOC = external loads of dissolved organic carbon (g C day"1).
The remainder of this section explains each term in Equations 4-20 to 4-22.
4.4.3 Effect of Algae on Organic Carbon
The terms within summation (£) in Equations 4-20 to 4-22 account for the effects of algae on organic
carbon through basal metabolism and predation.
4.4.3.1 Basal metabolism. Basal metabolism, consisting of respiration and excretion, returns algal
matter (carbon, nitrogen, phosphorus, and silica) back to the environment. Loss of algal biomass through
basal metabolism is (Eq. 4-18):
4-12
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
which indicates that the total loss of algal biomass due to basal metabolism is independent of ambient
dissolved oxygen concentration. In this model, it is assumed that the distribution of total loss between
respiration and excretion is constant as long as there is sufficient dissolved oxygen for algae to respire.
Under that condition, the losses by respiration and excretion may be written as:
(1 - FCDx)-BMxBx due to respiration (4-24)
FCDxBMxBx due to excretion (4-25)
where FCDX is a constant of value between 0 and 1. Algae cannot respire in the absence of oxygen,
however. Although the total loss of algal biomass due to basal metabolism is oxygen-independent
(Eq. 4-23), the distribution of total loss between respiration and excretion is oxygen-dependent. When
oxygen level is high, respiration is a large fraction of the total. As dissolved oxygen becomes scarce,
excretion becomes dominant. Thus, Eq. 4-24 represents the loss by respiration only at high oxygen
levels. In general, Eq. 4-24 can be decomposed into two fractions as a function of dissolved oxygen
availability:
(1 - FCD^ KHR + DOBMxB* due to respiration (4-26)
KHRx
(1 - FCDJ —— —BMxBx due to excretion (4-27)
KHRx + DO
Equation 4-26 represents the loss of algal biomass by respiration, and Eq. 4-27 represents additional
excretion due to insufficient dissolved oxygen concentration. The parameter KHR,,, which is defined as
the half-saturation constant of dissolved oxygen for algal dissolved organic carbon excretion in Eq. 4-22,
also be defined as the half-saturation constant of dissolved oxygen for algal respiration in Eq. 4-26.
Combining Equations 4-25 and 4-27, the total loss due to excretion is:
can
KHR
FCDx + (1 - FCDJ
KHR + DO
BM-Br (4-28)
Equations 4-26 and 4-28 combine to give the total loss of algal biomass due to basal metabolism, BMX-BX
(Eq. 4-23). The definition of FCDX in Eq. 4-22 becomes apparent in Eq. 4-28; i.e., fraction of basal
metabolism exuded as dissolved organic carbon at infinite dissolved oxygen concentration. At zero
oxygen level, 100% of total loss due to basal metabolism is by excretion regardless of FCDX. The end
carbon product of respiration is primarily carbon dioxide, an inorganic form not considered in the present
model, while the end carbon product of excretion is primarily dissolved organic carbon. Therefore, Eq.
4-28, that appears in Eq. 4-22, represents the contribution of excretion to dissolved organic carbon, and
4-13
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
there is no source term for participate organic carbon from algal basal metabolism in Equations 4-20 and
4-21.
4.4.3.2 Predation. Algae produce organic carbon through the effects of predation. Zooplankton take up
and redistribute algal carbon through grazing, assimilation, respiration, and excretion. Since zooplankton
are not included in the model, routing of algal carbon through zooplankton predation is simulated by
empirical distribution coefficients in Equations 4-20 to 4-22; FCRP, FCLP, and FCDP. The sum of these
three predation fractions should be unity.
4.4.4 Heterotrophic Respiration and Dissolution
The second term on the RHS of Equations 4-20 and 4-21 represents dissolution of particulate to
dissolved organic carbon and the third term in the second line of Eq. 4-22 represents heterotrophic
respiration of dissolved organic carbon. The oxic heterotrophic respiration is a function of dissolved
oxygen: the lower the dissolved oxygen, the smaller the respiration term becomes. Heterotrophic
respiration rate, therefore, is expressed using a Monod function of dissolved oxygen:
r - D0 r
HR = KHORDO + DO DOC (4'29)
KHORDO = oxic respiration half-saturation constant for dissolved oxygen (g O2 m"3)
KDOC = heterotrophic respiration rate of dissolved organic carbon at infinite dissolved oxygen
concentration (day"1).
Dissolution and heterotrophic respiration rates depend on the availability of carbonaceous substrate and
on heterotrophic activity. Algae produce labile carbon that fuels heterotrophic activity: dissolution and
heterotrophic respiration do not require the presence of algae though, and may be fueled entirely by
external carbon inputs. In the model, algal biomass, as a surrogate for heterotrophic activity, is
incorporated into formulations of dissolution and heterotrophic respiration rates. Formulations of these
rates require specification of algal -dependent and algal-independent rates:
KRPOC = (KRC + KRCalg x HDR ~ HDR (4-30)
KLPOC - (KLC +
KDOC - (KDC + KDCaig E - ~ wi (4-32)
x=c,d,g
KRC = minimum dissolution rate of refractory particulate organic carbon (day"1)
4-14
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
KLC = minimum dissolution rate of labile particulate organic carbon (day"1)
KDC = minimum respiration rate of dissolved organic carbon (day"1)
KRCaig, KLCalg = constants that relate dissolution of refractory and labile particulate organic carbon,
respectively, to algal biomass (day"1 per g C m"3)
KDCalg = constant that relates respiration to algal biomass (day"1 per g C m"3)
KTHDR = effect of temperature on hydrolysis of particulate organic matter (''C"1)
TRHDR = reference temperature for hydrolysis of particulate organic matter (°C)
= effect of temperature on mineralization of dissolved organic matter (''C"1)
= reference temperature for mineralization of dissolved organic matter (°C).
Equations 4-30 to 4-32 have exponential functions that relate rates to temperature.
In the present model, the term "hydrolysis" is defined as the process by which particulate organic matter
is converted to dissolved organic form, and thus includes both dissolution of particulate carbon and
hydrolysis of particulate phosphorus and nitrogen. Therefore, the parameters, KTHDR and TRHDR, are also
used for the temperature effects on hydrolysis of particulate phosphorus (Equations 4-28 and 4-29) and
nitrogen (Equations 4-54 and 4-55). The term "mineralization" is defined as the process by which
dissolved organic matter is converted to dissolved inorganic form, and thus includes both heterotrophic
respiration of dissolved organic carbon and mineralization of dissolved organic phosphorus and nitrogen.
Therefore, the parameters, KT^L and TR^L, are also used for the temperature effects on mineralization
of dissolved phosphorus (Eq. 4-46) and nitrogen (Eq. 4-56).
4.4.5 Effect of Denitrification on Dissolved Organic Carbon
As oxygen is depleted from natural systems, organic matter is oxidized by the reduction of alternate
electron acceptors. Thermodynamically, the first alternate acceptor reduced in the absence of oxygen is
nitrate. The reduction of nitrate by a large number of heterotrophic anaerobes is referred to as
denitrification, and the stoichiometry of this reaction is (Stumm and Morgan 1981):
4NO3~ + 4/T + 5CH2O - 2N2 + 7H2O + 5CO2 (4-33)
The last term in Eq. 4-22 accounts for the effect of denitrification on dissolved organic carbon. The
kinetics of denitrification in the model are first-order:
KHORnn NQ3
Denit = — -^=- AANOX-Knnr (4.34)
KHORDO + DO KHDNN + NO3 DOC ^ '
KHDNN = denitrification half-saturation constant for nitrate (g N m"3)
AANOX = ratio of denitrification rate to oxic dissolved organic carbon respiration rate.
4-15
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
In Eq. 4-34, the dissolved organic carbon respiration rate, KDOC, is modified so that significant
decomposition via denitrification occurs only when nitrate is freely available and dissolved oxygen is
depleted. The ratio, AANOX, makes the anoxic respiration slower than oxic respiration. Note that KDOC,
defined in Eq. 4-32, includes the temperature effect on denitrification.
4.5 Phosphorus
The present model has four state variables for phosphorus: three organic forms (refractory particulate,
labile particulate, and dissolved) and one inorganic form (total phosphate).
4.5.1 Particulate Organic Phosphorus
For refractory and labile particulate organic phosphorus, sources and sinks included in the model are
(Fig. 4-1):
• algal basal metabolism and predation
• dissolution to dissolved organic phosphorus
• settling
• external loads.
The kinetic equations for refractory and labile particulate organic phosphorus are:
dRPOP
(FPRXBMX + FPRP-PRx)APC-Bx -
CT x-c,d,g,m
+ -^(WSgp-RPOF) + WRpOP (4_35)
dLPOP = ^ (FPLXBMX + FPLP-PRx)APC-Bx - KLpQp-LPOP
(Jl x=c,d,g,m
WLPOP
(4~36)
RPOP = concentration of refractory particulate organic phosphorus (g P m"3)
LPOP = concentration of labile particulate organic phosphorus (g P m"3)
FPR,, = fraction of metabolized phosphorus by algal group x produced as refractory particulate organic
phosphorus
FPLX = fraction of metabolized phosphorus by algal group x produced as labile particulate organic
phosphorus
FPRP = fraction of predated phosphorus produced as refractory particulate organic phosphorus
FPLP = fraction of predated phosphorus produced as labile particulate organic phosphorus
APC = mean algal phosphorus-to-carbon ratio for all algal groups (g P per g C)
4-16
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
= hydrolysis rate of refractory particulate organic phosphorus (day"1)
KLPOP = hydrolysis rate of labile particulate organic phosphorus (day"1)
WRPOP = external loads of refractory particulate organic phosphorus (g P day"1)
WLPOP = external loads of labile particulate organic phosphorus (g P day"1).
4.5.2 Dissolved Organic Phosphorus
Sources and sinks for dissolved organic phosphorus included in the model are (Fig. 4-1):
• algal basal metabolism and predation
• dissolution from refractory and labile particulate organic phosphorus
• mineralization to phosphate phosphorus
• external loads.
The kinetic equation describing these processes is:
dDOP
(FPDx-BMx + FPDP-PRx)APC-Bx
Ot x=c,d,g,m
KLpQpLPOP - KDQpDOP + (4-37)
DOP = concentration of dissolved organic phosphorus (g P m"3)
FPDX = fraction of metabolized phosphorus by algal group x produced as dissolved organic phosphorus
FPDP = fraction of predated phosphorus produced as dissolved organic phosphorus
KDOP = mineralization rate of dissolved organic phosphorus (day"1)
WDOP = external loads of dissolved organic phosphorus (g P day"1).
4.5.3 Total Phosphate
For total phosphate that includes both dissolved and sorbed phosphate (Section 4.5.4), sources and sinks
included in the model are (Fig. 4-1):
• algal basal metabolism, predation, and uptake
• mineralization from dissolved organic phosphorus
• settling of sorbed phosphate
• sediment-water exchange of dissolved phosphate for the bottom layer only
• external loads.
The kinetic equation describing these processes is:
+ FPIP-PRx - Px}APC-Bx + KDOpDOP
x=c,d,g,m
4-17
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
BFPO4d WPO4t
z ^ + —— (4.38)
PO4t = total phosphate (g P m"3) = PO4d + PO4p (4-39)
PO4d = dissolved phosphate (g P m"3)
PO4p = particulate (sorbed) phosphate (g P m"3)
FPIX = fraction of metabolized phosphorus by algal group x produced as inorganic phosphorus
FPIP = fraction of predated phosphorus produced as inorganic phosphorus
WSTSS = settling velocity of suspended solid (m day"1), provided by the hydrodynamic model
BFPO4d = sediment-water exchange flux of phosphate (g P m"2 day"1), applied to the bottom layer only
WPO4t = external loads of total phosphate (g P day"1).
In Eq. 4-38, if total active metal is chosen as a measure of sorption site, the settling velocity of total
suspended solid, WSTSS, is replaced by that of particulate metal, WSS (Sections 4.5.4 and 4.10). The
remainder of this section explains each term in Equations 4-35 to 4-38, except BFPO4d (benthic flux of
dissolved orthophosphate), which is described in Chapter 5.
4.5.4 Total Phosphate System
Suspended and bottom sediment particles (clay, silt, and metal hydroxides) adsorb and desorb phosphate
in river and estuarine waters. This adsorption-desorption process has been suggested to buffer phosphate
concentration in water column and to enhance the transport of phosphate away from its external sources
(Carritt and Goodgal 1954; Froelich 1988; Lebo 1991). To ease the computational complication due to
the adsorption-desorption of phosphate, dissolved and sorbed phosphate are treated and transported as a
single state variable. Therefore, the model phosphate state variable, total phosphate, is defined as the
sum of dissolved and sorbed phosphate (Eq. 4-39), and the concentrations for each fraction are
determined by equilibrium partitioning of their sum.
In CE-QUAL-ICM, sorption of phosphate to particulate species of metals including iron and manganese
was considered based on a phenomenon observed in the monitoring data from the mainstem of the
Chesapeake Bay: phosphate was rapidly depleted from anoxic bottom waters during the autumn
reaeration event (Cerco and Cole 1993). Their hypothesis was that reaeration of bottom waters caused
dissolved iron and manganese to precipitate, and phosphate sorbed to newly formed metal particles and
rapidly settled to the bottom. One state variable, total active metal, in CE-QUAL-ICM was defined as the
sum of all metals that act as sorption sites, and the total active metal was partitioned into particulate and
dissolved fractions via an equilibrium partitioning coefficient (Section 4.10). Then phosphate was
assumed to sorb to only the particulate fraction of the total active metal.
4-18
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
In the treatment of phosphate sorption in CE-QUAL-ICM, the particulate fraction of metal hydroxides
was emphasized as a sorption site in bottom waters under anoxic conditions. Phosphorus is a highly
particle-reactive element, and phosphate in solution reacts quickly with a wide variety of surfaces, being
taken up by and released from particles (Froelich 1988). The present model has two options, total
suspended solid and total active metal, as a measure of a sorption site for phosphate, and dissolved and
sorbed fractions are determined by equilibrium partitioning of their sum as a function of total suspended
solid or total active metal concentration:
K -TSS K -TAMo
PO4p = ^ PO4t or PO4p = ^ ^PO4t (4.40)
1 , V T'OO 1 , V T AK jf~ \ '
1 + KP04p'TSS l + KP04p'TAMP
PO4d = PO4t nr PO4d = PO4t
1 + KpQ4pTSS 1 + KpQ4pTAMp
= PO4t - PO4p (4-41)
Kpo4P = empirical coefficient relating phosphate sorption to total suspended solid (per g m"3) or
particulate total active metal (per mol m"3) concentration
TAMp = particulate total active metal (mol m"3).
Dividing Eq. 4-40 by Eq. 4-41 gives:
-rr J- \-ftLJ J. -rr -i \-ftLJ J.
P04p = ~P04~dJsS °r P04p = ~P04~dTAMp (4'42)
where the meaning of KP04p becomes apparent, i.e., the ratio of sorbed to dissolved phosphate per unit
concentration of total suspended solid or particulate total active metal (i.e., per unit sorption site
available).
4.5.5 Algal Phosphorus-to-Carbon Ratio (APC)
Algal biomass is quantified in units of carbon per volume of water. In order to express the effects of
algal biomass on phosphorus and nitrogen, the ratios of phosphorus-to-carbon and nitrogen-to-carbon in
algal biomass must be specified. Although global mean values of these ratios are well known (Redfield
et al. 1963), algal composition varies especially as a function of nutrient availability. As phosphorus and
nitrogen become scarce, algae adjust their composition so that smaller quantities of these vital nutrients
are required to produce carbonaceous biomass (DiToro 1980; Parsons et al. 1984). Examining the field
data from the surface of upper Chesapeake Bay, Cerco and Cole (1993) showed that the variation of
nitrogen-to-carbon stoichiometry was small and thus used a constant algal nitrogen-to-carbon ratio,
ANCX. Large variations, however, were observed for algal phosphorus-to-carbon ratio indicating the
4-19
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
adaptation of algae to ambient phosphorus concentration (Cerco and Cole 1993): algal phosphorus
content is high when ambient phosphorus is abundant and is low when ambient phosphorus is scarce.
Thus, a variable algal phosphorus-to-carbon ratio, APC, is used in model formulation. A mean ratio for
all algal groups, APC, is described by an empirical approximation to the trend observed in field data
(Cerco & Cole 1994):
- (CPprml + CPprm2-^[-CPprm3-P04d\Y (4-43)
CPpnni = minimum carbon-to-phosphorus ratio (g C per g P)
CPpnn2 = difference between minimum and maximum carbon-to-phosphorus ratio (g C per g P)
CPprms = effect of dissolved phosphate concentration on carbon-to-phosphorus ratio (per g P m"3).
4.5.6 Effect of Algae on Phosphorus
The terms within summation (£) in Equations 4-35 to 4-38 account for the effects of algae on
phosphorus. Both basal metabolism (respiration and excretion) and predation are considered, and thus
formulated, to contribute to organic and phosphate phosphorus. That is, the total loss by basal
metabolism (BMX-BX in Eq. 4-5) is distributed using distribution coefficients; FPRX, FPLX, FPDX, and
FPIX. The total loss by predation (PRX-BX in Eq. 4-5), is also distributed using distribution coefficients;
FPRP, FPLP, FPDP, and FPIP. The sum of four distribution coefficients for basal metabolism should be
unity, and so is that for predation. Algae take up dissolved phosphate for growth, and algae uptake of
phosphate is represented by (- £ PX-APGBX) in Eq. 4-38.
4.5.7 Mineralization and Hydrolysis
The third term on the RHS of Equations 4-35 and 4-36 represents hydrolysis of particulate organic
phosphorus, and the last term in Eq. 3-7 represents mineralization of dissolved organic phosphorus.
Mineralization of organic phosphorus is mediated by the release of nucleotidase and phosphatase
enzymes by bacteria (Chrost and Overbek 1987) and algae (Boni et al. 1989). Since the algae themselves
release the enzymes and bacterial abundance is related to algal biomass, the rate of organic phosphorus
mineralization is related to algal biomass in model formulation. Another mechanism included in model
formulation is that algae stimulate production of an enzyme that mineralizes organic phosphorus to
phosphate when phosphate is scarce (Chrost and Overbek 1987; Boni et al. 1989). The formulations for
hydrolysis and mineralization rates including these processes are:
TJ~ rvr fLH± „
(4-44)
x-c,d,g
\-T ~ TRHDR^ (4-45)
4-20
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
KDOP = (KDP + d,g
KRP = minimum hydrolysis rate of refractory participate organic phosphorus (day"1)
KLP = minimum hydrolysis rate of labile particulate organic phosphorus (day"1)
KDP = minimum mineralization rate of dissolved organic phosphorus (day"1)
KRpaig, KLPalg = constants that relate hydrolysis of refractory and labile particulate organic phosphorus,
respectively, to algal biomass (day"1 per g C m"3)
KDpaig = constant that relates mineralization to algal biomass (day"1 per g C m"3)
KHP = mean half-saturation constant for algal phosphorus uptake (g P m"3).
= E
(4-47)
J x=c,d,g
When phosphate is abundant relative to KHP, the rates become close to the minimum values with little
influence from algal biomass. When phosphate becomes scarce relative to KHP, the rates increase with
the magnitude of increase depending on algal biomass. Equations 4-44 to 4-46 have exponential
functions that relate rates to temperature.
4.6 Nitrogen
The present model has five state variables for nitrogen: three organic forms (refractory particulate, labile
particulate, and dissolved) and two inorganic forms (ammonium and nitrate). The nitrate state variable in
the model represents the sum of nitrate and nitrite.
4.6.1 Particulate Organic Nitrogen
For refractory and labile particulate organic nitrogen, sources and sinks included in the model are (Figure
4-1):
• algal basal metabolism and predation
• dissolution to dissolved organic nitrogen
• settling
• external loads.
The kinetic equations for refractory and labile particulate organic nitrogen are:
= E (FNR-BMx + FNRP-PRx}ANC-Bx - K^-KPON
x=c,d,g,m
$ /7T7C DD^AA WRPON
+ ^(WSRp'RPON) + - Tr - (4-48)
oz V
4-21
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
- E (FNL-BMx + FNLP-PRx)ANC-Bx - KLPON-LPON
C x=c,d,g,m
$ /7T7C T D^AA WLPON
+ — (WSLp-LPON) + - - - (4-49)
oz V
RPON = concentration of refractory particulate organic nitrogen (g N m"3)
LPON = concentration of labile particulate organic nitrogen (g N m"3)
FNR,, = fraction metabolized nitrogen by algal group x as refractory particulate organic nitrogen
FNLX = fraction of metabolized nitrogen by algal group x produced as labile particulate organic
nitrogen
FNRP = fraction of predated nitrogen produced as refractory particulate organic nitrogen
FNLP = fraction of predated nitrogen produced as labile particulate organic nitrogen
ANCX = nitrogen-to-carbon ratio in algal group x (g N per g C)
= hydrolysis rate of refractory particulate organic nitrogen (day"1)
= hydrolysis rate of labile particulate organic nitrogen (day"1)
WRPON = external loads of refractory particulate organic nitrogen (g N day"1)
WLPON = external loads of labile particulate organic nitrogen (g N day"1).
4.6.2 Dissolved Organic Nitrogen
Sources and sinks for dissolved organic nitrogen included in the model are (Fig. 4-1):
• algal basal metabolism and predation
• dissolution from refractory and labile particulate organic nitrogen
• mineralization to ammonium
• external loads.
The kinetic equation describing these processes is:
(FNDxBMx + FNDP-PRx)ANCx-Bx
Oi x=c,d,g,m
+ KSKafRPON + KLPON-LPON - KDON-DON + (4.50)
DON = concentration of dissolved organic nitrogen (g N m"3)
FNDX = fraction of metabolized nitrogen by algal group x produced as dissolved organic nitrogen
FNDP = fraction of predated nitrogen produced as dissolved organic nitrogen
KDON = mineralization rate of dissolved organic nitrogen (day"1)
WDON = external loads of dissolved organic nitrogen (g N day"1).
4-22
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.6.3 Ammonium Nitrogen
Sources and sinks for ammonia nitrogen included in the model are (Fig. 4-1):
• algal basal metabolism, predation, and uptake
• mineralization from dissolved organic nitrogen
• nitrification to nitrate
• sediment-water exchange for the bottom layer only
• external loads.
The kinetic equation describing these processes is:
(FNI-BMx + FNIP-PRX - PN-PJANC^ + K^-DON
Ol x=c,d,g,m
xrv Art™ BFNH4 WNH4
- NU-NH4 + - + - (4-51)
FNI,, = fraction of metabolized nitrogen by algal group x produced as inorganic nitrogen
FNIP = fraction of predated nitrogen produced as inorganic nitrogen
PNX = preference for ammonium uptake by algal group x (0 < PNX < 1)
Nit = nitrification rate (day"1) given in Eq. 4-59
BFNH4 = sediment-water exchange flux of ammonium (g N m"2 day"1), applied to the bottom layer only
WNH4 = external loads of ammonium (g N day"1).
4.6.4 Nitrate Nitrogen
Sources and sinks for nitrate nitrogen included in the model are (Fig. 4-1):
• algal uptake
• nitrification from ammonium
• denitrification to nitrogen gas
• sediment-water exchange for the bottom layer only
• external loads.
The kinetic equation describing these processes is:
x + Nit-NH4 - ANDC-Denit-DOC
O x=c,d,g,m
BFNO3 WNO3
+ -AT + — <4-52)
ANDC = mass of nitrate nitrogen reduced per mass of dissolved organic carbon oxidized (0.933 g N
per g C from Eq. 4-33)
BFNO3 = sediment-water exchange flux of nitrate (g N m"2 day"1), applied to the bottom layer only
4-23
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
WNO3 = external loads of nitrate (g N day"1).
The remainder of this section explains each term in Equations 4-48 to 4-52, except BFNH4 and BFNO3
which are described in Chapter 5.
4.6.5 Effect of Algae on Nitrogen
The terms within summation (£) in Equations 4-48 to 4-52 account for the effects of algae on nitrogen.
As in phosphorus, both basal metabolism (respiration and excretion) and predation are considered, and
thus formulated, to contribute to organic and ammonium nitrogen. That is, algal nitrogen released by
both basal metabolism and predation are represented by distribution coefficients; FNR,,, FNLX, FNDX,
FNI,,, FNRP, FNLP, FNDP, and FNIP. The sum of four distribution coefficients for basal metabolism
should be unity; the sum of the predation distribution coefficients should also be unity.
Algae take up ammonium and nitrate for growth, and ammonium is preferred from thermodynamic
considerations. The preference of algae for ammonium is expressed as:
N03 KH^r
PN = NH4 - — - + NH4 - - - (4.53)
* (KHNx+NH4)(KHNx+NO3) (NH4 +NO3)(KHNx+NO3) v '
This equation forces the preference for ammonium to be unity when nitrate is absent, and to be zero
when ammonium is absent.
4.6.6 Mineralization and Hydrolysis
The third term on the RHS of Equations 4-48 and 4-49 represents hydrolysis of particulate organic
nitrogen and the last term in Eq. 4-50 represents mineralization of dissolved organic nitrogen. Including
a mechanism for accelerated hydrolysis and mineralization during nutrient-limited conditions (Section
4.5.7), the formulations for these processes are:
- TRHDR}) (4-54)
KDON = (KDN + KRN + NH4 +NQ3 KDNal^d S*)'eXPWwJr - TRMNlD (4-56)
Km = minimum hydrolysis rate of refractory particulate organic nitrogen (day"1)
KLN = minimum hydrolysis rate of labile particulate organic nitrogen (day"1)
4-24
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
KDN = minimum mineralization rate of dissolved organic nitrogen (day"1)
KiaMg, KLNalg = constants that relate hydrolysis of refractory and labile particulate organic nitrogen,
respectively, to algal biomass (day"1 per g C m~3)
KDNaig = constant that relates mineralization to algal biomass (day"1 per g C m"3)
KHN = mean half-saturation constant for algal nitrogen uptake (g N m"3).
= £
(4-57)
Equations 4-54 to 4-56 have exponential functions that relate rates to temperature.
4.6.7 Nitrification
Nitrification is a process mediated by autotrophic nitrifying bacteria that obtain energy through the
oxidation of ammonium to nitrite and of nitrite to nitrate. The stoichiometry of complete reaction is
(Bowie etal. 1985):
NH4+ + 2O2 - NO3~ + H2O + 2H+ (4-58)
The first term in the second line of Eq. 4-5 1 and its corresponding term in Eq. 4-52 represent the effect of
nitrification on ammonium and nitrate, respectively. The kinetics of complete nitrification process are
formulated as a function of available ammonium, dissolved oxygen and temperature:
AT-V DO NH4 AT-V f m
Nit = -- Nit-fvJT) C4-5Q^
KHNitDO + DO KHNitN + NH4 m m l4 ^'
fm(T) = exp(-KNitl[T - TNit]2) if T < TNit
= exp(-KNit2[TNit - 7]2) // T > TNit (4-60)
KHNitDO = nitrification half-saturation constant for dissolved oxygen (g O2 m"3)
KHNitN = nitrification half-saturation constant for ammonium (g N m"3)
Nitm = maximum nitrification rate at TNit (g N m"3 day"1)
TNit = optimum temperature for nitrification (°C)
KNitl = effect of temperature below TNit on nitrification rate (°C~2)
KNit2 = effect of temperature above TNit on nitrification rate (°C~2).
The Monod function of dissolved oxygen in Eq. 4-59 indicates the inhibition of nitrification at low
oxygen level. The Monod function of ammonium indicates that when ammonium is abundant, the
nitrification rate is limited by the availability of nitrifying bacteria. The effect of suboptimal temperature
is represented using Gaussian form.
4-25
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.6.8 Denitrification
The effect of denitrification on dissolved organic carbon was described in Section 4.4.5. Denitrification
removes nitrate from the system in stoichiometric proportion to carbon removal as determined by Eq.
4-33. The last term in the first line of Eq. 4-52 represents this removal of nitrate.
4.7 Silica
The present model has two state variables for silica: particulate biogenic silica and available silica.
4.7.1 Particulate Biogenic Silica
Sources and sinks for particulate biogenic silica included in the model are (Fig. 4-1):
• diatom basal metabolism and predation
• dissolution to available silica
• settling
• external loads
The kinetic equation describing these processes is:
^ = (FSPdBMd H- FSPP-PRd)ASCdBd - KSUA-SU + j-(WSd-SU) + -^ (4-61)
SU = concentration of particulate biogenic silica (g Si m"3)
FSPd = fraction of metabolized silica by diatoms produced as particulate biogenic silica
FSPP = fraction of predated diatom silica produced as particulate biogenic silica
ASCd = silica-to-carbon ratio of diatoms (g Si per g C)
KSUA = dissolution rate of particulate biogenic silica (day"1)
WSU = external loads of particulate biogenic silica (g Si day"1).
4.7.2 Available Silica
Sources and sinks for available silica included in the model are (Fig. 4-1):
• diatom basal metabolism, predation, and uptake
• settling of sorbed (particulate) available silica
• dissolution from particulate biogenic silica
• sediment-water exchange of dissolved silica for the bottom layer only
• external loads.
The kinetic equation describing these processes is:
= (FSIdBMd + FSIP-PR, - Pd}ASCdBd + KSUA-SU + -jj-(WSTSS-SAp)
4-26
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
BFSAd WSA
<4-62)
SA = concentration of available silica (g Si m"3) = SAd + SAp (4-63)
SAd = dissolved available silica (g Si m"3)
SAp = participate (sorbed) available silica (g Si m"3)
FSId = fraction of metabolized silica by diatoms produced as available silica
FSIP = fraction of predated diatom silica produced as available silica
BFSAd = sediment-water exchange flux of available silica (g Si m"2 day"1), applied to bottom layer only
WSA = external loads of available silica (g Si day"1).
In Eq. 4-62, if total active metal is chosen as a measure of sorption site, the settling velocity of total
suspended solid, WSTSS, is replaced by that of particulate metal, WSS (Sections 4.7.3 and 4.10).
4.7.3 Available Silica System
Analysis of Chesapeake Bay monitoring data indicates that silica shows similar behavior as phosphate in
the adsorption-desorption process (Cerco and Cole 1993). As in phosphate, therefore, available silica is
defined to include both dissolved and sorbed fractions (Eq. 4-63). Treatment of available silica is the
same as total phosphate, and the same method to partition dissolved and sorbed phosphate is used to
partition dissolved and sorbed available silica:
K^-TSS K^-TAMp
SAp = - ^ - SA or SAp = - ^ - — - SA (4-64)
1 + KSApTSS 1 + KSApTAMp ( >
SAd = - - - SA or SAd = - - - SA
1 H- K^-TSS 1 H- K^-TAMp
= SA - SAp (4-65)
KSAp = empirical coefficient relating available silica sorption to total suspended solid (per g m"3) or
particulate total active metal (per mol m"3) concentration.
As in KP04p in Section 4.5.4, KSAp is the ratio of sorbed to dissolved available silica per unit sorption site
available.
4.7.4 Effect of Diatoms on Silica
In Equations 4-62 and 4-63, those terms expressed as a function of diatom biomass (Bd) account for the
effects of diatoms on silica. As in phosphorus and nitrogen, both basal metabolism (respiration and
excretion) and predation are considered, and thus formulated, to contribute to particulate biogenic and
4-27
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
available silica. That is, diatom silica released by both basal metabolism and predation are represented
by distribution coefficients; FSPd, FSId, FSPP, and FSIP. The sum of two distribution coefficients for
basal metabolism should be unity and so is that for predation. Diatoms require silica as well as
phosphorus and nitrogen, and diatom uptake of available silica is represented by
(-Pd-ASCd-Bd)mEq.4-63.
4.7.5 Dissolution
The term (- KSUA-SU) in Eq. 4-62 and its corresponding term in Eq. 4-63 represent dissolution of
particulate biogenic silica to available silica. The dissolution rate is expressed as an exponential function
of temperature:
KSUA - KSU-^(KTSUA[T - TRSUA}} (4-66)
Ksu = dissolution rate of particulate biogenic silica at TRSUA (day"1)
KTSUA = effect of temperature on dissolution of particulate biogenic silica ("C"1)
TRSUA = reference temperature for dissolution of particulate biogenic silica (°C).
4.8 Chemical Oxygen Demand
In the present model, chemical oxygen demand is the concentration of reduced substances that are
oxidizable through inorganic means. The source of chemical oxygen demand in saline water is sulfide
released from sediments. A cycle occurs in which sulfate is reduced to sulfide in the sediments and
reoxidized to sulfate in the water column. In fresh water, methane is released to the water column by the
sediment process model. Both sulfide and methane are quantified in units of oxygen demand and are
treated with the same kinetic formulation. The kinetic equation, including external loads, if any, is:
DO T^~~ ^~ BFCOD WCOD
KCOD.COD
dt KHCOD + DO Az V
(4-67)
COD = concentration of chemical oxygen demand (g O2-equivalents m"3)
KHCOD = half-saturation constant of dissolved oxygen required for oxidation of chemical oxygen
demand (g O2 m"3)
KCOD = oxidation rate of chemical oxygen demand (day"1)
BFCOD = sediment flux of chemical oxygen demand (g O2-equivalents m"2 day"1), applied to bottom
layer only
WCOD = external loads of chemical oxygen demand (g O2-equivalents day"1).
An exponential function is used to describe the temperature effect on the oxidation rate of chemical
oxygen demand:
4-28
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
KCOD = KCD-Qxp(KTCOD[T - TRCOD\) (4-68)
KCD = oxidation rate of chemical oxygen demand at TR^n (day"1)
KTCOD = effect of temperature on oxidation of chemical oxygen demand (°C-1)
TRcOD = reference temperature for oxidation of chemical oxygen demand (°C).
4.9 Dissolved Oxygen
Sources and sinks of dissolved oxygen in the water column included in the model are (Fig. 4-1):
• algal photosynthesis and respiration
• nitrification
• heterotrophic respiration of dissolved organic carbon
• oxidation of chemical oxygen demand
• surface reaeration for the surface layer only
• sediment oxygen demand for the bottom layer only
• external loads.
The kinetic equation describing these processes is:
(1-3 - 0.3-/WJP - (1 - FCD) - — - BM\AOCR-BX
x x xJ x
dt x-g,m x x xKHRx + DO x
- AONT-NU-NH4 - AOCR-Km-DOC - - — - KCOD-COD
HR KHCQD + DO
r- /™ 7^ SOD WDO
+ Kr(D°S ~ D°) + -£- + —^~ (4-69
AONT = mass of dissolved oxygen consumed per unit mass of ammonium nitrogen nitrified (4.33 g O2
per g N; see Section 4.9.2)
AOCR = dissolved oxygen-to-carbon ratio in respiration (2.67 g O2 per g C; see Section 4.9.1)
Kr = reaeration coefficient (day"1): the reaeration term is applied to the surface layer only
DOS = saturated concentration of dissolved oxygen (g O2 m"3)
SOD = sediment oxygen demand (g O2 m"2 day"1), applied to the bottom layer only; positive is to the
water column
WDO = external loads of dissolved oxygen (g O2 day"1).
The two sink terms in Eq. 4-69, heterotrophic respiration and chemical oxygen demand, are explained in
Section 4.4.4 (Eq. 4-29) and Section 4.8 (Eq. 4-67), respectively. The remainder of this section explains
the effects of algae, nitrification, and surface reaeration.
4-29
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
4.9.1 Effect of Algae on Dissolved Oxygen
The first line on the RHS of Eq. 4-69 accounts for the effects of algae on dissolved oxygen. Algae
produce oxygen through photosynthesis and consume oxygen through respiration. The quantity produced
depends on the form of nitrogen utilized for growth. Equations describing production of dissolved
oxygen are (Morel 1983):
1Q6CO2 + 16NH4+ + H2PO4 + 1Q6H2O - protoplasm + 1Q6O2 + 15H+ (4-70)
106C02 + \6NO3 + H2PO4 + 122H2O + 17H+ - protoplasm + 138O2 (4-71)
When ammonium is the nitrogen source, one mole of oxygen is produced per mole of carbon dioxide
fixed. When nitrate is the nitrogen source, 1.3 moles of oxygen are produced per mole of carbon dioxide
fixed. The quantity, (1.3 - 0.3-PNX), in the first term of Eq. 4-69 is the photosynthesis ratio and
represents the molar quantity of oxygen produced per mole of carbon dioxide fixed. It approaches unity
as the algal preference for ammonium approaches unity.
The last term in the first line of Eq. 4-69 accounts for the oxygen consumption due to algal respiration
(Eq. 4-26). A simple representation of respiration process is:
CH2O + O2 = CO2 + H2O (4-72)
from which, AOCR = 2.67 g O2 per g C.
4.9.2 Effect of Nitrification on Dissolved Oxygen
The stoichiometry of nitrification reaction (Eq. 4-58) indicates that two moles of oxygen are required to
nitrify one mole of ammonium into nitrate. However, cell synthesis by nitrifying bacteria is
accomplished by the fixation of carbon dioxide so that less than two moles of oxygen are consumed per
mole ammonium utilized (Wezernak and Gannon 1968), i.e., AONT = 4.33 g O2 per g N.
4.9.3 Effect of Surface Reaeration on Dissolved Oxygen
The reaeration rate of dissolved oxygen at the air-water interface is proportional to the oxygen gradient
across the interface, (DOS - DO), when assuming the air is saturated with oxygen. The saturated
concentration of dissolved oxygen, which decreases as temperature and salinity increase, is specified
using an empirical formula (Genet et al. 1974):
DOs = 14.5532 - 0.38217-T + 5.4258xlO~3-T2
- CZ-(1.665xl(r4 - 5.866 xl(T6-r + 9.796 x l(T8-r2) (4-73)
CL = chloride concentration (mg/L) = S/l.80655.
4-30
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
The reaeration coefficient includes the effect of turbulence generated by bottom friction (OConnor and
Dobbins 1958) and that by surface wind stress (Banks and Herrera 1977):
K =
Kro
ueq
heq
B,
Wrei
Uw
KT
+ w
J_
Az
•KT
T - 20
= proportionality constant = 3.933 in MKS unit
= weighted velocity over cross-section (m sec"1) = E(ukVk)/E(Vk)
= weighted depth over cross-section (m) = ^(V^/B^
= width at the free surface (m)
= wind-induced reaeration (m day"1)
= 0.728 Uw* - 0.317 Uw + 0.0372 Uw2
= wind speed (m sec"1) at the height of 10 m above surface
= constant for temperature adjustment of DO reaeration rate.
(4-74)
(4-75)
4.10 Total Active Metal
The present model requires simulation of total active metal for adsorption of phosphate and silica if that
option is chosen (Fig. 4-1). The total active metal state variable is the sum of iron and manganese
concentrations, both particulate and dissolved. In the model, the origin of total active metal is benthic
sediments. Since sediment release of metal is not explicit in the sediment model (see Chapter 5), release
is specified in the kinetic portion of the water column model. The only other term included is settling of
the particulate fraction. Then the kinetic equation for total active metal, including external loads, if any,
may be written as:
dTAM = KHbmf BFTAM Ktan(T. Ttam} + d_
dt
KHbmf + DO Az
dz
-TAMp)
V
(4-76)
(4-77)
TAM = total active metal concentration (mol m"3) = TAMd + TAMp
TAMd = dissolved total active metal (mol m"3)
TAMp = particulate total active metal (mol m"3)
KHbmf = dissolved oxygen concentration at which total active metal release is half the anoxic release
rate (g O2 m"3)
BFTAM = anoxic release rate of total active metal (mol m"2 day"1), applied to the bottom layer only
Ktam = effect of temperature on sediment release of total active metal (^'C'1)
Ttam = reference temperature for sediment release of total active metal (°C)
WSS = settling velocity of particulate metal (m day"1)
WTAM = external loads of total active metal (mol day"1).
4-31
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
In estuaries, iron and manganese exist in particular and dissolved forms depending on dissolved oxygen
concentration. In the oxygenated water, most of the iron and manganese exist as particulate while under
anoxic conditions, large fractions are dissolved, although solid-phase sulfides and carbonates exist and
may predominate. The partitioning between particulate and dissolved phases is expressed using a
concept that total active metal concentration must achieve a minimum level, which is a function of
dissolved oxygen, before precipitation occurs:
TAMd = mmimum{TAMdmx-exp(- Kdotam-DO) , TAM} (4-78)
TAMp = TAM - TAMd (4-79)
TAMdmx = solubility of total active metal under anoxic conditions (mol m"3)
Kdotam = constant that relates total active metal solubility to dissolved oxygen (per g O2 m"3).
4.11 Fecal Coliform Bacteria
Fecal coliform bacteria are indicative of organisms from the intestinal tract of humans and other animals
and can be used as an indicator bacteria as a measure of public health (Thomann and Mueller 1987). In
the present model, fecal coliform bacteria have no interaction with other state variables, and have only
one sink term, die-off The kinetic equation, including external loads, may be written as:
9FCB = - KFCB-TFCBT ~ 20-FCB + WFCB
dt v
FCB = bacteria concentration (MPN per 100 ml)
KFCB = first order die-off rate at 20°C (day1)
TFCB = effect of temperature on decay of bacteria (°C-1)
WFCB = external loads of fecal coliform bacteria (MPN per 100 ml m3 day"1).
4.12 Method of Solution
The kinetic equations for the 21 state variables in the EFDC water column water quality model can be
expressed in a 21 x 21 matrix after linearizing some terms, mostly Monod type expressions:
(4-80)
-[q = [K\-[C] + [R\ (4-81)
ot
where [C] is in mass volume"1, [K] is in time"1, and [R] is in mass volume"1 time"1. Since the settling of
particulate matter from the overlying cell acts as an input for a given cell, when Eq. 4-8 1 is applied to a
cell of finite volume, it may be expressed as:
l + [R]k (4-82)
4-32
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
where the four matrices [C], [Kl], [K2], and [R] are defined in Appendix A of Park et al. (1995). The
subscript k designates a cell at the k* vertical layer. The layer index k increases upward with KC vertical
layers; k = 1 is the bottom layer and k = KC is the surface layer. Then A = 0 for k = KC; otherwise,
A = 1 . The matrix [K2] is a diagonal matrix, and the non-zero elements account for the settling of
particulate matter from the overlying cell.
Equation 4-82 is solved using a second-order accurate trapezoidal scheme over a time step of 0, which
may be expressed as:
(4-83)
where 0 = 2-m- At is the time step for the kinetic equations; [I] is a unit matrix; [C]A = [C]N + [C]°; the
superscripts O and N designate the variables before and after being adjusted for the relevant kinetic
processes. Since Eq. 4-83 is solved from the surface layer downward, the term with [C]k+1A is known for
the k* layer and thus placed on the RHS. In Eq. 4-83, inversion of a matrix can be avoided if the 21 state
variables are solved in a proper order. The kinetic equations are solved in the order of the variables in
the matrix [C] defined in Appendix A of Park et al. (1995).
4.13 Macroalgae (Periphyton) State Variable
The EFDC water quality model was augmented to represent benthic attached algae (often referred to as
macroalgae in estuarine waters and periphyton in fresh waters) using the existing framework for
phytoplankton growth kinetics. Mathematical relationships based on the impacts of temperature,
available light, available nutrients, stream velocity, and density-dependent interactions were incorporated
into the algae growth kinetics framework within EFDC. The major differences between modeling
techniques for attached and free-floating algae are: (1) attached algae are expressed in terms of areal
densities rather than volumetric concentrations; (2) attached algae growth can be limited by the
availability of bottom substrate; (3) the availability of nutrients to the macroalgae matrix can be
influenced by stream velocity; and (4) macroalgae are not subject to hydrodynamic transport. A good
description of periphyton kinetics as it relates to water quality modeling can be found in Warwick et al.
(1997) and has been used to develop this section of the report.
A mass-balance approach was used to model macroalgae growth, with carbon serving as the measure of
standing crop size or biomass. For each model grid cell the equation for macroalgae growth is slightly
different than the one for free-floating algae (Eq. 4-6):
Pm = PMm -f^N) -f2(I) -/3(7) /4(F) '/5(Z>) (4-84)
4-33
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
where
PMm = maximum growth rate under optimal conditions for macroalgae
f[(N) = effect of suboptimal nutrient concentration (0 < f, < 1)
f2(I) = effect of suboptimal light intensity (0 < f2 < 1)
f3(T) = effect of suboptimal temperature (0 < f3 < 1)
f4(V) = velocity limitation factor (0 < f4 < 1)
f5(D) = density-dependent growth rate reduction factor (0 < f5 < 1).
The basic growth kinetics for macroalgae were developed from those supplied by EFDC and others
developed by Runke (1985). The macroalgae population as a whole is characterized by the total biomass
present without considering the different species and their associated environmental processes. The
optimum growth for the given temperature is adjusted for light, nutrients, velocity, and density-dependent
limitations. Each growth limitation factor can vary from 0 to 1. A value of 1 indicates the factor does
not limit growth, and a value of 0 means the factor is so severely limiting that growth is stopped entirely
(Bowie etal. 1985).
Stream velocity has a twofold effect on periphyton productivity in freshwater streams: velocity increases
to a certain level to enhance biomass accrual, but further increases result in substantial scouring (Horner
et al. 1990). A benthic algal population is typified as a plant community with an understory and an
overstory. The entire community is called a matrix. As the matrix develops, the periphyton community
is composed of an outer layer of photosynthetically active cells and inner layers of senescent and
decomposing cells. Layering can also develop among different species of periphyton. Environmental
conditions within the matrix are altered by the physical structure of the periphyton. This influences
nutrient uptake and primary production rates of the algae (Sand-Jensen 1983). Above a certain level,
current has a simulating effect on periphyton metabolism by mixing the overlying waters with nutrient-
poor waters that develop around cells (Whitford and Schumacher 1964). The physical structure of the
periphyton community and nutrient uptake by periphyton interfere with nutrient flux through the
microbial matrix (Stevenson and Glover 1993).
Current is constantly scouring periphyton from its substrate. At high enough velocities, shear stress can
result in substantial biomass reduction. Even at low velocities, sudden increases in velocity raise
instantaneous loss rates substantially, but these high rates persist only briefly (Horner et al. 1990). An
increase in velocity above that to which benthic algae are accustomed leads to increased loss rates and
temporarily reduced biomass. However, recolonization and growth after biomass reduction are usually
rapid. The effects of suboptimal velocity upon growth rate are represented in the model by a velocity
limitation function. Two options are available in the model for specifying the velocity limitation: (1) a
4-34
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Michaelis-Menton (or Monod) equation (4-85) and (2) a five-parameter logistic function (4-86). The
Monod equation limits macroalgae growth due to low velocities, whereas the five-parameter logistic
function can be configured to limit growth due to either low or high velocities (Figure 4-2).
Velocity Limitation Function
for Control of Periphyton Growth
o
CO
O
W-»
CO
O
0.25
0.50
0.75 1.00
Velocity (ft/sec)
1.25
1.50
1.75
2.00
Logistic function
Monod Function
Figure 4-2. Velocity limitation function for (Option 1) the Monod equation where KMV = 0.25
m/sec and KMVmin=0.15 m/sec, and (Option 2) the five-parameter logistic function
where a=1.0, b=12.0, c=0.3, c/=0.35, and e=3.0 (high velocities are limiting)
Velocity limitation option 1, the Michaelis-Menton equation, is written as follows:
_ U
~ KMV + U
(4-85)
where
U = stream velocity (m/sec)
KMV = half-saturation velocity (m/sec)
Velocity limitation option 2, the five-parameter logistic function is as follows:
a - d
- d
(->* r
c
(4-86)
where
U = stream velocity (m/sec)
a = asymptote at minimum x
4-35
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
b = slope after asymptote a
c = x-translation
d = asymptote at maximum x
e = slope before asymptote d
The half-saturation velocity in Eq. 4-85 is the velocity at which half the maximum growth rate occurs.
This effect is analogous to the nutrient limitation because the effect of velocity at suboptimal levels on
periphyton growth is due to increasing the exchange of nutrients between the algal matrix and the
overlying water (Runke 1985). However, this formula can be too limiting at low velocities. This
function does not allow periphyton growth in still waters, but periphyton does grow in still waters such as
lakes. Therefore, the function is applied only at velocities above a minimum threshold level (KMVmin).
When velocities are at or below this lower level, the limitation function is applied at the minimum level.
Above this velocity, the current produces a steeper diffusion gradient around the periphyton (Whitford
and Schumacher 1964). A minimum formulation is used to combine the limiting factors for nitrogen,
phosphorus, and velocity. The most severely limiting factor alone limits periphyton growth. Note that
Eq. 4-86 can be configured so that low velocities are limiting by setting parameter d greater than
parameter a, and vice versa to limit growth due to high velocities. In waters that are rich in nutrients, low
velocities will not limit growth. However, high velocities may cause scouring and detachment of the
macroalgae, resulting in a reduction in biomass. The five -parameter logistic function can be configured
to approximate this reduction by limiting growth at high velocities.
Macroalgae (periphyton) growth can also be limited by the availability of suitable substrate (Ross 1983).
Macroalgae communities reach maximum rates of primary productivity at low levels of biomass
(Mclntire 1973; Pfeifer and McDiffett 1975). The relationship between standing crop and production
employs the Michaelis-Menton kinetic equation:
KBP+P (4'87)
tn
where
KBP = half-saturation biomass level (g C/m2)
Pm = macroalgae biomass level (g C/m2).
The half-saturation biomass level (KBP) is the biomass at which half the maximum growth rate occurs.
Caupp et al. (1991) used a KBP value of 5.0 g C/m2 (assuming 50% of ash free dry mass is carbon) for a
region of the Truckee River system in California. The function in Eq. 4-87 allows maximum rates of
primary productivity at low levels of biomass with decreasing rates of primary productivity as the
community matrix expands.
4-36
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
5-EFDC SEDIMENT PROCESS MODEL
A sediment process model developed by DiToro and Fitzpatrick (1993; hereinafter referred to as D&F)
and was coupled with CE-QUAL-ICM for Chesapeake Bay water quality modeling (Cerco and Cole
1993). The sediment process model was slightly modified and incorporated into the EFDC water quality
model to simulate the processes in the sediment and at the sediment-water interface. The description of
the EFDC sediment process model in this section is from Park et al. (1995). The sediment process model
has 27 water-quality related state variables and fluxes (Table 5-1).
Table 5-1. EFDC sediment process model state variables and flux terms
(1) particulate organic carbon Gl class in layer 2 (15) nitrate nitrogen in layer 1
(2) parti culate organic carbon G2 class in layer 2 (16) nitrate nitrogen in layer 2
(3) participate organic carbon G3 class in layer 2 (17) phosphate phosphorus in layer 1
(4) particulate organic nitrogen Gl class in layer 2 (18) phosphate phosphorus in layer 2
(5) particulate organic nitrogen G2 class in layer 2 (19) available silica in layer 1
(6) particulate organic nitrogen G3 class in layer 2 (20) available silica in layer 2
(7) particulate organic phosphorus Gl class in layer 2 (21) ammonia nitrogen flux
(8) particulate organic phosphorus G2 class in layer 2 (22) nitrate nitrogen flux
(9) particulate organic phosphorus G3 class in layer 2 (23) phosphate phosphorus flux
(10) particulate biogenic silica in layer 2 (24) silica flux
(11) sulfide/methane in layer 1 (25) sediment oxygen demand
(12) sulfide/methane in layer 2 (26) release of chemical oxygen demand
(13) ammonia nitrogen in layer 1 (27) sediment temperature
(14) ammonia nitrogen in layer 2
The nitrate state variables, (15), (16), and (22), in the model represent the sum of nitrate and nitrite
nitrogen. The three G classes for particulate organic matter (POM) in Layer 2 and the two layers for
inorganic substances are described below.
In the sediment model, benthic sediments are represented as two layers (Fig. 5-1). The upper layer
(Layer 1) is in contact with the water column and may be oxic or anoxic depending on dissolved oxygen
concentration in the overlying water. The lower layer (Layer 2) is permanently anoxic. The upper layer
depth, which is determined by the penetration of oxygen into the sediments, is at its maximum only a
small fraction of the total depth. Because Hj (~ 0.1 cm) « H2,
H = Hj + H2 * H2 (5-1)
5-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
WATER COLUMN
DE
01
LU
SEDIMEN
,2 L/
UL
LLJ
^
POSITION A L fd-v
1
SURFACE MASS TRANSFER (s): ., , fct]
PARTITIONING: fd ^ ^fn
\ \
REACTIONS: KI
SED
r ^
DIAC
(MENTATION (W):
PARTICLE DIFFUSION
MIXING A (KL)
(5)) t
1 I
3ENESIS: J ^ Ct
PARTITIONING: fdj* ^ fp~
REACTIONS: Kj
SEDIMENTATION (W):
1
Figure 5-1. Sediment layers and processes included in sediment process model.
where H is the total depth (approximately 10 cm), Hj is the upper layer depth and H2 is the lower layer
depth.
The model incorporates three basic processes (Fig. 5-2): (1) depositional flux of POM, (2) the diagenesis
of POM, and (3) the resulting sediment flux. The sediment model is driven by net settling of particulate
organic carbon, nitrogen, phosphorus, and silica from the overlying water to the sediments (depositional
flux). Because of the negligible thickness of the upper layer (Eq. 5-1), deposition is considered to
proceed from the water column directly to the lower layer. Within the lower layer, the model simulates
the diagenesis (mineralization or decay) of deposited POM, which produces oxygen demand and
inorganic nutrients (diagenesis flux). The third basic process is the flux of substances produced by
diagenesis (sediment flux). Oxygen demand, as sulfide (in salt water) or methane (in fresh water), takes
three paths out of the sediments: (1) oxidation at the sediment-water interface as sediment oxygen
demand, (2) export to the water column as chemical oxygen demand, or (3) burial to deep, inactive
5-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
sediments. Inorganic nutrients produced by diagenesis take two paths out of the sediments: (1) release to
the water column or (2) burial to deep, inactive sediments (Fig. 5-2).
NET SETTLING
OF POM
Figure 5-2. Schematic diagram for sediment process model
This section describes the three basic processes with reactions and sources/sinks for each state variable.
The method of solution includes finite difference equations, solution scheme, boundary, and initial
conditions. Complete model documentation can be found in D&F (1993).
5.1 Depositional Flux
Deposition is one process that couples the water column model with the sediment model. Consequently,
deposition is represented in both the water column and sediment models. In the water column model, the
governing mass-balance equations for the following state variables contain settling terms, which
represent the depositional fluxes:
• three algal groups, cyanobacteria, diatoms and green algae (Eq. 4-5)
• refractory and labile particulate organic carbon (Equations 4-20 and 4-21)
• refractory and labile particulate organic phosphorus (Equations 4-35 and 4-36) and total
phosphate (Eq. 4-38)
• refractory and labile particulate organic nitrogen (Equations 4-48 and 4-49)
• particulate biogenic silica (Eq. 4-61) and available silica (Eq. 4-62).
5-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
The sediment model receives these depositional fluxes of particulate organic carbon (POC), particulate
organic nitrogen (PON), particulate organic phosphorus (POP), and particulate biogenic silica (PSi).
Because of the negligible thickness of the upper layer (Eq. 5-1), deposition is considered to proceed from
the water column directly to the lower layer. Since the sediment model has three G classes of POM, G, (i
= 1, 2, or 3), depending on the time scales of reactivity (Section 5.2), the POM fluxes from the water
column should be mapped into three G classes based on their reactivity. Then the depositional fluxes for
the ith G class (i = 1, 2, or 3) may be expressed as:
Jpoc,i = FCLP-WSLP-LPOCN + FCRP-WSRP-RPOCN + £ FCBX,{WSX'B* (5-2)
x=c,dg
JpoN,i = FNLPiWSLp-LPONN + FNRPiWSRp-RPONN+ £ FNBxiANCxWSxBxN (5.3)
x=c,dg
Jpopt = FPLPi-WSLp-LPOPN + FPRPi-WSRp-RPOPN + £ FPBxiAPC-WSxBxN
(5-4)
JPSi = WSd-SUN + ASCd-WSd-Bd + WSTSS-SApN (5-5)
JPOM.I = depositional flux of POM (M = C, N or P) routed into the i*11 G class (g m"2 day"1)
JPSl = depositional flux of PSi (g Si m"2 day"1)
FCLP15 FNLP15 FPLP, = fraction of water column labile POC, PON, and POP, respectively, routed into
the i* G class in sediment
FCRP15 FNRP15 FPRP, = fraction of water column refractory POC, PON, and POP, respectively, routed
into the i*11 G class in sediment
FCBX 15 FNBX 15 FPBX , = fraction of POC, PON, and POP, respectively, in the algal group x routed into
the i* G class in sediment
Yi = 1 for i = 1
0 fori = 2or3.
In the source code, the sediment process model is solved after the water column water quality model, and
the calculated fluxes using the water column conditions at t = tn are used for the computation of the water
quality variables at t = tn+ 0. The superscript N indicates the variables after being updated for the kinetic
processes, as defined in Eq. 4-82.
5-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
The settling of sorbed phosphate is considered to contribute to the labile Gj pool in Eq. 5-4, and settling
of sorbed silica contributes to JPSl in Eq. 5-5 to avoid creation of additional depositional fluxes for
inorganic particulates. The sum of distribution coefficients should be unity: £; FCLP, = £; FNLP, = £;
FPLP, = £, FCRP, = I, FNRP, = I, FPRP, = JV FCBX, = JV FNBX, = JV FPBX, = 1. The settling
velocities, WSLP, WSRP, WSX, and WSTSS, as defined in the EFDC water column model (Section 4), are
net settling velocities. If total active metal is selected as a measure of sorption site, WSTSS is replaced by
WSS in Equations 5-4 and 5-5 (see Sections 4.5 and 4.7).
5.2 Diagenesis Flux
Another coupling point of the sediment model to the water column model is the sediment flux, which is
described in Section 5.3. The computation of sediment flux requires that the magnitude of the diagenesis
flux be known. The diagenesis flux is explicitly computed using mass-balance equations for deposited
POC, PON, and POP. (Dissolved silica is produced in the sediments as the result of the dissolution of
PSi. Since the dissolution process is different from the bacterial-mediated diagenesis process, it is
presented separately in Section 5.4.) In the mass-balance equations, the depositional fluxes of POM are
the source terms and the decay of POM in the sediments produces the diagenesis fluxes. The integration
of the mass-balance equations for POM provides the diagenesis fluxes that are the inputs for the mass-
balance equations for ammonium, nitrate, phosphate, and sulfide/methane in the sediments (Section 5.3).
The difference in decay rates of POM is accounted for by assigning a fraction of POM to various decay
classes (Westrisch and Berner 1984). POM in the sediments is divided into three G classes, or fractions,
representing three scales of reactivity. The Gl (labile) fraction has a half life of 20 days, and the G2
(refractory) fraction has a half life of one year. The G3 (inert) fraction is nonreactive, i.e., it undergoes
no significant decay before burial into deep, inactive sediments. The varying reactivity of the G classes
controls the time scale over which changes in depositional flux will be reflected in changes in diagenesis
flux. If the G! class would dominate the POM input into the sediments, then there would be no
significant time lag introduced by POM diagenesis and any changes in depositional flux would be readily
reflected in diagenesis flux.
Because the upper layer thickness is negligible (Eq. 5-1) and thus depositional flux is considered to
proceed directly to the lower layer (Equations 5-2 to 5-5), diagenesis is considered to occur in the lower
layer only. The mass-balance equations are similar for POC, PON, and POP, and for different G classes.
5-5
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
The mass-balance equation in the anoxic lower layer for the iih G class (i = 1, 2, or 3) may be expressed
as:
POMj _ _ -TS- (±T 20 ^i TT _ TTT7 f~1 , J
~ ~ JPOM,i
GPOM,i = concentration of POM (M = C, N, or P) in the i* G class in Layer 2 (g m"3)
KPOMl = decay rate of the 1th G class POM at 20°C in Layer 2 (day1)
QPOM.I = constant for temperature adjustment for KPOHi
T = sediment temperature (°C)
W = burial rate (m day"1).
Since the G3 class is inert, KPOM3 = 0.
Once the mass-balance equations for GPOM1 and GPOM2 are solved, the diagenesis fluxes are computed
from the rate of mineralization of the two reactive G classes:
2
7 = Y^ K flr"20 C ff C5 7^
M ^ POM,, ,1 POM,,
JM = diagenesis flux (g ~2 day"1) of carbon (M = C), nitrogen (M = N), or phosphorus (M = P).
5.3 Sediment Flux
The mineralization of POM produces soluble intermediates, which are quantified as diagenesis fluxes in
the previous section. The intermediates react in the oxic and anoxic layers, and portions are returned to
the overlying water as sediment fluxes. Computation of sediment fluxes requires mass-balance equations
for ammonium, nitrate, phosphate, sulfide/methane, and available silica. This section describes the flux
portion for ammonium, nitrate, phosphate, and sulfide/methane of the model. Available silica is
described in Section 5.4.
In the upper layer, the processes included in the flux portion are (Fig. 5-1)
• exchange of dissolved fraction between Layer 1 and the overlying water
• exchange of dissolved fraction between Layer 1 and 2 via diffusive transport
• exchange of particulate fraction between Layer 1 and 2 via particle mixing
• loss by burial to the lower layer (Layer 2)
• removal (sink) by reaction
• internal sources.
5-6
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Since the upper layer is quite thin, Hj ~ 0.1 cm (Eq. 5-1) and the surface mass transfer coefficient (s) is
on the order of 0.1 m day"1, then the residence time in the upper layer is H^s ~ 10"2 days. Hence, a
steady-state approximation is made in the upper layer. Then the mass-balance equation for ammonium,
nitrate, phosphate, or sulfide/methane in the upper layer is:
dCt,
^ = 0 = s(fd0-Ct0 -fd.-Ct,} + KL(fd2-Ct2 -fd.-Ct,}
at
K2
+ 6)(^2-a2 - fpv-Ct^ - W-Ctv - — Qj + Jv (5-8)
s
Ctj & Ct2 = total concentrations in Layer 1 and 2, respectively (g m"3)
Ct0 = total concentration in the overlying water (g m"3)
s = surface mass transfer coefficient (m day"1)
KL = diffusion velocity for dissolved fraction between Layer 1 and 2 (m day"1)
fr> = particle mixing velocity between Layer 1 and 2 (m day"1)
fd0 = dissolved fraction of total substance in the overlying water (0 < fd0 < 1)
fdj = dissolved fraction of total substance in Layer 1 (0 < fdj < 1)
fpj = particulate fraction of total substance in Layer 1 (= 1 - fdj)
fd2 = dissolved fraction of total substance in Layer 2 (0 < fd2 < 1)
fp2 = particulate fraction of total substance in Layer 2 (= 1 - fd2)
KJ = reaction velocity in Layer 1 (m day"1)
Jj = sum of all internal sources in Layer 1 (g m"2 day"1).
The first term on the RHS of Eq. 5-8 represents the exchange across sediment-water interface. Then the
sediment flux from Layer 1 to the overlying water, which couples the sediment model to the water
column model, may be expressed as:
Jaq-s(fd,-Ct, -fdo-Cto) (5-9)
Jaq = sediment flux of ammonium, nitrate, phosphate, or sulfide/methane to the overlying water
(g m"2 day"1).
The convention used in Eq. 5-9 is that positive flux is from the sediment to the overlying water.
In the lower layer, the processes included in the flux portion are (Fig. 5-1)
• exchange of dissolved fraction between Layer 1 and 2 via diffusive transport
• exchange of particulate fraction between Layer 1 and 2 via particle mixing
• deposition from Layer 1 and burial to the deep inactive sediments
5-7
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
• removal (sink) by reaction
• internal sources including diagenetic source.
The mass-balance equation for ammonium, nitrate, phosphate or sulfide/methane in the lower layer is:
d a,
H2— ± = - KL(fd2-Ct2 -fd.-Ct,) - K(fp2-Ct2 - fa-CtJ
at
+ w(cti - a2) - K2-a2 + J2 (5-io)
K2 = reaction velocity in Layer 2 (m day"1)
J2 = sum of all internal sources including diagenesis in Layer 2 (g m"2 day"1).
The substances produced by mineralization of POM in sediments may be present in both dissolved and
particulate phases. This distribution directly affects the magnitude of the substance that is returned to the
overlying water. In Equations 5-8 to 5-10, the distribution of a substance between the dissolved and
particulate phases in a sediment is parameterized using a linear partitioning coefficient. The dissolved
and particulate fractions are computed from the partitioning equations:
(5-iD
2 = TT fa = l-U (5-12)
m1; m2 = solid concentrations in Layer 1 and 2, respectively (kg L"1)
TCJ, K2 = partition coefficients in Layer 1 and 2, respectively (per kg L"1).
The partition coefficient is the ratio of particulate to dissolved fraction per unit solid concentration (i.e.,
per unit sorption site available).
All terms, except the last two terms, in Equations 5-8 and 5-10 are common to all state variables and are
described in Section 5.3.1. The last two terms represent the reaction and source/sink terms, respectively.
These terms, which take different mathematical formulations for different state variables, are described
in Sections 5.3.2 to 5.3.5 for ammonium, nitrate, phosphate, and sulfide/methane, respectively.
5.3.1 Common Parameters for Sediment Flux
Parameters that are needed for the sediment fluxes are s, So, KL, W, H2, mls m2, TCJ, K2, KI; K2, J1; and J2 in
Equations 5-8 to 5-12. Of these, KI; K2, J1; and J2 are variable-specific. Among the other common
5-8
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
parameters, W, H2, m1; and m2, are specified as input. The modeling of the remaining three parameters,
s, fr>, and KL, is described in this section.
5.3.1.1 Surface mass transfer coefficient. Owing to the observation that the surface mass transfer
coefficient, s, can be related to the sediment oxygen demand, SOD (DiToro et al. 1990), s can be
estimated from the ratio of SOD and overlying water oxygen concentration:
flj DO,
Dj = diffusion coefficient in Layer 1 (m2 day"1).
Knowing s, it is possible to estimate the other model parameters.
5.3.1.2 Particulate phase mixing coefficient. The particle mixing velocity between Layer 1 and 2 is
parameterized as:
'~2° D00
Dp = apparent diffusion coefficient for particle mixing (m2 day"1)
0Dp = constant for temperature adjustment for Dp
GPQC.R = reference concentration for GPOC,i (g C m"3)
KMDp = particle mixing half-saturation constant for oxygen (g O2 m"3).
The enhanced mixing of sediment particles by macrobenthos (bioturbation) is quantified by estimating
Dp. The particle mixing appears to be proportional to the benthic biomass (Matisoff 1982), which is
correlated to the carbon input to the sediment (Robbins et al. 1989). This is parameterized by assuming
that benthic biomass is proportional to the available labile carbon, GPOCjl, and GPOC,R is the reference
concentration at which the particle mixing velocity is at its nominal value. The Monod-type oxygen
dependency accounts for the oxygen dependency of benthic biomass.
It has been observed that a hysteresis exists in the relationship between the bottom water oxygen and
benthic biomass. Benthic biomass increases as the summer progresses. However, the occurrence of
anoxia/hypoxia reduces the biomass drastically and also imposes stress on benthic activities. After full
overturn, the bottom water oxygen increases, but the population does not recover immediately. Hence,
the particle mixing velocity, which is proportional to the benthic biomass, does not increase in response
5-9
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
to the increased bottom water oxygen. Recovery of benthic biomass following hypoxic events depends
on many factors including severity and longevity of hypoxia, constituent species, and salinity (Diaz and
Rosenberg 1995).
This phenomenon of reduced benthic activities and hysteresis is parameterized based on the idea of stress
that low oxygen imposes on the benthic population. It is analogous to the modeling of the toxic effect of
chemicals on organisms (Mancini 1983). A first order differential equation is employed, in which the
benthic stress (1) accumulates only when overlying oxygen is below KMDp and (2) is dissipated at a first
order rate (Fig. 5-3a):
dST
- K^-ST +
dt •"J"z)P
1 -
D0n
// DO, < KMDp
(5-15)
if DO, > KMDp
ui
ST = accumulated benthic stress (day)
KST = first order decay rate for ST (day"1).
The behavior of this formulation can be understood by evaluating the steady-state stresses at two extreme
conditions of overlying water oxygen, DO0:
as DO0 = 0 KST- ST = 1 f(ST) = (1 - KST- ST) = 0
as DO0 > KMDp KST- ST = 0 f(ST) = (1 - KST- ST) = 1
The dimensionless expression, f(ST) = 1 - KST-ST, appears to be the proper variable to quantify the effect
of benthic stress on benthic biomass and thus particle mixing (Fig. 5-3b).
The final formulation for the particle mixing velocity, including the benthic stress, is:
H —^ (5-16)
H2
Dpmm = minimum diffusion coefficient for particle mixing (m2 day"1).
The reduction in particle mixing due to the benthic stress, f(ST), is estimated by employing the following
procedure. The stress, ST, is normally calculated with Eq. 5-15. Once DO0 drops below a critical
concentration, DOST c, for NChypoxm consecutive days or more, the calculated stress is not allowed to
decrease until tMBS days of DO0 > DOST c. That is, only when hypoxic days are longer than critical
5-10
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
hypoxia days (NChypoxia), the maximum stress, or minimum (1 - KST-ST), is retained for a specified period
(tj^g days) after DO0 recovery (Fig. 5-3). No hysteresis occurs if DO0 does not drop below DOST c or if
hypoxia lasts less than NChypoxia days. When applying maximum stress for tMBS days, the subsequent
hypoxic days are not included in t^g. This parameterization of hysteresis essentially assumes seasonal
hypoxia, i.e., one or two major hypoxic events during summer, and might be unsuitable for systems with
multiple hypoxic events throughout a year.
25 -i
1.0
0.8 -
0.6 -
0.2 -
0.0
(b)
TIME
Figure 5-3. Benthic stress (a) and its effect on particle mixing (b) as a function of overlying
water column dissolved oxygen concentration
Three parameters relating to hysteresis, DOST c, NChypoxia, and t^g, are functions of many factors including
severity and longevity of hypoxia, constituent species, and salinity, and thus have site-specific
variabilities (Diaz and Rosenberg 1995). The critical overlying oxygen concentration, DOSTc, also
depends on the distance from the bottom of the location of DO0 The critical hypoxia days, NChypoxm,
depend on tolerance of benthic organisms to hypoxia and thus on benthic community structure (Diaz and
Rosenberg 1995). The time lag for the recovery of benthic biomass following hypoxic events, t^g, tends
to be longer for higher salinity. The above three parameters are considered to be spatially constant input
parameters.
5-11
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
5.3.1.3 Dissolved phase mixing coefficient. Dissolved phase mixing between Layer 1 and 2 is via
passive molecular diffusion, which is enhanced by the mixing activities of the benthic organisms (bio-
irrigation). This is modeled by increasing the diffusion coefficient relative to the molecular diffusion
coefficient:
D • QT ~ 20
KL = d ™ + R -to (5-17)
H2
Dd = diffusion coefficient in pore water (m2 day"1)
0Dd = constant for temperature adjustment for Dd
RBI.BT = ratio of bio-irrigation to bioturbation.
The last term in Eq. 5-17 accounts for the enhanced mixing by organism activities.
5.3.2 Ammonia Nitrogen
Diagenesis is assumed not to occur in the upper layer because of its shallow depth, and ammonium is
produced by diagenesis in the lower layer:
JN (ft"0™ El- 5~7) (5-18)
Ammonium is nitrified to nitrate in the presence of oxygen. A Monod-type expression is used for the
ammonium and oxygen dependency of the nitrification rate. Then the oxic layer reaction velocity in
Eq. 5-8 for ammonium may be expressed as:
2 .07-20
KNH4 VNH4 (5-19)
- „, , „„ „, , ,„,,
2'KMNH4,02 + D00 KMNH4 + NH4\
and then the nitrification flux becomes:
, (5-20)
5
KMNH4 02 = nitrification half-saturation constant for dissolved oxygen (g O2 m"3)
NH4j = total ammonium nitrogen concentration in Layer 1 (g N m"3)
KMNH4 = nitrification half-saturation constant for ammonium (g N m"3)
KNH4 = optimal reaction velocity for nitrification at 20°C (m day"1)
0NH4 = constant for temperature adjustment for KNH4
JNlt = nitrification flux (g N m"2 day"1).
Nitrification does not occur in the anoxic lower layer:
5-12
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
° (5-21)
Once Equations 5-8 and 5-10 are solved for NH4j and NH42, the sediment flux of ammonium to the
overlying water, JaqNH4, can be calculated using Eq. 5-9. Note that it is not NH4j and NH42 that
determine the magnitude of JaqNH4 (Section X-B-2 in D&F 1993). The magnitude is determined by
(1) the diagenesis flux, (2) the fraction that is nitrified, and (3) the surface mass transfer coefficient (s)
that mixes the remaining portion.
5.3.3 Nitrate Nitrogen
Nitrification flux is the only source of nitrate in the upper layer, and there is no diagenetic source for
nitrate in both layers:
Jijf03 = JNit (from Eq. 5-19) J^O3 = 0 (5-22)
Nitrate is present in sediments as dissolved substance, i.e., TTljN03 = K2jN03 = 0, making fdljN03 = fd2N03 = 1
(Equations 5-1 1 and 5-12): it also makes So meaningless, hence So = 0. Nitrate is removed by
denitrification in both oxic and anoxic layers with the carbon required for denitrification supplied by
carbon diagenesis. The reaction velocities in Equations 5-8 and 5-10 for nitrate may be expressed as:
No3 (5-23)
T — ^0
NO3 (5-24)
and the denitrification flux out of sediments as a nitrogen gas becomes:
2
(5-25)
%03 1 = reaction velocity for denitrification in Layer 1 at 20°C (m day"1)
%03,2 = reaction velocity for denitrification in Layer 2 at 20°C (m day"1)
QNOS = constant for temperature adjustment for KN03 : and KN03 2
JN2(g) = denitrification flux (g N m"2 day"1)
NO3j = total nitrate nitrogen concentration in Layer 1 (g N m"3)
NO32 = total nitrate nitrogen concentration in Layer 2 (g N m"3).
Once Equations 5-8 and 5-10 are solved for NO3j and NO32, the sediment flux of nitrate to the overlying
water, JaqN03, can be calculated using Eq. 5-9. The steady-state solution for nitrate showed that the
nitrate flux is a linear function of NO30 (Eq. 111-15 in D&F 1993): the intercept quantifies the amount of
5-13
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
ammonium in the sediment that is nitrified but not denitrified (thus releases as JaqN03), and the slope
quantifies the extent to which overlying water nitrate is denitrified in the sediment. It also revealed that
if the internal production of nitrate is small relative to the flux of nitrate from the overlying water, the
normalized nitrate flux to the sediment, - Jaq)N03/NO30, is linear in s for small s and constant for large s
(Section III-C in D&F 1993). For small s (~ 0.01 m day1), Hj is large (Eq. 5-13) so that oxic layer
denitrification predominates and JaqN03 is essentially zero independent of NO30 (Fig. III-4 in D&F 1993).
5.3.4 Phosphate Phosphorus
Phosphate is produced by the diagenetic breakdown of POP in the lower layer:
Jlf04 = 0 J2f04 = Jp (from Eq. 5-7) (5-26)
A portion of the liberated phosphate remains in the dissolved form and a portion becomes particulate
phosphate, either via precipitation of phosphate-containing minerals (Troup 1974), e.g., vivianite,
Fe3(PO4)2(s), or by partitioning to phosphate sorption sites (Lijklema 1980; Barrow 1983; Giordani and
Astorri 1986). The extent of particulate formation is determined by the magnitude of the partition
coefficients, TTljP04 and TT2jP04, in Equations 5-1 1 and 5-12. Phosphate flux is strongly affected by DO0,
the overlying water oxygen concentration. As DO0 approaches zero, the phosphate flux from the
sediments increases. This mechanism is incorporated by making KljP04 larger, under oxic conditions,
than TC2p04. In the model, when DO0 exceeds a critical concentration, (DO0)critP04, sorption in the upper
layer is enhanced by an amount AnP041:
(5-27)
When oxygen falls below (DO0)cntP04, then:
(5-28)
which smoothly reduces TCljP04 to TC2P04 as DO0 goes to zero. There is no removal reaction for phosphate
in both layers:
° (5-29)
Once Equations 5-8 and 5-10 are solved for PO4j and PO42, the sediment flux of phosphate to the
overlying water, JaqP04, can be calculated using Eq. 5-9.
5.3.5 Sulfide/Methane and Oxygen Demand
5-14
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
5.3.5.1 Sulfide. No diagenetic production of sulfide occurs in the upper layer. In the lower layer,
sulfide is produced by carbon diagenesis (Eq. 5-7) decremented by the organic carbon consumed by
denitrification (Eq. 5-25). Then:
(5-30)
a02C = stoichiometric coefficient for carbon diagenesis consumed by sulfide oxidation (2.6667 g O2-
equivalents per g C)
ao2,N03 = stoichiometric coefficient for carbon diagenesis consumed by denitrification (2.8571 g O2-
equivalents per g N).
A portion of the dissolved sulfide that is produced in the anoxic layer reacts with the iron to form
particulate iron monosulfide, FeS(s) (Morse et al. 1987). The particulate fraction is mixed into the oxic
layer where it can be oxidized to ferric oxyhydroxide, Fe2O3(s). The remaining dissolved fraction also
diffuses into the oxic layer where it is oxidized to sulfate. Partitioning between dissolved and particulate
sulfide in the model represents the formation of FeS(s), which is parameterized using partition
coefficients, K1>H2S and n2H2S, in Equations 5-1 1 and 5-12.
The present sediment model has three pathways for sulfide, the reduced end product of carbon
diagenesis: (1) sulfide oxidation, (2) aqueous sulfide flux, and (3) burial. The distribution of sulfide
among the three pathways is controlled by the partitioning coefficients and the oxidation reaction
velocities (Section V-E in D&F 1993). Both dissolved and particulate sulfide are oxidized in the oxic
layer, consuming oxygen in the process. In the oxic upper layer, the oxidation rate that is linear in
oxygen concentration is used (Cline and Richards 1969; Millero 1986; Boudreau 1991). In the anoxic
lower layer, no oxidation can occur. Then the reaction velocities in Equations 5-8 and 5-10 may be
expressed as:
j 2 .. aT ~ 20
» „,, (5-31)
Z'KMH2S,02
(5-32)
KH2s,di = reaction velocity for dissolved sulfide oxidation in Layer 1 at 20°C (m day"1)
KH2sPi = reaction velocity for particulate sulfide oxidation in Layer 1 at 20°C (m day"1)
9ms = constant for temperature adjustment for KH2S dl and KH2S pl
KMH2S 02 = constant to normalize the sulfide oxidation rate for oxygen (g O2 m"3).
5-15
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
The constant, KMH2S 02, which is included for convenience only, is used to scale the oxygen
concentration in the overlying water. At DO0 = KMH2S 02, the reaction velocity for sulfide oxidation rate
is at its nominal value.
The oxidation reactions in the oxic upper layer cause oxygen flux to the sediment, which exerts SOD. By
convention, SOD is positive: SOD = -Jaq02- Th£ SOD in the model consists of two components,
carbonaceous sediment oxygen demand (CSOD) due to sulfide oxidation and nitrogenous sediment
oxygen demand (NSOD) due to nitrification:
SOD = CSOD + NSOD = -^^H2Sl + aO2jm4-Jm (5-33)
jj
H2Sj = total sulfide concentration in Layer 1 (g O2-equivalents m"3)
ao2,NH4 = stoichiometric coefficient for oxygen consumed by nitrification (4.33 g O2 per g N).
Equation 4-29 is nonlinear for SOD because the RHS contains s (= SOD/DO0) so that SOD appears on
both sides of the equation: note that JNlt (Eq. 5-20) is also a function of s. A simple back substitution
method is used, as explained in Section 5.6.1.
If the overlying water oxygen is low, then the sulfide that is not completely oxidized in the upper layer
can diffuse into the overlying water. This aqueous sulfide flux out of the sediments, which contributes to
the chemical oxygen demand in the water column model, is modeled using
Jaqfl2S = s(fdifl2SH2Si ~ COD) (5-34)
The sulfide released from the sediment reacts very quickly in the water column when oxygen is available,
but can accumulate in the water column under anoxic conditions. The COD, quantified as oxygen
equivalents, is entirely supplied by benthic release in the water column model (Eq. 3-16). Since sulfide
also is quantified as oxygen equivalents, COD is used as a measure of sulfide in the water column in
Eq. 5-34.
5.3.5.2 Methane. When sulfate is used up, methane can be produced by carbon diagenesis and methane
oxidation consumes oxygen (DiToro et al. 1990). Owing to the abundant sulfate in the saltwater, only
the aforementioned sulfide production and oxidation are considered to occur in the saltwater. Since the
sulfate concentration in fresh water is generally insignificant, methane production is considered to
replace sulfide production in fresh water. In fresh water, methane is produced by carbon diagenesis in
5-16
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
the lower layer decremented by the organic carbon consumed by denitrification, and no diagenetic
production of methane occurs in the upper layer (Eq. 5-30):
(5-35)
The dissolved methane produced takes two pathways: (1) oxidation in the oxic upper layer causing
CSOD or (2) escape from the sediment as aqueous flux or as gas flux:
J2,CH4 = CSOD + Jaq,CH4 + -W) (5-36)
JaqcH4 = aqueous methane flux (g O2-equivalents m"2 day"1)
JcH4(g) = gaseous methane flux (g O2-equivalents m"2 day"1).
A portion of dissolved methane that is produced in the anoxic layer diffuses into the oxic layer where it is
oxidized. This methane oxidation causes CSOD in the freshwater sediment (DiToro et al. 1990):
- 20
CSOD = CSOD_
1 _ sech[KcH4'ecH4 ] I (5-37)
= mmimum2-KL-CH4sat-J2CH4, J^H4 (5-38)
CH4 , = 100 1 + 2 1.02420 - T (5-39)
sat ( 10 )
CSODmax = maximum CSOD occurring when all the dissolved methane transported to the oxic layer is
oxidized
KCH4 = reaction velocity for dissolved methane oxidation in Layer 1 at 20°C (m day"1)
0CH4 = constant for temperature adjustment for KCH4
CH4sat = saturation concentration of methane in the pore water (g O2-equivalents m"3).
The term, (h + H2)/10 where h and H2 are in meters, in Eq. 5-39 is the depth from the water surface that
corrects for the in situ pressure. Equation 5-39 is accurate to within 3% of the reported methane
solubility between 5 and 20°C (Yamamoto et al. 1976).
If the overlying water oxygen is low, the methane that is not completely oxidized can escape the sediment
into the overlying water either as aqueous flux or as gas flux. The aqueous methane flux, which
contributes to the chemical oxygen demand in the water column model, is modeled using (DiToro et al.
1990):
5-17
-------
-er"20
Model Report for Christina River Basin, Nutrient and DO TMDL
(5-40)
Methane is only slightly soluble in water. If its solubility, CH4sat given by Eq. 5-39, is exceeded in the
pore water, it forms a gas phase that escapes as bubbles. The loss of methane as bubbles, i.e., the
gaseous methane flux, is modeled using Eq. 5-36 with J2CH4 from Eq. 5-35, CSOD from Eq. 5-37, and
Jaq,cH4 from Eq. 5-40 (DiToro et al. 1990).
5.4 Silica
The production of ammonium, nitrate, and phosphate in sediments is the result of the mineralization of
POM by bacteria. The production of dissolved silica in sediments is the result of the dissolution of
particulate biogenic or opaline silica, which is thought to be independent of bacterial processes.
The depositional flux of particulate biogenic silica from the overlying water to the sediments is modeled
using Eq. 5-5. With this source, the mass-balance equation for particulate biogenic silica may be written
as:
H2^- - ~ SSi-H2 - W-PSi + Jm + J^ (5-41)
PSi = concentration of particulate biogenic silica in the sediment (g Si m"3)
SSl = dissolution rate of PSi in Layer 2 (g Si m"3 day"1)
JPSl = depositional flux of PSi (g Si m"2 day"1) given by Eq. 5-5
JDSl = detrital flux of PSi (g Si m"2 day"1) to account for PSi settling to the sediment that is not
associated with the algal flux of biogenic silica.
The processes included in Eq. 5-41 are dissolution (i.e., production of dissolved silica), burial, and
depositional and detrital fluxes from the overlying water. Equation 5-41 can be viewed as the analog of
the diagenesis equations for POM (Eq. 5-6). The dissolution rate is formulated using a reversible
reaction that is first order in silica solubility deficit and follows a Monod-type relationship in particulate
silica:
SVi = K^-Qt;: - - - •
s, s, s,
KSl = first order dissolution rate for PSi at 20°C in Layer 2 (day"1)
0Sl = constant for temperature adjustment for KSl
KMPSl = silica dissolution half-saturation constant for PSi (g Si m"3)
5-18
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Sisat = saturation concentration of silica in the pore water (g Si m"3).
The mass-balance equations for mineralized silica can be formulated using the general forms, Equations
5-8 and 5-10. There is no source/sink term and no reaction in the upper layer:
Ji# = Ki^ = ° (5-43)
In the lower layer, silica is produced by the dissolution of particulate biogenic silica, which is modeled
using Eq. 5-42. The two terms in Eq. 5-42 correspond to the source term and reaction term in Eq. 5-10:
T - V O^ ~ 20
J'y K - A.r,,'Uw
- 20 PSl
i PSi
A portion of silica dissolved from particulate silica sorbs to solids and a portion remains in the dissolved
form. Partitioning using the partition coefficients, TCljSi and n2 Sl, in Equations 5-1 1 and 5-12 controls the
extent to which dissolved silica sorbs to solids. Since silica shows similar behavior as phosphate in the
adsorption-desorption process, the same partitioning method as applied to phosphate (Section 5.3.4) is
used for silica. That is, when DO0 exceeds a critical concentration, (DO0)cnt Sl, sorption in the upper layer
is enhanced by an amount AnSll:
ni# = w2,s/-(A^,i) D0o > (^oW (5-46)
When oxygen falls below (DO0)cntSl, then:
"I* = *2y(A^ur°/(DOoW D0» * WW (5-47)
which smoothly reduces nljSi to K2 Sl as DO0 goes to zero.
Once Equations 5-8 and 5-10 are solved for Sij and Si2, the sediment flux of silicate the overlying water,
JaqSl, can be calculated using Eq. 5-9.
5.5 Sediment Temperature
All rate coefficients in the aforementioned mass-balance equations are expressed as a function of
sediment temperature, T. The sediment temperature is modeled based on the diffusion of heat between
the water column and sediment:
f = IT>(TW ~ ^ (5~48)
5-19
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
DT = heat diffusion coefficient between the water column and sediment (m2 sec"1)
Tw = temperature in the overlying water column (°C) calculated by Eq. 4-82.
The model application in D&F and Cerco and Cole (1993) used DT = 1.8 x 10~7 m2 sec"1.
5.6 Method of Solution
5.6.1 Finite-Difference Equations and Solution Scheme
An implicit integration scheme is used to solve the governing mass-balance equations. The finite
difference form of Eq. 5-8 may be expressed as:
0 = sd-Ct -d-Ct
W-Ct[ - — Ct[
(5-49)
where the primed variables designate the values evaluated at t+ (9 and the unprimed variables are those at
t, where 0 is defined in Eq. 4-82. The finite difference form of Eq. 5-10 may be expressed as:
0 = - KL
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
(2) Solve Eq. 5-51 for ammonium, nitrate, and sulfide/methane.
(3) Compute the SOD using Eq. 5-33.
(4) Refine the estimate of SOD: a root finding method (Brent's method in Press et al. 1986) is used
to make the new estimate.
(5) Go to (2) if no convergence.
(6) Solve Eq. 5-51 for phosphate and silica.
For the sake of symmetry, the equations for diagenesis, particulate biogenic silica and sediment
temperature are also solved in implicit form. The finite difference form of the diagenesis equation
(Eq. 5-6) may be expressed as:
GPOMJ = (GPOMJ + TT-W) f1 + Q-KPOM/VPOM™ + ^H ' (5-53)
I ^2 )( H1 )
The finite difference form of the PSi equation (Eq. 5-41) may be expressed as:
PSi' = PSi
*2
PSi + KM^ H,
2
-i
(5-54)
using Eq. 5-36 for the dissolution term, in which PSi in the Monod-type term has been kept at time level t
to simplify the solution. The finite difference form of the sediment temperature equation (Eq. 5-48) may
be expressed as:
-i
(5-55)
5.6.2 Boundary and Initial Conditions
The above finite difference equations constitute an initial boundary-value problem. The boundary
conditions are the depositional fluxes (JPOM,i and JPSi) and the overlying water conditions (Ct0 and Tw) as a
function of time, which are provided from the water column water quality model. The initial conditions
are the concentrations at t = 0, GPOMl(0), PSi(O), Ct^O), Ct2(0), and T(0), to start the computations.
Strictly speaking, these initial conditions should reflect the past history of the overlying water conditions
and depositional fluxes, which is often impractical because of lack of field data for these earlier years.
5-21
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
5-22
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
6 - EFDC WATER QUALITY MODEL CALIBRATION
The EFDC hydrodynamic and water quality model was used to determine the receiving water quality
conditions in the tidal and non-tidal streams in the Christina River Basin. Nutrient loads were input to
the EFDC model by means of linkage to the HSPF watershed loading models and the XP-SWMM CSO
simulation flow model. Flows and loads from over 100 NPDES facilities were also included in the
EFDC model.
6.1 Modeling Assumptions
The main objective of applying the EFDC model was to develop hydrodynamic and water quality
information for the primary stream channels throughout Christina River Basin. Specifically, it was
necessary to accurately understand the variability of flow throughout the stream network under variable
flow conditions. Major assumptions that contributed to the final approach taken included:
The waterbody was well mixed laterally and vertically, therefore a longitudinal one-dimensional
configuration was appropriate for the freshwater stream channels.
Thermal stratification was not likely due to the shallow and narrow characteristics of the creek,
thus temperature is not an important driving force for flow and transport.
Wind effects on flow and transport were not a critical factor due to the one-dimensional flow
pattern.
The impact of groundwater interaction on flow and transport was minimal during low flow
conditions, thus flow distribution can be obtained through directly balancing upstream and
downstream flow rates.
6.2 Model Configuration
The general procedure for application of the EFDC model to the Christina River Basin followed a
sequence of steps beginning with model configuration and continued through model execution of the
calibration time period. Model configuration involved the construction of the horizontal grid for the
waterbodies in the basin, interpolation of bathymetric data to the grid, construction of EFDC input files,
and compilation of the Fortran source code with appropriate parameter specification of array dimensions.
The model included 120 NPDES point-source discharges and 28 consumptive use water withdrawals.
The locations of the NPDES discharges are shown in Figure 6-1. Schematic drawings of the EFDC grid
configuration are presented in Appendix C. The locations of the NPDES discharges relative to EFDC
grid cells are shown in Figure C-l. The locations of the water withdrawals are shown in Figure C-2. The
EFDC model also included flows and loads from 38 CSO discharges and was linked to the HSPF
watershed loading models to incorporate nonpoint source flows and loads.
6-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
* 1 Fi1,,64 L >-"1"CS2
KILOMETERS
••I
24 6 8 10
Figure 6-1. Locations of NPDES discharges in Christina River Basin.
6-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
6.2.1 Segmentation
The numerical model domain includes the tidal Delaware River from Reedy Point on the south to Chester
on the north. Both the tidal and nontidal Christina River reaches are included in the model. The lower
Christina River is directly connected to the Delaware River. The nontidal Christina River is connected to
the tidal portion by a dam control structure at Smalleys Pond. The tidal Brandywine Creek is connected
to the tidal Christina River by means of a tidal inlet control structure. The tidal White Clay Creek is also
connected to the tidal Christina River via a tidal inlet control structure.
The basic equations in EFDC were solved using the finite-difference method. The grid was designed to
resolve velocity shears both axially and laterally, and at the same time allow a time step suitable for
efficient computation. Solutions to the hydrodynamics were obtained using a 60-second time step. The
spatial domain of the study area was divided into a grid of discrete cells. To achieve close conformance
of the grid to the estuary geometry, the cells in the Delaware River were represented using curvilinear
horizontal grid cells constructed using an orthogonal mapping procedure (Ryskin and Leal 1983) to form
a 2-D grid domain. The cells in the narrow tidal and nontidal streams were represented by a 1-D
Cartesian coordinate system (see Figure C-l). To obtain adequate resolution in the streams, longitudinal
cells were configured according to lengths ranging from 500 to approximately 1,000 meters. Cell widths
were adjusted according to estimated wetted stream channel widths under low-flow conditions.
Velocities were computed on the boundaries between cells, and temperature, salinity, and density were
computed at the center of each cell. The numerical grid consisted of 406 cells in the horizontal plane and
a single vertical layer. A single layer was chosen because the estuary and streams are well mixed,
thereby implying that stratification would not be an issue. In addition, field data available from STORET
and from Davis (1998) did not distinguish vertical sample depths.
6.2.2 Streamflow Estimation
Variable streamflow discharge was estimated using flows from the HSPF model for the calibration period
1994-1998. The streamflow was validated using observed daily average flows at several USGS stream
gages throughout the Christina River Basin (Senior and Koerkle, (2003a, 2003b, 2003c, and 2003d).
6.2.3 Atmospheric and Tidal Boundary Conditions
Atmospheric nutrient loads are typically divided into wet and dry deposition. Wet deposition is
associated with dissolved substances in rainfall. The settling of particulate matter during non-rainfall
events contributes to dry deposition. Observations of concentrations in rainwater are frequently
available, and dry deposition is usually estimated as a fraction of the wet deposition. The atmospheric
deposition rates reported in the Long Island Sound Study (HydroQual, 1991) and the Chesapeake Bay
6-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Model Study (Cerco and Cole, 1993) as well as information provided by DNREC for Lewes, Delaware,
were used to develop both dry and wet deposition loads for the EFDC model of the Christina River Basin
(see Tables 6-1 and 6-2). Meteorological information (i.e., atmospheric pressure, temperature, relative
humidity, wind speed and direction, rainfall, cloud cover, and solar radiation) was obtained from the
NOAA National Climatic Data Center weather station (WBAN 13781) at the New Castle County Airport
near Wilmington, Delaware.
Table 6-1. Atmospheric dry deposition rates used in Christina River Basin EFDC model.
Parameter
Refractory Part. Organic Carbon
Labile Part. Organic Carbon
Dissolved Organic Carbon
Dissolved Organic Phosphorus
Orthophosphate
Available Silica
Deposition Rate
(g/m2/day)
0.000387
0.000387
0.000773
0.000054
0.000019
0.000247
Parameter
Refractory Part. Organic Nitrogen
Labile Part. Organic Nitrogen
Dissolved Organic Nitrogen
Ammonia Nitrogen
Nitrate+Nitrite Nitrogen
Deposition Rate
(g/m2/day)
0.000530
0.000530
0.000771
0.000214
0.000393
Table 6.2. Atmospheric wet deposition concentrations used in Christina River Basin EFDC model.
Parameter
Refractory Part. Organic Carbon
Labile Part. Organic Carbon
Dissolved Organic Carbon
Dissolved Organic Phosphorus
Orthophosphate
Available Silica
Concentration
(mg/L)
0.325
0.325
0.650
0.045
0.016
0.0
Parameter
Refractory Part. Organic Nitrogen
Labile Part. Organic Nitrogen
Dissolved Organic Nitrogen
Ammonia Nitrogen
Nitrate+Nitrite Nitrogen
Concentration
(mg/L)
0.0
0.0
0.140
0.222
0.332
Tides were specified at the north and south boundaries in the Delaware River based on the astronomical
harmonic constants for the NOAA subordinate tide stations at Reedy Point, Delaware (south boundary)
and Chester, Pennsylvania (north boundary). The predicted tides from the harmonic constants do not
include any low-frequency influences due to storms or regional low-pressure conditions (NOAA, 1998)..
The specification of boundary conditions was required at the model north and south interface with the
Delaware River. The EFDC water quality model accommodates 21 boundary variables, each specified in
an individual time-series data file of concentrations. Advective boundary conditions in the Christina
River model were of the "upwind" type. Evaluation of the boundary concentration depended on the
direction of flow at the boundary. When flow was out of the model, the boundary concentration was
6-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
assigned the concentration in the model cell immediately upstream of the boundary. When the tidal flow
was into the model, the boundary concentration was assigned a specified, time-varying value
representative of conditions outside the model domain. To estimate recirculation at the boundary near
the time of flow reversal from outgoing to incoming tide, the last outgoing concentration at the boundary
is used as the incoming concentration for a certain amount of time specified by the user. This
concentration linearly approaches the specified outside boundary concentration over that time period.
For the Christina River model, the recirculation time interval was specified as 60 minutes based on
experience gained from previous water quality model applications of the EFDC model.
Delaware River boundary conditions for salinity, temperature, total suspended sediment, algae, organic
carbon, dissolved oxygen, nitrogen, phosphorus, silica, and fecal coliform bacteria were specified based
on available STORET data at stations in the Delaware River. The boundary time-series were created
using observations that were averaged by month over the simulation period. If data for a parameter were
not available for any given month, then the long-term average (over the period 1988-1998) for that month
was used instead.
6.2.4 Initial Conditions
Initial conditions for freshwater streams in the EFDC model at the starting time of October 1, 1994, were
estimated using the simulated flows and nutrient concentrations calculated by the HSPF model. Initial
water quality concentration conditions in the tidal Delaware River and tidal Christina River were
estimated using the ending conditions from the 1995 low-flow validation run (September 30, 1995).
These initial conditions allow the model to begin its simulation at a stable numeric state. The impacts of
initial conditions diminish quickly with time.
6.2.5 Point and Nonpoint Source Representation
External flows and loads of nutrients and oxygen demand were divided into four categories: (1) nonpoint
source loads (i.e., diffuse sources) including tributary sources and groundwater sources, (2) point-source,
(3) water withdrawals, and (4) atmospheric deposition. Nonpoint source loads were carried by
freshwater flows and groundwater entering the main stream reaches. Point-source loads were discharges
from the NPDES facilities and CSOs in the study area. Consumptive use water withdrawals were
removed from the model system at the appropriate grid cell. Atmospheric loads were transfers from the
atmosphere to the water surface via rainfall (wet deposition) and other processes (dry deposition).
Atmospheric deposition is not a significant source in the narrow stream channels, but may be more
important in the open estuary waterbodies in the lower Christina River and Delaware River because of
the larger water surface area in those regions.
6-5
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Nonpoint sources were estimated by the delineation of subbasins and land use categories in the HSPF
watershed loading models. The nonpoint source loads generated by the watershed models provided
predictive nutrient loads to the receiving waters reflective of variable meteorological (rainfall-runoff)
characteristics.
Discharge Monitoring Records (DMRs) for various NPDES point sources in the Brandywine Creek
watershed were provided in hard copy form by the Brandywine Valley Association. Other DMRs were
provided in electronic format by PADEP and DNREC. The hard-copy data were keypunched and the
electronic data were reformatted into a database file for use in developing point source loads for the
water quality model. A list of all 120 NPDES discharges included in the model is given in Appendix C
(Table C-l). The August 1997 field monitoring study (Davis 1998) included seven NPDES discharges
that were monitored for flow and water quality parameters (see Appendix C, Figure C-3). Loading
values for the various water quality constituents were computed based on the flow rates and
concentrations provided on the DMRs or measured during the August 1997 study.
The NPDES discharges included single residence discharges (SRD) that are not required to submit DMR
data. For purposes of model calibration, it was assumed that these SRD discharges operated at their
permit discharge limits. Characteristic concentrations for the various water quality parameters were
assigned to the NPDES source based on the type of discharge, and the loading in kg/day for each
constituent was computed for input to the EFDC model. The characteristic effluent concentrations used
for this study are listed in Appendix C, Table C-2, and the characteristic effluent parameter ratios are
listed in Table C-3. The characteristic effluent concentrations and parameter ratios were derived from
effluent monitoring data collected by Davis (1998) in August 1998 and from literature values reported in
the Technical Guidance Manual for Developing TMDLs (USEPA, 1995).
For model calibration, a time-series of monthly average loads for the 1994-1998 simulation period was
developed for nutrients, dissolved oxygen, CBOD, and total suspended solids for each NPDES point
source based on available discharge monitoring records (DMRs). The methodology for estimating the
various species of nitrogen and phosphorus is outlined in Table 6-3 and was described in the low-flow
modeling report (USEPA, 2000). The ratios for converting CBOD5 to organic carbon for the model were
determined based on data collected during a special study conducted in August-September 1999 from
several of the larger WWTPs in the basin (USEPA, 2000).
6-6
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Table 6-3. Methodology for developing EFDC point source loads from DM R data.
Water Quality Parameter
CBQD-5-day
CBOD-ultimate
Total organic carbon
Dissolved organic carbon
Refractory particulate organic carbon
Labile particulate organic carbon
Total phosphorus
Total organic phosphorus
Refractory particulate organic phosphorus
Labile particulate organic phosphorus
Dissolved organic phosphorus
Total orthophosphate
Total nitrogen
Nitrite nitrogen
Total organic nitrogen
Refractory particulate organic nitrogen
Labile particulate organic nitrogen
Dissolved organic nitrogen
Ammonia nitrogen
Nitrate nitrogen
Unavailable biogenic silica
Dissolved available silica
Chemical oxygen demand
Dissolved oxygen
Total active metal
Fecal coliform bacteria
EFDC Code
TOC
DOC
RPOC
LPOC
RPOP
LPOP
OOP
PO4T
RPON
LPON
DON
NH3
NO3
SUU
SAA
COD
DOO
TAM
FCB
Calculation
CBOD5 = BQD5 * (CBOD5:BOD5 ratio)
CBODu = CBOD5* (CBODu:CBOD5 ratio)
TOC = CBODu * (TOC:CBODu ratio)
DOC = TOC * (DOC:TOC ratio)
0.5 * (TOC - DOC)
0.5 * (TOC - DOC)
If TP not reported on DMR, use default TP from Table C-2
TOP = TP - (TP * (OPO4TP ratio))
0.25 TOP
0.25 TOP
0.50 TOP
TP * (OPO4:TP ratio)
TN = NH3-N * (TN:NH3 ratio)
NO2-N = NH3-n * (NO2:NO3 ratio)
TON = TN - NO2-N - NO3-N - NH3-N
0.25 TON
0.25 TON
0.50 TON
reported on DMR (or use default NH3-N from Table C-2)
NO3-N = NH3 * (NO3:NH3 ratio)
0.10 mg/L (default value)
1.00 mg/L (default value)
9.6 * CBOD5
reported on DMR (or use default value from
Table C-2)
0.0 (not simulated)
reported on DMR (or use default value from
Table C-2)
CSO flows were estimated using XP-SWMM and were provided by the City of Wilmington. Nutrient
loads from CSO outfalls were estimated using the XP-SWMM flow rates and event mean concentrations
based on storm event monitoring conducted by the City of Wilmington and Delaware DNREC (see
Tables 2-1 a, b, and c).
6.2.6 Time Step and Simulation Duration
The EFDC model was executed at a time step of 60 seconds and the calibration simulated four
consecutive water years covering the period from October 1, 1994 to October 1, 1998. A listing of the
key EFDC input data files is presented in Appendix D.
6-7
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
6.3 Model Calibration Results
Model calibration involves the adjustment of certain model input quantities in an attempt to achieve a
specified level of model performance. An extensive set of field data were gathered, processed, and
displayed for modeling hydrodynamics and water quality transport in the Christina River Basin. The data
set included database files containing more than 40,000 records at about 200 stations scattered
throughout the interior of the basin as well as in the Delaware River itself. This section presents the
results of the calibration of the EFDC hydrodynamic and water quality model. Parameters considered for
calibration include flow rate and a suite of water quality parameters including nitrogen, phosphorus,
carbon, and dissolved oxygen.
6.3.1 Tide Elevation and Phase
Calibration of the model with respect to water surface elevation was accomplished by analysis of
observed and model predicted time-series data at two interior tide stations. For tidal waters, least squares
harmonic analysis is the most commonly utilized procedure (Oey, Mellor, and Hires, 1985; Cheng et al.,
1993; Shen et al., 1999). Tide elevation data were obtained from the USGS tide stations on the Christina
River at the Port of Wilmington near the mouth and at Newport about 7.0 miles upstream of the mouth.
These data were compared with surface elevations computed by the model at cell 56,13 (Port of
Wilmington) and 45,13 (Newport). The time-series of tide elevations for the month of August 1997 for
both the field data and model results were subjected to a harmonic analysis. The five most important
astronomical harmonic constituents (M2, S2, N2, Kl, and Ol) were computed for both the field data and
model simulation results. The harmonic analysis results, shown in Table 6-4, indicate the model is in
good agreement with the measured tide data for both amplitude and phase. The model-data amplitudes
for the M2 harmonic constituent agree within 5 cm (6%) and the phases agree to within 4 degrees (3%).
Time-series graphs (Appendix C, Figure C-4) of the observed and model tide elevations at both the Port
of Wilmington and Newport covering a 15-day period (August 1 - 15, 1997) provide a visual means of
assessing the skill of the model in simulating tidal elevations. The model tides are forced at the north and
south boundaries in the Delaware River based on the NOAA predictions at the Reedy Point, DE, and
Chester, PA, subordinate stations (NOAA, 1998).
6-8
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Table 6-4. Harmonic analysis of tides at Port of Wilmington and Newport
Harmonic Constant
M2 - observed
M2 - model
Difference
S2 - observed
S2 - model
Difference
N2 - observed
N2 - model
Difference
K1 - observed
K1 - model
Difference
O1 - observed
O1 - model
Difference
Port of Wilmington
Amplitude (m)
0.7594
0.7135
0.0459
0.0894
0.1001
-0.0107
0.1271
0.1383
-0.0112
0.0802
0.0633
0.0169
0.0626
0.0546
0.0080
Phase (degrees)
130.382
134.180
-3.798
20.621
30.806
-10.185
323.153
336.181
-13.028
174.059
178.335
-4.276
316.879
326.765
-9.886
Newport
Amplitude (m)
0.6901
0.6768
0.0133
0.0900
0.0890
0.0010
0.1275
0.1240
0.0035
0.0615
0.0606
0.0009
0.0546
0.0514
0.0032
Phase (degrees)
153.634
155.560
-1.926
36.374
59.180
-22.806
345.054
3.603
-18.549
184.740
190.948
-6.208
332.386
337.937
-5.551
6.3.2 Water Depth and Stream Velocity
Measurements of flow, water depth, and stream velocity were made at eight locations during the August
1997 field survey (Davis, 1998). The field measurements were made on the following dates: East Branch
Brandywine Creek (08/12 - 08/14/1997), West Branch Brandywine Creek (08/19 - 08/20/1997), West
Branch Red Clay Creek (08/05 - 08/07/1997 and 08/12 - 08/14/1997), and East Branch White Clay Creek
(08/26 - 08/28/1997). A comparison of these measurements with the model results at the appropriate grid
cell (I,J) location is given in Table 6-5.
Table 6-5. Model-data comparison of velocity, flow, and geometry (August 1997 data).
Stream Reach
East Branch Brandywine Creek
East Branch Brandywine Creek
EFDC
Cell
54,61
54,56
Velocity (fps)
Field
0.33
0.85
EFDC
0.48
0.56
Depth (ft)
Field
0.82
1.02
EFDC
0.87
1.11
Flow (cfs)
Field
14.5
34.3
EFDC
25.6
34.5
Channel Width (ft)
Field
53.6
39.6
EFDC
52.5
52.5
i/Vest Branch Brandywine Creek
i/Vest Branch Brandywine Creek
19,79
26,79
0.40
0.41
0.41
0.36
1.09
0.70
0.94
0.82
9.5
32.0
14.9
32.9
45.0
111.5
42.6
111.5
East Branch White Clay Creek
East Branch White Clay Creek
19,31
19,29
0.44
0.42
0.40
0.41
0.93
0.85
0.96
0.86
5.30
7.35
5.33
7.34
13.0
20.6
12.8
20.3
West Branch Red Clay Creek
i/Vest Branch Red Clay Creek
29,43
33,43
0.35
0.49
0.44
0.52
0.75
0.90
0.78
0.94
3.55
5.45
3.35
4.92
13.5
12.4
13.5
12.4
6-9
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
6.3.3 Sediment Oxygen Demand and Benthic Nutrient Flux Rates
The need for a predictive benthic sediment processes model for water quality modeling projects has been
apparent for some time. When using a water quality model for management scenario analysis, one of the
biggest sources of uncertainty involves what to use for the future sediment flux rates after a proposed
management control has been implemented. The predictive sediment submodel in EFDC helps address
this uncertainty with two fundamental capabilities: (1) the ability to predict effects of management
alternatives on sediment-water exchange processes and (2) the ability to predict the time scale for
alterations in the sediment-water exchange processes. To meet these requirements, a predictive sediment
process model was incorporated into the EFDC model framework and was based on DiToro and
Fitzpatrick (1993). The sediment submodel is driven by net settling of organic matter from the water
column to the sediments. In the benthos, the sediment submodel simulates the decay (diagenesis) of
organic matter, which produces oxygen demand and inorganic nutrients. Oxygen demand takes three
paths out of the sediments: (1) export to the water column as chemical oxygen demand, (2) oxidation at
the sediment-water interface as sediment oxygen demand, or (3) burial to a deep, inactive sediment layer.
The inorganic nutrients produced by diagenesis can take two pathways out of the bottom sediment: (1)
release back to the overlying water column or (2) burial to the deep, inactive sediment layer.
In the predictive sediment submodel, benthic sediments are represented as two layers with a total depth
of 10 cm. The upper benthic layer is in contact with the water column and may be oxic or anoxic
depending on the dissolved oxygen concentration in the water. The lower benthic layer is permanently
anoxic. The thickness of the upper benthic layer is determined by the penetration of oxygen into the
sediments, and at its maximum thickness, the oxic layer depth is a small fraction of the total thickness.
The sediment submodel consists of three basic processes:
• Particulate organic matter settles from the water column to the sediments. Because of the
negligible thickness of the upper benthic layer, deposition proceeds from the water column
directly to the lower anoxic layer.
• Within the lower layer, organic matter is subject to decay (diagenesis).
• The flux of substances produced by diagenesis moves to the upper benthic layer, to the water
column, and to the deep, inactive benthic layer (burial). The flux portion of the sediment
submodel is the most complex. The computation of flux requires consideration of
(1) reactions in both benthic layers, (2) sedimentation from the upper to lower benthic layer
as well as from the lower benthic layer to the deep inactive sediments, (3) particle mixing
between layers, (4) diffusion between layers, and (5) mass transfer between the upper layer
and the water column.
6-10
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Very limited field data were available during the calibration period to verify the flux rates computed by
the predictive sediment submodel. SOD rates were measured in July and August 1996, at three locations
in the tidal Christina River and Brandywine Creek. An SOD rate of 0.5 g/m2/day was used in the tidal
Delaware River in another model study commissioned by the Delaware River Basin Commission
(DRBC) and was used as the basis for comparison to predicted SOD rates from this study. The simulated
SOD rates were converted to rates at 20°C and are compared with the measured data in Table 6-6. The
relative errors were less than 13% at all locations, which is considered to be a very good model-data skill
assessment.
Table 6-6. Model-data comparison of sediment oxygen demand rates (g/m2/day)
Location
Christina River at 1-495 bridge
Christina River at Newport, Rt. 141
bridge
Brandywine Creek, 0.6 mi. from mouth
Delaware River (from HydroQual study)
Sampling
Date
Aug 12, 1996
Jul 10, 1996
Aug 12, 1996
-
Monitored
SOD at 20°C
0.81
1.67
1.23
0.50
EFDC Model
SOD at 20°C
0.91
1.56
1.19
0.46
Relative
Error
12.9%
6.5%
3.4%
8.8%
6.3.4 Water Quality Results
Each field observation was collected at an instant in time and at a single point in space. Time scales
realistically represented in the EFDC model were determined by time scales of primary forcing functions:
60-second tidal hydrodynamics time-step, hourly meteorological inputs, monthly ocean boundary
conditions, daily nonpoint source loads, monthly point source loads, daily CSO loads, constant
atmospheric dry deposition, and hourly atmospheric wet deposition during rain events. The minimum
model spatial scales were determined by the size of the grid cells, ranging from 500 to about 1,000
meters in the longitudinal direction along the streams. The disparity in the temporal and spatial scales
between the model and prototype, especially for the nonpoint and point source loads, meant that
individual observations may not be directly comparable with model prediction at a specific time in a
given model grid cell.
Model-data comparisons were be made qualitatively (time-series graphics) and quantitatively (model-
data statistics). The time-series graphics are provided in Appendix A and cover the entire 4-year
calibration period beginning Oct 1, 1994 and continuing to Oct 1, 1998. The model-data time-series
comparison graphics were made at 27 monitoring locations on various streams in the study area (see the
map in Appendix A, Figure A-0).
The graphical model-data time-series comparisons in Appendix A provide a qualitative evaluation of
model performance. A seasoned modeler can examine the plots and form an experience-based judgment
6-11
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
on the status of model calibration and verification. The model-data statistical analysis provides a
different perspective on model-data comparison that numerically quantifies the state of model
calibration/verification (sometimes referred to as model "skill assessment").
Although numerous methods exist for analyzing and summarizing model performance, there is no
consensus in the modeling community on a standard analytical suite. A set of basic statistical methods
were used to compare model predictions and sampling observations which included the mean error
statistic, the absolute mean error, the root-mean-square error, and the relative error. Statistics for the
observations and model predictions were calculated over the period Oct 1, 1994 to Oct 1, 1998 at 24
monitoring locations in the Christina River Basin (see Table 6-7 and the map in Appendix A,
Figure A-0).
Table 6-7. Monitoring stations used for time-series model-data statistical analysis
Station
104011
104021
104051
WQN0105
103041
103061
103031
103011
WQN0149
105031
105011
105131
105071
106191
106141
106031
106021
106011
106291
106281
BCWB05
BCWB04
BCEB02
RCWB02
EFDCgrid cell (I,J)
54,20
54,23
54,32
54,36
43,38
48,52
43,30
43,24
19,18
21,18
41,18
31,34
31,40
14,13
22,13
32,13
47,13
53,13
55,13
43,55
27,79
21,79
54,55
29,43
Stream and Location
Brandywine Creek at Brandywine Park
Brandywine Creek at Road 279
Brandywine Creek at Smith Bridge
Brandywine Creek
Red Clay Creek at Ashland, DE
Burroughs Run at Rt. 241
Red Clay Creek at Woodale, DE
Red Clay Creek at Stanton, DE
White Clay Creek
White Clay Creek at Road 329 near Thompson
White Clay Creek at Rt. 7 in Stanton
Muddy Run at Road 303
Mill Creek at Road 282
Christina River above Newark at Rt. 273
Christina River at Road 26
Christina River at Smalleys Pond
Christina River at Rt. 141 in Newport
Christina River at US Rt. 1 3
Christina River at RR Bridge near Port of Wilmington
Little Mill Creek at Atlantic Avenue
Brandywine Creek West Branch at Modena, PA
Brandywine Creek West Branch at Coatesville, PA
Brandywine Creek East Branch below Downingtown, PA
Red Clay Creek West Branch near Kennett Square, PA
6-12
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
6.3.4.1 Mean Error Statistic. The mean error between model predictions and observations is defined in
Eq. 6- 1 . A mean error of zero is ideal. A non-zero value is an indication that the model may be biased
toward either over- or underprediction. A positive mean error indicates that on average the model
predictions are less than the observations. A negative mean error indicates that on average the model
predictions are greater than the observed data. The mean error statistic may give a false ideal value of
zero (or near zero) if the average of the positive deviations between predictions and observations is about
equal to the average of the negative deviations in a data set. Because of that possibility, it is never a good
idea to rely solely on this statistic as a measure of performance. Instead, it should be used in tandem with
the other statistical measures that are described in this section.
P S (O - P)
- -L—>
where:
E = mean error
O = observation, aggregated by month and over the water column
p = model prediction, aggregated by month and over vertical layers
n = number of observed-predicted pairs
6.3.4.2 Absolute Mean Error Statistic. The absolute mean error between model predictions and
observations is defined in Eq. 6-2. An absolute mean error of zero is ideal. The magnitude of the
absolute mean error indicates the average deviation between model predictions and observed data.
Unlike the mean error, the absolute mean error cannot give a false zero.
_ S \0 - P\
LabS -- (6-2)
where:
Eabs = absolute mean error.
6.3.4.3 Root-Mean-Square Error Statistic. The root-mean-square error (Ermv) is defined in Eq. 6-3. A
root-mean-square error of zero is ideal. The root-mean-square error is an indicator of the deviation
between model predictions and observations. The Erms statistic is an alternative to (and is usually larger
than) the absolute mean error.
P
rms
\
(11-3)
where:
6-13
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
= root-mean-square error
6.3.4.4 Relative Error Statistic. The relative error between model predictions and observations is
defined in Eq. 6-4. A relative error of zero is ideal. The relative error is the ratio of the absolute mean
error to the mean of the observations and is expressed as a percent.
_
rd
s o
(6-4)
where:
Erel = relative error.
6.3.4.5 Statistics Results. A summary of the error statistics for eight key water quality parameters of the
Christina River Basin model calibration simulation is given in Table 6-8. The relative error statistic
permits comparisons between the various water quality substances. Temperature and dissolved oxygen
were the parameters with the smallest relative error. The results for temperature indicate a relative error
of about 5.7%, and the relative error for dissolved oxygen was less than 8.4%. The relative error for total
nitrogen was less than 15%, ammonia nitrogen was less than 40%, total phosphorus was about 30%, total
organic carbon was less than 17%, and dissolved organic carbon was about 27%.
Table 6-8. Statistical summary of EFDC water quality model 1994-1998 calibration results
Parameter
Dissolved Oxygen (mg/L)
Total Organic Carbon (mg/L)
Diss. Organic Carbon (mg/L)
Total Nitrogen (mg/L)
Ammonia Nitrogen (mg/L)
Nitrate Nitrogen (mg/L)
Total Phosphorus (mg/L)
Temperature (degC)
Mean Error
-0.2632
0.4766
0.9848
0.1644
0.0142
0.0329
0.0230
-0.2813
Absolute Mean Error
0.7732
1.0698
1.4431
0.4340
0.0276
0.3405
0.0357
0.7487
RMS Error
1.1747
1.9450
2.1061
0.7479
0.0553
0.5089
0.0803
1.2794
Relative Error
8.35%
16.65%
26.94%
14.41%
39.91%
35.30%
30.20%
5.74%
Number of
Samples
859
820
818
778
774
812
785
862
According to the Technical Guidance Manual for Performing Waste Load Allocations (USEPA 1990),
acceptable relative error statistic criteria are 15% for dissolved oxygen and 45% for nutrient parameters
(nitrogen, phosphorus, and carbon). The overall relative error statistics for the Christina River model
were 8.35% for dissolved oxygen, 14.4% for total nitrogen, 30.2% for total phosphorus, and 16.7% for
total organic carbon. Since the relative error statistics for the Christina River EFDC water quality model
meet the general guidance criteria published in USEPA (1990), and the model is considered acceptable
for conducting TMDL allocation analyses.
6-14
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
7 - REFERENCES
Ambrose, R.B., T.A. Wool, and J.L. Martin. 1993. The water quality analysis and simulation program, WASPS:
Part A, model documentation version 5.1. U.S. Environmental Protection Agency, Environmental Research
Laboratory, Athens, GA, 210 pp.
Andrews, D.G., and M.E. Mclntyre. 1978. An exact theory for of nonlinear waves on a Lagrangian flow. J. Fluid
Mech. 89:609-646.
Banks, R.B. and F.F. Herrera. 1977. Effect of wind and rain on surface re aeration. ASCE J. of the Environ. Engr.
Div. 103(EE3):489-504.
Bennett, A.F. 1976. Open boundary conditions for dispersive waves. J. Atmos. Sci. 32:176-182.
Bennett, A.F., and P.C. Mclntosh. 1982. Open ocean modeling as an inverse problem: tidal theory. J. Phys. Ocean.
12:1004-1018.
Bennett, J.R., and A.H. Clites. 1987. Accuracy of trajectory calculation in a finite-difference circulation model. J.
Comp. Phys. 68:272-282.
Blumberg, A.F., and L.H. Kantha. 1985. Open boundary condition for circulation models. J. Hydr. Engr.
111:237-255.
Blumberg, A.F., and G.L. Mellor. 1987. A description of a three-dimensional coastal ocean circulation model.
Three-Dimensional Coastal Ocean Models, Coastal andEstuarine Science, Vol. 4, ed. N.S. Heaps, pp. 1-19.
American Geophysical Union.
Boni, L., E. Carpene, D. Wynne, and M. Reti. 1989. Alkaline phosphatase activity inProtogonyaulax tamarensis.
J. of Plankton Research 11(5):879-885.
Boudreau, B.P. 1991. Modelling the sulfide-oxygen reaction and associated pH gradients in porewaters.
Geochimica et Cosmochimica Acta, 55:145-159.
Bowie, G.L., W.B. Mills, D.B. Porcella, C.L. Campbell, J.R. Pagenkopf, G.L. Rupp, K.M. Johnson, P.W.H.Chan,
S.A. Gherini, and C.E. Chamberlin. 1985. Rates, constants and kinetics formulations in surface water quality
modeling (2nd edition). EPA/600/3-85/040, U.S. Environmental Protection Agency, Environmental Research Lab.,
Athens, GA, 455 pp.
Burban, P.Y., W. Lick, and J. Lick. 1989. The flocculation of fine-grained sediments in estuarine waters. J.
Geophys. Res. 94:8323-8330.
Burban, P.Y., Y.J. Xu, J. McNeil, and W. Lick. 1990. Settling speeds offices in fresh and seawater. J. Geophys.
Res. 95:18,213-18,220.
Carritt, D.E. and S. Goodgal. 1954. Sorption reactions and some ecological implications. Deep-Sea Research
1:224-243.
Caupp, C.L., Brock, J.T., and Runke, H.M. 1991. Application of the dynamic stream simulation and assessment
model (DSSAM III) to the Truckee River below Reno, Nevada: Model formulation and program description. Report
prepared by Rapid Creek Water Works for Nevada Division of Environmental Protection, Carson City and Washoe
County Dept. Of Comprehensive Planning, Reno, NV.
Cerco, C.F., and T. Cole. 1993. Three-dimensional eutrophication model of Chesapeake Bay. J. Environ. Engnr.
119:1006-1025.
7-1
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Cheng, R.T., V. Casulli, and J.W. Gartner. 1993. Tidal, residual, intertidal mudflat (TRIM) model and its
applications to San Francisco Bay, California. Estuarine, Coastal, and Shelf Science 36: 235-280.
Chrost, RJ. and J. Overbek. 1987. Kinetics of alkaline phosphatase activity and phosphorus availability for
phytoplankton and bacterioplankton in Lake PluBsee (North German eutrophic lake). Microbial Ecology 13:229-
248.
Cline, J.D. and F.A. Richards. 1969. Oxygenation of hydrogen sulfide in seawater at constant salinity, temperature
and pH. Environmental Science & Technology, 3(9):838-843.
Cole, T.M., and E.M. Buchak. 1994. CE-QUAL-W2: A two-dimensional laterally averaged, hydrodynamic and
water quality model, Version 2.0, Report ITL-93-7. U. S. Army Corps of Engineers, Waterway Experiment Station,
Vicksburg, MS.
Davis, J. 1998. Measurement of community photosynthesis and respiration rates for selected reaches of the
Christina Watershed. Prepared for Pennsylvania Department of Environmental Protection and Delaware Department
of Natural Resources and Environmental Control, March 1998.
Diaz, RJ. and R. Rosenberg. 1995. Marine benthic hypoxia: a review of its ecological effects and the behavioural
responses of benthic macrofauna. Oceanography and Marine Biology: an Annual Review 33:245-303.
DiToro, D.M. 1980. Applicability of cellular equilibrium and Monod theory to phytoplankton growth kinetics.
Ecological Modelling 8:201-218.
DiToro, D.M., P.R. Paquin, K. Subburamu, and D.A. Gruber. 1990. Sediment oxygen demand model: methane and
ammonia oxidation. ASCE J. Environ. Engr. 116(5):945-986.
DiToro, D.M., and J.F. Fitzpatrick. 1993. Chesapeake Bay sediment flux model. Contract Report EL-93-2.
Prepared by HydroQual, Inc. for U. S. EPA Chesapeake Bay Program, U. S. Army Engineer District, Baltimore,
MD, and U.S. Army Engineer Waterways Exp. Station.
Frick, W.E. 1984. Non-empirical closure of the plume equations, Atmos. Environ. 18:653-662.
Froelich, P.N. 1988. Kinetic control of dissolved phosphate in natural rivers and estuaries: a primer on the
phosphate buffer mechanism. Limnol. and Oceanogr. 33(4, part 2):649-668.
Galperin, B., L.H. Kantha, S. Hassid, and A. Rosati. 1988. A quasi-equilibrium turbulent energy model for
geophysical flows. J. Atmos. Sci. 45:55-62.
Genet, L.A., D.J. Smith, and M.B. Sonnen, M.B. 1974. Computer program documentation for the dynamic estuary
model. Prepared for U.S. Environmental Protection Agency, Systems Development Branch, Washington, D.C.
Grant, W.D., and O.S. Madsen. 1986. The continental-shelf bottom boundary layer. Annual Review of Fluid
Mechanics, ed. M. Van Dyke et al., pp. 365-306, Annual Review, Inc.
Giordani, P. and M. Astorri. 1986. Phosphate analysis of marine sediments. Chemistry in Ecology 2:103-112.
Hamrick, J.M. 1992a. A Three-Dimensional Environmental Fluid Dynamics Computer Code: Theoretical and
Computational Aspects, Special Report 317. The College of William and Mary, Virginia Institute of Marine Science.
63 pp.
Hamrick, J.M. 1992b. Estuarine environmental impact assessment using a three-dimensional circulation and
transport model. Estuarine and Coastal Modeling, Proceedings of the 2nd International Conference, M. L. Spaulding
et al, Eds., American Society of Civil Engineers, New York, 292-303.
7-2
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Hamrick, J. M. 1994: Linking hydrodynamic and biogeochemcial transport models for estuarine and coastal waters.
Estuarine and Coastal Modeling, Proceedings of the 3rd International Conference, ed. M.L. Spaulding et al.,
pp. 591-608. American Society of Civil Engineers, New York.
Hamrick, J.M. 1996. A User's Manual for the Environmental Fluid Dynamics Computer Code (EFDC), Special
Report 331. The College of William and Mary, Virginia Institute of Marine Science. 234 pp.
Hamrick, J.M., and T.S. Wu. 1996. Computational design and optimization of the EFDC/HEM3D surface water
hydrodynamic and eutrophication models. Computational Methods for Next Generation Environmental Models,
ed. G. Delich. Society of Industrial and Applied Mathematics, Philadelphia.
Homer, R.R., E.B. Welch, M.R. Seeley, and J.M. Jacoby. 1990. Responses of periphyton to changes in current
velocity, suspended sediment, and phosphorus concentration. Freshwater Biol. 24(2):215-232.
HydroQual. 1991. Draft Water Quality Modeling Analysis of Hypoxia in Long Island Sound. Prepared for
Management Committee Long Island Sound Estuary Study and New England Interstate Water Pollution Control
Commission. Prepared by HydroQual, Inc., Mahwah, NJ. Job Number: NENG0012. July 1991.
Jirka, G. H., and R.L. Doneker. 1991. Hydrodynamic classification of submerged single-port discharges, J. ofliydr.
Engr. 117:1095-1112.
Jirka, G.H., and P.J. Akar. 1991. Hydrodynamic classification of submerged multiport-diffuser discharges, J. of
Hydr. Engr. 117:1113-1128.
Kang, I.S., and L.G. Leal. 1992. Orthogonal grid generation in a 2D domain via the boundary integral technique. J.
Com p. Phys. 102:78-87.
Kremer, J.N. and S.W . Nixon. 1978. A coastal marine ecosystem: simulation and analysis. Ecological Studies 24,
Springer-Verlag, New York. 217 pp.
Lebo, M.E. 1991. Particle-bound phosphorus along an urbanized coastal plain estuary. Marine Chemistry 34:225-
246.
Lee, J.H.W ., and V. Cheung. 1990. Generalized Lagrangian model for buoyant jets in a current. J. Environ. Engrg.
116:653-662.
Lijklema, L. 1980. Interaction of ortho-phosphate with iron (III) and aluminum hydroxides. Environmental Science
& Technology, 14(5):537-541.
Madsen, P.A., and J. Larsen. 1987. An efficient finite-difference approach to the mild-slope equation. Coastal Engr.
11:329-351.
Mancini, J.L. 1983. A method for calculating effects, on aquatic organisms, of time varying concentrations. Water
Research 17(10):1355-1362.
Matisoff, G. 1982. Mathematical models of bioturbation, ed. P.L. McCall and M.J.S. Tevesz, pp. 289-330. Animal-
Sediment Relations: The Biogenic Alteration of Sediments, Plenum Press, NY.
Mellor, G.L., and T. Yamada. 1982. Development of a turbulence closure model for geophysical fluid problems.
Rev. Geophys. Space Phys., 20:851-875.
Millero, F.J. 1986. The thermodynamics and kinetics of the hydrogen sulfide system in natural waters. Marine
Chemistry 18:121-147.
7-3
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Mobley, C.D., and R.J. Stewart. 1980. On the numerical generation of boundary-fitted orthogonal Curvilinear
coordinate systems. J. Comp. Phys. 34:124-135.
Morel, F. 1983. Principles of Aquatic Chemistry. John Wiley & Sons, New York, NY. 446pp.
Morse, J.W., FJ. Millero, J.C. Cornwell, and D. Rickard, D. 1987. The chemistry of the hydrogen sulfide and iron
sulfide systems in natural waters. Earth-Science Reviews 24:1-42.
Moustafa, M.Z., and J.M. Hamrick. 1994. Modeling circulation and salinity transport in the Indian River Lagoon.
Estuarine and Coastal Modeling, Proceedings of the 3rd International Conference, ed. M. L. Spaulding et al.,
pp. 381-395. American Society of Civil Engineers, New York.
Mclntire, C.C. 1973. Periphyton dynamics in laboratory streams: a simulation model and its implications.
Ecological Monographs. 43:399-420.
Nixon, S. 1981. Remineralization and nutrient cycling in coastal marine ecosystems. Estuaries and Nutrients, ed.
Neilson and Cronin. pp.111-138. Humana Press, Clifton, NJ.
NOAA. 1998. Tide Tables 1998. National Oceanographic and Atmospheric Administration, National Ocean
Service, Silver Spring, MD.
O'Connor, D.J. and W .E. Dobbins. 1958. Mechanism of reaeration in natural streams. Transactions of the
American Society of Civil Engineers, 123(2934):641-684.
Odum, E.P. 1971. Fundamentals of ecology (3rd edition). W.B. Saunders Co., Philadelphia, PA. 574pp.
Oey, L.-Y., G.L. Mellor, and R.I. Hires. 1985. A three-dimensional simulation of the Hudson Raritan estuary. J.
Phys. Oceanogr. 15:1676-1720.
Park, K., A.Y. Kuo, J. Shen, and J.M. Hamrick. 1995. A three-dimensional hydrodynamic-eutrophication model
(HEM3D): description of water quality and sediment processes submodels, Special Report 327. The College of
William and Mary, Virginia Institute of Marine Science. 113 pp.
Parsons, T.R., M. Takahashi, and B. Hargrave. 1984. Biological Oceanographic processes (3rd edition). Pergamon
Press. 330 pp.
Pfeifer, R.F., and McDiffett, W. 1975. Some factors affecting primary production of stream communities. Archives
ofHydrobiol. 75:306-317.
Redfield, A.C., B.H. Ketchum, and F.A. Richards. 1963. The influence of organisms on the composition of sea-
water, ed. M.N. Hill, Chapter 2, pp 26-77. The Sea - Ideas and Observations on Progress in the Study of the Seas:
Vol. 2, Composition of Sea Water, Comparative and Descriptive Oceanography, Interscience Publishers.
Rennie, S., and J.M. Hamrick. 1992. Techniques for visualization of estuarine and coastal flow fields. Estuarine and
Coastal Modeling, Proceedings of the 2nd International Conference, ed. M. L. Spaulding etal., pp. 48-55.
American Society of Civil Engineers, New York.
Robbins, J.A., T. Keilty, D.S. White, and D.N. Edgington. 1989. Relationships among tubificid abundances,
sediment composition and accumulation rates in Lake Erie. Canadian J. of Fisheries & Aquatic Sciences,
46(2):223-231.
Rosati, A.K., and K. Miyakoda. 1988. A general circulation model for upper ocean simulation. J. Phys. Ocean,
18:1601-1626.
7-4
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Ross, P.J. 1983. Dynamics of periphyton communities. Periphyton in Freshwater Ecosystems, ed. R.G. Wetzel, pp.
5-10. Junk, Boston, MA, 5-10.
Runke, H.M. 1985. Simulation of the lotic periphyton community of a small mountain stream by digital computer.
PhD thesis, Utah State University, Logan, Utah.
Ryskin, G. and L.G. Leal. 1983. Orthogonal mapping. J. Comp. Phys. 50:71-100.
Sand-Jensen, K. 1983. Physical and chemical parameters regulating the growth of periphyton communities.
Periphyton in Freshwater Ecosystems, ed. R.G. Wetzel, pp. 63-71. Junk, Boston, MA.
Senior and Koerkle. 2003a. Simulation of streamflow and water quality in the Brandywine Creek subbasin of the
Christina River Basin, Pennsylvania and Delaware, 1994-98. U.S. Geological Survey Water-Resources
Investigations Report 02-4279, 207pp.
Senior and Koerkle. 2003b. Simulation of streamflow and water quality in the White Clay Creek subbasin of the
Christina River Basin, Pennsylvania and Delaware, 1994-98. U.S. Geological Survey Water-Resources
Investigations Report 03-4031, 242pp.
Senior and Koerkle. 2003c. Simulation of streamflow and water quality in the Red Clay Creek subbasin of the
Christina River Basin, Pennsylvania and Delaware, 1994-98. U.S. Geological Survey Water-Resources
Investigations Report 03-4138, 119pp.
Senior and Koerkle. 2003d. Simulation of streamflow and water quality in the Christina River subbasin and
overview of simulations in other subbasins of the Christina River Basin, Pennsylvania and Delaware, 1994-98. U.S.
Geological Survey Water-Resources Investigations Report 03-4193, 144pp.
Shen, J., J. Boon, and A.Y. Kuo. 1999. A numerical study of a tidal intrusion front and its impact on larval
dispersion in the James River estuary, Virginia. Estuary 22(3A):681-692.
Smolarkiewicz, P.K., and T.L. Clark. 1986. The multidimensional positive definite advection transport algorithm:
further development and applications. J. Comp. Phys. 67:396-438.
Smolarkiewicz, P.K., and W.W. Grabowski. 1990. The multidimensional positive definite advection transport
algorithm: nonoscillatory option. J. Comp. Phys. 86:355-375.
Smolarkiewicz, P.K., and L.G. Margolin. 1993. On forward-in-time differencing for fluids: extension to a
curvilinear framework. Man. Weather Rev. 121:1847-1859.
Steele, J.H. 1962. Environmental control of photosynthesis in the sea. Limnol. and Oceanogr. 7(2):137-150.
Stevenson, R.J., and Glover, R. 1993. Effects of algal density and current on ion transport through periphyton
communities. Limnol. and Oceanogr., 38(6): 1276-1281.
Stumm, W. and J.J. Morgan. 1981. Aquatic chemistry, an introduction emphasizing chemical equilibria in natural
waters (2nd edition). John Wiley & Sons, Inc. 780 pp.
Thomann, R.V., N.J. Jaworski, S.W. Nixon, H.W. Paerl, and J. Taft. 1985. The 1983 algal bloom in the Potomac
Estuary. The Algal Bloom Expert Panel, pepared for the Potomac Strategy State/EPA Management Committee, US
Environmental Protection Agency, Region III, Philadelphia, PA.
Thomann, R.V. and J.A. Mueller. 1987. Principles of surface water quality modeling and control. Harper & Row,
Publishers, Inc. 644 pp.
7-5
-------
Model Report for Christina River Basin, Nutrient and DO TMDL
Troup, R. 1974. The interaction of iron with phosphate, carbonate and sulfide in Chesapeake Bay interstitial waters:
a thermodynamic interpretation. Ph.D. Dissertation, Johns Hopkins University, MD. 114 p.
USEPA. 1990. Technical Guidance Manual for Performing Waste Load Allocations, Book 111 Estuaries, Part 2,
Application ofEstuarine Waste Load Allocation Models. EPA823-R-92-003. U.S. Environmental Protection
Agency, Office of Water. May 1990.
USEPA. 1995. Technical Guidance Manual for Developing Total Maximum Daily Loads, Book 11: Streams and
Rivers, Part 1: Biochemical Oxygen Demand/Dissolved Oxygen and Nutrients/Eutrophication. EPA 823-B-95-007.
U.S. Environmental Protection Agency, Office of Water. September 1995.
USEPA. 2000. Hydrodynamic and water quality model of Christina River Basin. U.S. Environmental Protection
Agency, Region III, Philadelphia, PA. December 2000.
USEPA. 2004. Christina River Basin High-Flow TMDLs for Nutrients, Low Dissolved Oxygen, and Bacteria: Data
Report (Draft). U.S. Environmental Protection Agency, Region III, Philadelphia, PA. April 30, 2004.
USEPA. 2005. Total Maximum Daily Loads for Nutrients and Low Dissolved Oxygen in the Christina River Basin,
Pennsylvania, Delaware, and Maryland. USEPA Region III, Philadelphia, PA. April 8, 2005.
Warwick, J.J., D. Cockrum, andM. Horvath. 1997. Estimating non-point source loads and associated water quality
impacts. ASCE J. Water Res. Plan, and Manage. 123(5):302-310.
Westrich, J.T. and B.A. Berner. 1984. The role of sedimentary organic matter in bacterial sulfate reduction: the G
model tested. Limnol. andOceanogr. 29(2):236-249.
Wezernak, C.T. and J.J. Gannon. 1968. Evaluation of nitrification in streams. ASCE J. Sanitary Engr. Div.
94(5):883-895.
Whitford, L.A. and G.J. Schumacher. 1964. Effect of current on respiration and mineral uptake of Spirogyra and
Oedogonium. Ecology 45:168-170.
Wu, T.S., J.M. Hamrick, S.C. McCutechon, and R.B. Ambrose. 1996. Benchmarking the EFDC/HEM3D surface
water hydrodynamic and eutrophication models. Computational Methods for Next Generation Environmental
Models, ed. G. Delich. Society of Industrial and Applied Mathematics, Philadelphia.
Yamamoto, S., J.B. Alcauskas, and T.E. Crozler. 1976. Solubility of methane in distilled water and seawater. J. of
Chemical and Engineering Data 21(1):78-80.
7-6
------- |