&EPA
United States
Environmental Protection
Agency
EPA/600/R15-152 June 2015
www.epa.gov/ord
THREE-DIMENSIONAL
MODELING
OF
HYDRODYNAMICS
AND
TRANSPORT
IN
NARRAGANSETT BAY
Office of Research and Development
National Health and Environmental Effects Research Laboratory
-------
-------
vvEPA
United States EPA/600/R-15/152
Environmental Protictior
Agency www.epa.gov/ord
Three-Dimensional Modeling of
Hydrodynamics and Transport
in Narragansett Bay
Mohamed A. Abdelrhman
US Environmental Protection Agency
Office of Research and Development
National Health and Environmental Effects Research Laboratory
Atlantic Ecology Division,
27 Tarzwell Drive,
Narragansett, RI 02882 USA
-------
DISCLAIMER
This document is a final draft. It has not been formally released by the U.S. Environmental
Protection Agency and should not at this stage be construed to represent Agency policy. It is
being circulated for comments on its technical merit and policy implications. Although the
material described here has been funded by the U.S. Environmental Protection Agency, it has not
been subject to Agency-level review and therefore does not necessarily reflect the views of the
Agency, nor does mentioning trade names or commercial products endorse or recommend them.
This report has the ORD Tracking Number ORD-008162 of USEPA Office of Research and
Development, National Health and Environmental Effects Research Laboratory, Atlantic
Ecology Division.
-------
CONTENTS
DISCLAIMER ii
CONTENTS iii
FIGURES v
TABLES x
ABBREVIATIONS AND ACRONYMS xi
ACKNOWLEDGMENTS xi
EXECUTIVE SUMMARY xii
1. Introduction 1
1.1. Background 1
1.2. Objectives 2
1.3. Approach 2
1.4. Physical Setting 3
1.5. Quality Assurance and Quality Control 4
1.6. Data 6
2. The Hydrodynamic and Transport Model 31
2.1. Model Equations 31
2.2. Model Configuration 35
2.2.1. Numerical grid 36
2.2.2. Boundary forcing functions 36
2.2.3. Initial conditions for temperature and salinity 37
2.2.4. Time step and duration 38
3. Model Calibration, Validation, and Skill 47
3.1. Water Surface Elevation 47
3.2. Temperature 48
3.3. Salinity 49
3.4. Velocity 49
3.5. Model Skill Analysis 50
3.6. The Overall Behavior of Narragansett Bay 50
4. Sample Results 97
4.1. Water Surface Elevation Distribution 97
4.2. Velocity Distribution 97
4.3. Temperature Distribution 97
4.4. Salinity Distribution 97
4.5. Water Density and Stratification Strength 98
4.6. Dye Distribution 98
4.6.1. Dye building-up experiment 98
4.6.2. Dye flushing experiment 98
4.6.2.1. Flushing of estuarine water 99
4.6.2.2. Flushing of freshwater 99
4.6.2.3. Flushing of saltwater 100
5. Summary and Conclusion 101
5.1. More on Calibration 101
5.2. Hydrodynamics and Water Quality 102
5.2.1. Water quality with WASP7 102
5.2.2. Water quality with EFDC 102
5.3. Future Scenarios 103
-------
Appendices
Appendix A. Definition of Skill Parameters 107
Appendix B. Input files 109
Appendix C. Elevation Contours 110
Appendix D. Horizontal Velocity Vectors 112
Appendix E. Temperature Contours 116
Appendix F. Salinity Contours 119
Appendix G. Temperature Stratification 122
Appendix H. Salinity Stratification 129
Appendix I. Vertical Density Gradient 136
Appendix J. Profiles of Temperature, Salinity, Density, and Stratification Strength 142
Appendix K. Dye Flushing 158
Appendix L. Flushing Profiles 166
Appendix M. Dye Build up 168
REFERENCES 173
IV
-------
FIGURES
Figure 1. General layout of Narragansett Bay with its major islands, passages, and rivers
with locations of data stations and WWTPs 17
Figure 2. Narragansett Bay watershed 18
Figure 3. Bathymetry (feet) of Narragansett Bay 20
Figure 4. Navigation charts for Narragansett Bay showing all station locations 21
Figure 5. Digital bathymetric data raster with depths in meters (left) and sample of actual
sounding locations (right) 22
Figure 6. Freshwater inflow rates from all rivers and sub-watersheds 23
Figure 7. Freshwater effluent from WWTPs discharging directly into Narragansett Bay 24
Figure 8. River freshwater temperature and air temperature resembling a sine wave 24
Figure 9. Precipitation rate at station PR 25
Figure 10. Wind data at station PR: (A) Wind speed and direction, (B) Wind power (U2),
and (C) Wind rose 25
Figure 11. Incident shortwave solar radiation on the earth's surface at the airport station TFG.. 27
Figure 12. Cloud cover at the airport station TFG 27
Figure 13. Relative humidity at station PR 28
Figure 14. Atmospheric pressure at station PR 28
Figure 15. Water temperature at station NP 29
Figure 16. Salinity at station GD 29
Figure 17. Water surface elevation at station NP 30
Figure 18. Digital representation of numerical grid 40
Figure 19. Numerical grid and station locations for Narragansett Bay: (A) Data Stations,
(B) WWTPs, (C) Riverine and Riparian flow 41
Figure 20. Bottom elevation for the numerical grid 44
-------
Figure 21. Quasi-steady state distribution of dye buildup to approximate initial distributions.... 45
Figure 22. Time series of observed and predicted water surface elevation at station PR 55
Figure 23. Scatter plot between predicted and observed elevations at station PR 55
Figure 24. Time series of predicted and observed sub-tidal water surface variations
at station PR 56
Figure 25. Scatter plot between predicted and observed sub-tidal elevations at station PR 56
Figure 26. Time series of predicted and observed water surface elevation at station FR1 57
Figure 27. Scatter plot of predicted and observed water surface elevation at station FR1 57
Figure 28. Time series of sub-tidal variations at station FR1 58
Figure 29. Scatter plot between predicted and observed sub-tidal elevations at station FR1 58
Figure 30. Time series of predicted and observed surface water temperature at station PR 59
Figure 31. Scatter plot of predicted and observed surface water temperature at station PR 59
Figure 32. Time series of predicted and observed surface temperature at station FR1 60
Figure 33. Scatter plot of predicted and observed surface water temperature at station FR1 60
Figure 34. Time series of predicted and observed surface water temperature at stationNP 61
Figure 35. Scatter plot of predicted and observed surface water temperature at station NP 61
Figure 36. Time series of predicted and observed surface (A) and bottom (B)
water temperature at stationPD 62
Figure 37. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station BR 63
Figure 38. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station CP 64
Figure 39. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station GB 65
Figure 40. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station SR 66
vi
-------
Figure 41. Time series of predicted and observed surface (A) and bottom (B)
water temperature at stationNPI 67
Figure 42. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station MV 68
Figure 43. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station QP 69
Figure 44. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station TW 70
Figure 45. Time series of predicted and observed surface (A) and bottom (B)
water temperature at station PP 71
Figure 46. Time series of predicted and observed surface (A) and bottom (B)
water temperature at stationMH 72
Figure 47. Time series of predicted and observed surface water temperature at station GD 73
Figure 48. Time series of predicted and observed surface salinity at station PR 73
Figure 49. Time series of predicted and observed surface (A) and bottom (B)
salinity at stationPD 74
Figure 50. Time series of predicted and observed surface (A) and bottom (B)
salinity at stationBR 75
Figure 51. Time series of predicted and observed surface (A) and bottom (B)
salinity at station CP 76
Figure 52. Time series of predicted and observed surface (A) and bottom (B)
salinity at station GB 77
Figure 53. Time series of predicted and observed surface (A) and bottom (B)
salinity at station SR 78
Figure 54. Time series of predicted and observed surface (A) and bottom (B)
salinity at stationNPI 79
Figure 55. Time series of predicted and observed surface (A) and bottom (B)
salinity at stationMV 80
Figure 56. Time series of predicted and observed surface (A) and bottom (B)
salinity at station QP 81
vii
-------
Figure 57. Time series of predicted and observed surface (A) and bottom (B)
salinity at station TW 82
Figure 58. Time series of predicted and observed surface (A) and bottom (B)
salinity at station PP 83
Figure 59. Time series of predicted and observed surface (A) and bottom (B)
salinity at station MH 84
Figure 60. Time series of predicted and observed surface salinity at station GD 85
Figure 61. Time series of observed and predicted current speed at station QP2
in February 2009 86
Figure 62. Time series of observed and predicted current direction at station QP2
in February 2009 86
Figure 63. Time series of observed and predicted eastward (u) velocity at Station QP2
in February 2009 87
Figure 64. Time series of observed and predicted northward (v) velocity at station QP2
in February 2009 87
Figure 65. Time series of observed and predicted current speed at station FR2
in February 2009 88
Figure 66. Time series of observed and predicted current direction at station FR2
in February 2009 88
Figure 67. Time series of observed and predicted eastward (u) velocity at station FR2
in February 2009 89
Figure 68. Time series of observed and predicted northward (v) velocity at station FR2
in February 2009 89
Figure 69. Time series of predicted and observed surface water temperatures for
the whole Bay 90
Figure 70. Scatter plot of predicted and observed surface water temperatures for
the whole Bay 90
Figure 71. Time series of predicted and observed bottom water temperatures for
the whole Bay 91
VIII
-------
Figure 72. Scatter plot of predicted and observed bottom water temperatures for
the whole Bay 91
Figure 73. Time series of observed surface and bottom water temperatures to illustrate
temperature stratification for the whole Bay 92
Figure 74. Time series of predicted surface and bottom water temperatures to illustrate
temperature stratification for the whole Bay 92
Figure 75. Time series of predicted and observed surface water salinities for the whole Bay 93
Figure 76. Scatter plot of predicted and observed surface water salinities for the whole Bay 93
Figure 77. Time series of predicted and observed bottom water salinities for the whole Bay 94
Figure 78. Scatter plot of predicted and observed bottom water salinities for the whole Bay 94
Figure 79. Time series of observed surface and bottom salinities to illustrate salinity
stratification for the whole Bay 95
Figure 80. Time series of predicted surface and bottom salinities to illustrate salinity
stratification for the whole Bay 95
Figure 81. Observed water temperatures at PR, QP1, CP, and NP during 2003 105
IX
-------
TABLES
Table 1. Sub-watersheds and main rivers in Narragansett Bay watershed 11
Table 2. Sub-watersheds of Greenwich Bay 12
Table 3. Eleven WWTPs with direct discharge into Narragansett Bay 13
Table 4. Stations and locations used for historical data 14
Table 5. Tidal stations used to force, calibrate, and validate hydrodynamics
in Narragansett Bay 16
Table 6. Open boundary node string 39
Table 7. Skill parameters for calibration and validation 52
TableS. Local flushing time at stations PR, CP,GB,TW, and FR1 100
Table 9. Tidal constituents at Newport 105
-------
ABBREVIATIONS AND ACRONYMS
AED Atlantic Ecology Division
EFDC Environmental Fluid Dynamics Code
DMR Discharge Monitoring Report
GEODAS National Geophysical Data Center
GIS Global Information System
LFT Local Flushing Time
MLLW Mean Lowest Low Water
MSL Mean Sea Level
NAD83 North American Datum of 1983
NB Narragansett Bay
NCDC National Climate Data Center
NOAA National Ocean and Atmospheric Administration
NOS National Ocean Service
NSRDB National Solar Radiation Data Base
QA/QC Quality Assurance/Quality Control
RENL Renewable Energy National Laboratory
RIDEM Rhode Island Department of Environmental Management
RIGIS Rhode Island Geographic Information System
ROMS Regional Ocean Modeling System
SMS Surface-water Modeling System
USEPA United States Environmental Protection Agency
USGS United States Geological Survey
UTM Universal Transverse Mercator
WASP Water Quality Analysis Simulation Program
WWTP Wastewater Treatment Plant
ACKNOWLEDGMENTS
The author thanks Mr. Brandon Jarvis (USEPA-GED) for providing most of the preprocessing
software to prepare input files for EFDC. Thanks are also due to Drs. John Hamrick (Tetra Tech.,
Inc.) and Tim Wool (USEPA-Region 4) for providing the executables for EFDC and WASP7
models. The author appreciates the GIS mapping assistance provided by Jane Copeland
(USEPA-AED). The author is grateful to Dr. Earl J. Hayter (US Army Corps of Engineers,
Engineer Research and Development Center) and Dr. Hugo Rodrigues (Tetra Tech, Inc.) for their
support with EFDC technical issues. The author wishes to thank the reviewers of this report,
including Drs. Edward Dettmann, Henry Walker, Brenda Rashleigh, and Glen Thursby (USEPA-
AED) as well as Mr. Brandon Jarvis (USEPA-GED) for their technical reviews, insights, and
constructive comments.
XI
-------
EXECUTIVE SUMMARY
The work presented here addresses the specific needs of physical information required by tasks
within the two SSWR projects:
• SSWR 2.3. A (Nutrient management for sustainability of upland and coastal ecosystems:
Building a locally applicable management tool box for application across the US); and
• SSWR 6.1 (Narragansett Bay and Watershed Sustainability - Demonstration Project).
The main objective of this work is to develop the methodology to generate the hydrodynamics
needed to predict the transport of constituents and properties that govern water quality and
ecology in estuarine and coastal systems. The specific objective is to apply the developed
methods to identify hydrodynamics and transport mechanisms in the Narragansett Bay, RI, to
serve as a prototype for future implementation on other systems.
This report presents the methodology to apply, calibrate, and validate a three-dimensional
hydrodynamic model to provide the advection and dispersion mechanisms required by water
quality and ecological models. The methodology is applied to Narragansett Bay, RI to generate
sample predictions of hydrodynamics and transport for the year 2009. The same methodology
can be applied to other years and periods of interest. The EPA-recommended Environmental
Fluid Dynamics Code (EFDC) is chosen for this work.
The approach presented here attempts to structure the model horizontal and vertical grid
resolution to resolve the hydrodynamic, water quality, and ecological needs using computations
with the same spatial and temporal resolution. The temporally resolved sub-daily variations are
intended to adequately address water quality needs (e.g., for temperature, light, oxygen, and
phytoplankton growth). The process of calibration and validation of model results did not include
any newly implemented field monitoring programs; existing historical data were utilized. The
generated hydrodynamics would be used with the water quality modules in the EPA's Water
Quality Analysis Simulation Program (WASP) to provide predictions for present and future
scenarios of interest. The quality of the hydrodynamic predictions provided to WASP was
maintained by a thorough calibration, validation, and performance evaluation of the
hydrodynamic model. These data were internally calculated by EFDC and they included volume
of segment, depth of segment, velocity through segment, flow along each flow path between
segments, and dispersion along each flow path between segments.
The report includes sections with detailed information, tables, and figures related to the physical
setting and available data, model equations, and model configuration for Narragansett Bay,
model calibration and validation, a summary of the skill parameters used to evaluate model
performance, and a sample of model results. The summary and conclusion section includes a
brief discussion of the use of the hydrodynamic model results with the WASP water quality
model and some suggestions to conduct hydrodynamic simulations of future scenarios.
XII
-------
1. INTRODUCTION
Practical and realistic modeling of hydrodynamics is the first step towards proper modeling of
ecological phenomena in estuaries because the hydrodynamics control the movement and
dispersion of constituents important to ecology. The practicality of a hydrodynamic model lies in
the reasonable spatial and temporal resolutions which dictate the ease and timely production of
flow as well as the transport and dispersion of constituents. These results can be easily checked
against observed data. A key factor in this process is the choice of the spatial resolution for the
hydrodynamic model.
This report presents the methodology to apply, calibrate, and validate a three-dimensional (3-D)
hydrodynamic model to provide the advection and dispersion mechanisms required by water
quality and ecological models. The methodology is applied to Narragansett Bay (NB, Bay, or
system) in Rhode Island (RI) to generate sample predictions of hydrodynamics and transport for
the year 2009.
The organization of this report includes five sections with relevant tables and figures at the end
of each. Section 1 (this section) presents the background, objectives, approach, physical setting,
and available data. Section 2 presents information about model equations and model
configuration for NB (numerical grid, forcing functions, initial conditions, time step, and input
files). Section 3 covers model calibration and validation for water surface elevation, temperature,
salinity, and velocity, as well as a summary of skill parameters used to assess model
performance. Section 4 presents a sample of the predicted horizontal distributions of water
surface elevation, velocity vectors, temperature, salinity, and dye, in addition to analysis of
flushing times of freshwater, seawater, and estuarine water within the Bay. Section 5 presents the
summary and conclusion with brief discussion of the use of the hydrodynamic model results with
the water quality model and some suggestions to conduct hydrodynamic simulations of future
scenarios to examine effects of changes in land use and pollution control activities as well as
anticipated effects of climate change (e.g., changes in precipitation patterns and amounts as well
as sea level rise). Definitions of the skill parameters are presented in Appendix A. Input files
required to run the model are listed in Appendix B. More results are presented in appendices
C-M.
1.1. Background
The Environmental Protection Agency's (EPA's) Safe and Sustainable Water Resources (SSWR)
program seeks to ensure that clean, adequate, and equitable supplies of water are available to
support the well-being of humans and aquatic ecosystems. The EPA requires the development of
numeric nutrient criteria (i.e., allowable concentrations of nutrients) for estuarine and coastal
waters to reduce undesirable impacts to beneficial uses. Changes in nutrient concentrations have
to be predicted in both space and time. Critical relationships between nutrient concentrations and
biological responses have to be evaluated through the application of hydrodynamic, water
quality, and ecological models.
-------
The work presented here addresses the specific needs of hydrodynamic information required by
tasks within the two SSWR projects: SSWR 2.3.A (Nutrient management for sustainability of
upland and coastal ecosystems: Building a locally applicable management tool box for
application across the US); and SSWR 6.1 (Narragansett Bay and Watershed Sustainability -
Demonstration Project).
1.2. Objectives
The main objective of this work is to develop the methodology to generate the hydrodynamic
information needed to predict transport of constituents and properties that govern water quality
and ecology in estuarine and coastal systems. The methodology should include the process used
to communicate hydrodynamic information to water quality and ecological models. The specific
objective is to apply the developed methods to identify hydrodynamic and transport mechanisms
to guide decision-making in the NB watershed and as a prototype for future implementation in
other systems.
The inherent hydrodynamic complexities in coastal and estuarine systems arise from salinity and
temperature stratification, estuarine circulation, meteorological (e.g., wind) impact, and
astronomical (e.g., tide) forcing. These complexities require a 3-D model to provide adequate
predictions offerees on and movements of water masses (with their associated constituents) in
both space and time. The EPA-recommended Environmental Fluid Dynamics Code (EFDC) is
chosen for this work.
1.3. Approach
Multidimensional hydrodynamic models usually use fine grid (on the order of meters or tens of
meters) to resolve the effect of geometrical changes in bathymetry and shoreline on the flow and
circulation. Model stability requires a very short time step (on the order of seconds) to
accommodate such fine spatial resolution. For example, Zhao et al. (2006) applied the Finite
Volume Coastal Ocean Model (FVCOM) (Chen et al. 2006) to NB and Mount Hope Bay. The
grid resolution was 50 m and the time step was on the order of one second. The same model was
applied to Breton Sound estuary, Mississippi (Huang et al. 2011) using a spatial resolution and
time step of 20 m and 0.5 s, respectively. The Regional Ocean Modeling System (ROMS)
(Haidvogel et al., 2000; Shchepetkin and McWilliams, 2005) was also applied to NB with a
spatial resolution of 50 m and a time step of 8 s (personal communication, David Ullman,
Graduate School of Oceanography, University of Rhode Island).
The above-mentioned models can simulate hydrodynamics and water quality at the same spatial
and temporal resolutions used by the hydrodynamic model. However, this approach becomes
prohibitive when large regions and long simulations are needed. To overcome this difficulty,
Kremer et al. (2010) used a computational scheme that implemented ROMS to simulate the fine-
scale hydrodynamics. They aggregated the dense spatial and temporal information into a set of
matrices to identify daily exchanges between the coarser water quality segments throughout a
full year. This approach required extensive hydrodynamic runs to cover all the possible
-------
exchanges between the water quality segments (with their nested hydrodynamic cells). Such runs
have to be repeated with every change in the hydrodynamic setting, e.g., to accommodate future
scenarios of global warming, climate change, sea-level rise, and river inflow.
The approach presented here attempts to structure the model horizontal and vertical grid
resolution to resolve the hydrodynamic, water quality, and ecological needs without any spatial
or temporal aggregation. Thus, the spatial resolution should adequately resolve the
hydrodynamics; meanwhile, it should be reasonably coarse to accommodate the dimensionality
limitations for water quality and ecological calculations (See Section 2.2.1 - Numerical Grid). In
addition to insuring the stability of the hydrodynamic model, the temporal resolution should
resolve sub-daily variations to adequately address water quality assessment needs (e.g., for
temperature, light, oxygen, and phytoplankton growth). The water quality modules in the EPA's
Water Quality Analysis Simulation Program (WASP) can then be used to provide predictions for
present and future scenarios of interest. Wool et al. (2002) applied the same approach with
EFDC and WASP to study hydrodynamics and water quality in the Neuse River Estuary, NC.
The following sections include a description of the hydrodynamic and transport model, EFDC,
and its application to the NB system. The process of calibration and validation of model results
did not include any field monitoring programs; existing literature and field based datasets were
utilized.
1.4. Physical Setting
The NB has three major islands (Conanicut, Prudence, and Aquidneck), which confine the flow
into the Sakonnet River, the East Passage, and the West Passage (Fig. 1). In addition, there are
many small islands that were ignored since they have insignificant (localized) effects on
circulation and mixing, as noticed in a separate modeling effort using the FVCOM with very fine
segment size (~ 25 m, not presented).
The Bay has four boundaries: the landward boundary with its adjacent watershed, the seaward
(open) boundary with RI Sound, the surface boundary with the overlying atmosphere, and the
bottom boundary with the underlying bed sediment. The landward boundary is defined by the
shoreline of NB. The watershed of NB has an area of 4,353.49 km2 (RIGIS, personal
communication: Jane Copeland, USEPA-Atlantic Ecology Division, Narragansett, RI). The
watershed is located within the states of RI and Massachusetts (MA). It is composed of nine sub-
watersheds (Table 1 and Fig. 2 A) with eight of which draining surface and ground water through
gauged rivers. The ninth sub-watershed (498.03 km2) lies along the Bay's shore line and it drains
directly into the Bay. This sub-watershed is referred to as "riparian" in this report. Nutrients and
contaminants from point sources such as inland wastewater treatment plants (WWTPs) and from
nonpoint sources enter the Bay with the freshwater inflow from these rivers. In addition, eleven
WWTPs are located along the shore line and discharge their loads directly into the Bay.
The Greenwich Bay is a sub-embayment of NB with seven small sub-watersheds having a total
area of 52 km2 (Fig 2B) (RIDEM, 2005). These sub-watersheds were gathered into three major
-------
groups which discharge their freshwater flows at the grid locations shown in Table (2). The
freshwater flow for each group was prorated from the Hunt River discharge multiplied by the
ratio between the watershed area of the group and that of the Hunt River (64.09 km2, Table 1).
The calculated flow was increased by 5% to compensate for ground water inflow to Greenwich
Bay (see Section 1.6 C). The Greenwich Bay watershed was subtracted from the riparian area
assigned to the Hunt River (62.25 Km2, Table 2) to reduce the net riparian area assignment to
10.25 km2.
The seaward open boundary of the Bay is forced by the semidiurnal tides in the RI Sound. Also,
water temperature and salinity in the RI Sound control and moderate their respective values
inside the Bay. Similarly, concentrations of other water quality parameters (e.g., nutrients) in the
Bay are affected by their values at the seaward open boundary.
Meteorological forcing and atmospheric deposition impact the Bay through its surface area,
which is 385.91 km2 (RIGIS, personal communication: Jane Copeland, USEPA-Atlantic Ecology
Division, Narragansett, RI). The major meteorological forces include precipitation, evaporation,
solar radiation, wind speed and direction, air temperature, and atmospheric pressure.
The bottom boundary is defined by NB bathymetry (Fig. 3). The mean depth of the Bay is 8.0 m
(including Sakonnet River, Pilson 1985) with a maximum depth of 46 m in the east passage
(Zhao et al. 2006). In addition to friction and drag exerted by the bed on the overlying flow, the
bed acts as a source/sink to many of the dissolved and suspended materials in the overlying water
column.
1.5. Quality Assurance and Quality Control
The quality assurance and quality control (QA/QC) of the hydrodynamic and transport modeling
in NB was driven by the QA/QC of the hydrodynamic model (EFDC), the used historical data,
and the model predictions.
i. The EFDC model is supported by the USEP A, which states the following on its
website: (http://www2.epa.gov/exposure-assessment-models/efdc-read-me-epa-
version-101), "The Environmental Fluid Dynamics Code (EFDC) is a multifunctional
surface water modeling system, which includes hydrodynamic, sediment-contaminant,
and eutrophication components. The public domain EFDC model was developed by
Dr. John M. Hamrick at the Virginia Institute of Marine Science and is currently
maintained by Tetra Tech, Inc. with support from the U.S. EPA. EFDC has been used
for more than 80 modeling studies of rivers, lakes, estuaries, coastal regions and
wetlands in the United States and abroad." Thus, no further modeling QA/QC effort
for the model itself was pursued in this work. However, a recommendation for
separating horizontal diffusion of mass and momentum was made to improve
calibration and validation of model results (Section 5.1). In addition, it was noticed
that the FORTRAN code by-passed the main subroutine which calculates horizontal
-------
diffusion of mass and property. This deficiency was fixed and the separation of
momentum and mass diffusivities was introduced to the code at the USEPA-AED.
However, the work presented in this report was based on the original unmodified
code.
ii. Available historical data were used to force, calibrate, and validate the predicted
hydrodynamics in NB. Data sources included federal, state, and academic institutions
including: The National Oceanic and Atmospheric Administration (NOAA), the
United States Geological Surveys (USGS), the USEPA Data Monitoring Report
(DMR), the National Renewable Energy Laboratory (NREL), and the Rhode Island
Department of Environmental Management (RIDEM) with the Narragansett Bay
Commission (NBC) and the University of Rhode Island-Graduate School of
Oceanography (URI-GSO). These institutions have their own QA/QC requirements,
which were not available. However, the following information was provided on the
institution's website or with the posted data:
• The NOAA data included: water data (e.g., surface elevation, temperature,
conductivity, velocity), meteorological data (e.g., air temperature, pressure, wind
speed and direction), and bathymetric maps. The NOAA's Center for Operational
Oceanographic Products and Services (CO-OPS) environmental measurement
systems sensor specifications and measurement algorithm resided in
http://tidesandcurrents.noaa.gov/publications/CO-
OPS_Measurement_SpecUpdated_4.pdf Due to the low quality of NOAA
conductivity/salinity data (presented in the first draft of this report) the data were
excluded (except for data at Newport) in this final report and they were replaced
with data from the twelve buoys mentioned below.
• The RIDEM/NBC/URI-GSO buoy data included water temperature and salinity.
These data have been subject to corrections after post-deployment and re-
calibration of sensors. Notes provided with data from all the NBC's stations
indicated that data spikes that exceeded the 2 standard-deviation rule and/or the 2
X the nearest neighbor rule were deleted as well as data which demonstrated
fouling effect or calibration offset. In addition, a disclaimer indicated that some
data are preliminary and subject to changes based on QA/QC procedures and
contacts were provided for additional information on the use of the data.
• The USGS data included riverine discharge rates. Quality guidelines were posted
on http://www.usgs.gov/info qual. Under typical scenarios, uncertainty for stream
flow measurements ranges from 6% to 19% (Harmel et al. 2006)
http://pubs.er.usgs.gov/publication/70030364
• The USEPA-DMR data included effluent from WWTPs. Contacts for QA of these
data were provided at http://www.eraqc.com/DMRQA/USEPA.
-------
• The NREL data included solar radiation and meteorological data. The incident
shortwave solar radiation had uncertainty of 8-25% (or more) for hourly values
(Wilcox, 2012). In this work, data for 2009 were reduced by 20% to reproduce
observed water temperatures at 12 buoy locations within the Bay.
All data were visually inspected and evaluated to identify data gaps as well as
erroneous and superfluous values. Gaps in time series of forcing functions were
interpolated between earlier and later values, while calibration and validation data
sets were left un-touched, as displayed in many figures. Dates of all missing data
were preserved in the time series. Erroneous and anomalous values were identified
visually and the full duration when they appeared was excluded. More investigation is
required to identify the cause of the predicted higher water temperatures and if it is
caused by the high uncertainty of the NREL's solar radiation data or a deficient
formulation of heat losses in EFDC. The following section provides the sources of all
the data that were used in this work. Brief statements about the quality of the data are
provided when needed.
iii. Quality of the predicted data was evaluated based on eight skill parameters for model
performance (Appendix A). The eight performance skill parameters were applied
with predictions at all stations (Table 6). In this work, the acceptance level for the
skill parameters was arbitrarily set at or above 80%. Many skill parameters exceeded
the 90% acceptance level (see red numbers in Table 6).
1.6. Data
Available information and historical data for the year 2009 were used to force, calibrate,
validate, and predict the hydrodynamics in NB.
A. System boundary. Three navigation charts were used to fully cover the extent of the NB
water surface: NOAA-NOS navigation charts No. 13221 (main area of the Bay), 13224
(Seekonk River), and 13226 (Taunton River) (e.g.,
http://www.charts.noaa.gov/OnLineViewer/13221.shtml) (Fig. 4). These maps were geo-
referenced in the Universal Transverse Mercator (UTM) projected coordinate system,
North American Datum of 1983 (NAD83), and grid Zone 19 (between longitudes 66 °W
to 72° W, http://en.wikipedia.Org/wiki/File:Utm-zones.svg). The maps were used to lay
the numerical grid and identify the landward boundaries as well as the southern seaward
open boundary with RI Sound. Locations of the eight rivers (Table 1) were identified on
the maps together with the locations of the stations used to obtain historical data for
calibration and validation (Table 3, Fig. 4).
B. Bottom bathymetry. Highly detailed bathymetry was obtained from the observed depth
data recorded on the National Geophysical Data Center (GEODAS) bathymetry CDs
(http://www.ngdc.noaa.gov/mgg/geodas/geodas.html). These data were referenced
-------
horizontally in UTM, NAD83, Zone 19. Bathymetric data were manipulated using GIS.
The mean depth option of the Point to Raster tool was used to assign the depth and the
horizontal coordinate location to each raster cell (Fig. 5). The raster was converted to
point shape file and a data table was generated. Depth data were referenced to the mean
lowest low water (MLLW) (Table 4). These data were converted to depth from the mean
sea level (MSL). Table 4 presents tidal information about the three NOAA stations NP,
PR, and FR (Table 3). Table 4 indicates that both MLLW and MSL are not the same at
NP and PR (MSLPR-MSLNp = 0.025 m, and MLLWNp-MLLWPR= 0.132 m). All depth
values were referenced to MSL at NP (MSLNP = 0.529 m) by adding a correction A =
(0.132+2(0.529))/2= 0.6 m to all reported bathymetric values. The grid generation routine
identified water depth at each model grid node as described below.
Four dredged navigation channels exist in the Bay with design depths from MLLW as
follows: Providence River channel 12.2 m (40 ft), Seekonk River channel 4.88 m (16 ft),
Mount Hope Bay-Fall River Harbor Channel 10.67 m (35 ft), and Quonset Point Channel
10.67 m (35 ft). An adjustment of 0.6 m is added to channel depths to reference them to
MSL at Newport for subsequent use at the respective grid locations.
C. Freshwater inflow. Freshwater can enter the estuary from three main sources: point
discharge from rivers and WWTPs, direct precipitation on the water surface, and the
freshwater load from the un-gauged area adjacent to the Bay (called here "riparian").
Groundwater is assumed to be implicitly included in the river and riparian flows. Rough
estimate of groundwater input is 5% of the total freshwater input to the Bay (Nowicki and
Gold 2008). Most groundwater comes from areas 1-2 km along the shoreline and
groundwater flow from the bottom is unavailable (Nowicki and Gold 2008). Evaporation
from the water surface is the only loss of freshwater. Using meteorological parameters
(see below), the hydrodynamic model internally calculates evaporation losses.
The riverine freshwater is delivered to NB through eight rivers (Table 1):
Woonasquatucket, Pawtuxet, Ten Mile, Taunton, Moshassuck, Blackstone, Hunt,
Warren/Palmer. Time series of freshwater inflow rates in cubic feet per second at daily
intervals were obtained for each river from the US Geological Surveys (USGS)
(e.g., http://waterdata.usgs.gov/RI/nwis/current?type=flow). Daily freshwater inflows
in m3 s"1 were calculated for 2009 and assigned to respective grid locations.
Flow in the Taunton River included flows from Mill and 3-mile Rivers. The
Warren/Palmer River did not have a USGS station and its flow was assumed proportional
to the flow in the adjacent Ten Mile River prorated by the ratio between the two drainage
areas (175.33 km2 / 144.29 km2=1.2749). Flow from each river was prorated (increased)
by the ratio between its watershed and its gauged drainage area (Table 1). Flow from the
-------
riparian sub-watershed (498.03 km2) was calculated from the average flow per unit area
from the total gauged area (3,108.53 km2). Based on visual inspection of the riparian area
(Fig. 2), half the riparian flow was distributed equally to the Pawtuxet, Taunton, Hunt,
and Warren/Palmer Rivers, and the other half was assigned to three locations centrally
located within the East Branch, West Branch, and Sakonnet River. The overall proration
factor for each river, which accounts for ungauged watershed and ground water inflow, is
presented in Table 1. Figure 6 shows freshwater inflow rates from all rivers. Salinity of
freshwater inflow was set to zero and its temperature was assumed to be the same as air
temperature (see below).
Effluents from the eleven WWTPs were obtained from historical records (Table 2 and
Fig. 7). Effluents from Fields Point and Bucklin Point were recorded three times per
week and they were obtained from the Narragansett Bay Commission (personal
communication, Edward Dettmann, USEPA-AED, Narragansett, RI). Effluents from the
other nine plants were downloaded from the US EPA's Integrated Compliance
Information System (ICIS) using the EPA's Discharge Monitoring Report (DMR)
pollutant loading tool (http://cfpub.epa.gov/dmr/). The additional freshwater supply to
WWTPs came from reservoirs and/or groundwater, not from rivers. The volume of the
flow through the Brayton Point power plant was not considered here because it was
withdrawn from and released to Mount Hope Bay. In addition, the power plant is
scheduled for closure in 2017 which will eliminate its impact on future scenarios and
decisions.
D. Meteorological and atmospheric parameters These parameters include: atmospheric
pressure (mb), dry air temperature (°C), relative humidity (%), precipitation rate (in h"1),
evaporation rate (calculated internally), hourly short-wave irradiation (Wh m"2),
cloud cover (%), wind speed (m s"1), wind direction of blowing-to (deg.) from North.
Meteorological data were obtained from the National Oceanic and Atmospheric
Administration (NOAA)-National Climate Data Center (NCDC) and National Solar
Radiation Data Base (NSRDB) at T.F. Green Airport (Station 725070 for solar radiation
and cloud cover, Station 8454000 for wind and air pressure and other meteorological
data) at hourly increments during the year 2009 and assumed to be uniformly distributed
over the entire surface area of NB. The NCDC-NSRBD tabular data included 49 columns
with Short-Wave irradiation (wh m"2) in column 7, total Sky Cover (%) in column 22,
Relative Humidity (%) in column 30, precipitation (mm) in column 42 and precipitation
duration (h) in column 44 (Wilcox, 2012). The model used hourly meteorological data for
wind speed (m s"1) and direction (degree from North), air temperature (°C) and pressure
(mb) from NOAA Station 8454000 (Providence, RI)
(http ://tidesandcurrents.noaa. gov/met.html?bdate=20090101 &edate=20091231 &units=m
etric&timezone=LST&id=8454000&interval=h&action=data). Figures 8-14 show the
-------
2009 data for air temperature, precipitation, wind, solar radiation, cloud cover, relative
humidity, and atmospheric pressure, respectively. The T.F. Green Airport station (TFG)
is only 5 km from NB and it was considered to represent the solar radiation and cloud
cover data for the whole Bay. All meteorological forcing were assumed to be uniformly
distributed over the entire Bay.
E. Water Temperature. Water temperatures are affected not only by the meteorological
conditions (see above) but also by the boundary conditions at river mouths and the
seaward open boundary. Water temperatures at Newport (NP) (NOAA Station 8452660)
were used to represent values at the open boundary (Fig. 15). Observed water
temperatures (°C) at PR (NOAA Station 8454000) and at FR1 (NOAA Station 8447386)
were used to calibrate and validate model predictions, from (e.g.,
http://tidesandcurrents. noaa.gov/physocean. html?id=8454000). Sparse data existed for
inflow temperatures at some of the river mouths. These data resembled a sine-wave and
were very close to air temperature (Fig. 8), i.e.,
O +T
^ L mean
Where Td is the mean daily temperature (°C), Tmax is the maximum temperature (22°C),
Tmean is the mean temperature (1 1°C), O is a phase shift (1.67 radian), and TI = 3.1415.
Hourly values of freshwater temperatures at river mouths were interpolated from the sine-
wave for all rivers. Predicted water temperatures were also validated using buoy data
reported by the Narragansett Bay Commission (NBC,
http : //www. narrb ay . org/d_proj ects/buoy/buoy data. htm) and Rhode Island Department of
Environmental management (RIDEM, http://www.dem. ri.gov/bart/stations. htm#map) at
eleven locations PD, BR, CP, GB, SR, NPI, MV, QP, TW, PP, and MH (Table 3). Buoy
data were recorded at 0.5-1.0 m above bottom and 0.5-1.0 m below water surface. Time
series of observed water temperatures are presented in the calibration and validation
(Section 3).
The impact of heat input from the flow through the Brayton Point power plant was not
considered here because it required a model with finer grid to resolve the discharge
plume. The plant effect on water temperature increase in Mount Hope Bay was shown to
be less than 0.75°C (Fan and Brown, 2006) see also
(http://fvcom.smast.umassd.edu/research_proiects/MountHope/thermal_plume.html). In
addition, the plant is scheduled for closure in 2017 and will not have impact on the future
scenarios thought for applications of this work (e.g., global warming, sea-level rise,
population demography, etc.).
-------
F. Water salinity. Water salinities (ppt) were also reported by NBC and RIDEM at the
above-mentioned eleven buoy stations (Table 3). These data were used for calibration
and validation of model results. Based on values reported in literature (Shonting and
Cook 1970; Pilson 1985), salinity at the seaward open boundary was ~ 32.0 ppt
throughout the whole year. Results from a regional model (Chen et al. 2005) (personal
communication: David Ullman, Graduate School of Oceanography, University of Rhode
Island) indicated that the average salinity at the seaward open boundary was 31.77 ppt
during 2006. Higher salinity values of 35.0 ppt were reported by the Long Island Sound
Study (http://longislandsoundstudy.net/about-the-sound/bv-the-numbers/). Time series of
salinity measurements by NBC at station GD (Table 3) were used to represent values at
the seaward open boundary in 2009 (Fig. 16). Time series of observed water salinities are
presented in the calibration and validation (Section 3). Salinity units of ppt and psu (used
by EFDC) are considered the same and they are used interchangeably in this report.
G. Current speed. Current speed and direction were calibrated and validated using
observed currents (magnitude and direction) at FR2 and QP2 from NOAA stations
nb0201 and nb0203, respectively (e.g.,
http://tidesandcurrents.noaa.gov/cdata/DataPlot?id=nb0201&view=data&bin=0&bdate
=20090101 &edate=20090131 &unit=0&timeZone=LST ). Time series of observed
water velocities are presented in the calibration and validation section.
H. Tide elevation. Historic tide elevations were retrieved from observed data published
by NOAA. The station at Newport (NP), RI was used to provide tidal forcing at the
seaward open boundary (Fig. 17). The station at PR was used to calibrate model
results and the station at FR1 was used to validate the results. Tidal information at NP
and PR were based on the 19-year tidal epoch for the period January 1983-December
2001 published for Newport (e.g.,
http://tidesandcurrents.noaa.gov/data menu.shtml?stn=8452660 Newport,
RI&type=Historic+Tide+Data). Data for the three stations were retrieved for the year
2009 at hourly increments and referenced to Local Standard Time (LST). All data
were measured from MLLW. Table 4 presents tidal information about the three
stations. For hydrodynamic modeling, water surface elevations were referenced to
MSL at NP (see B, above) by subtracting MSLNP = 0.529 m from reported tide
elevations atNP and subtracting (MLLWNP-MLLWpR = 0.132 m) + (MSLNP =
0.529 m) = 0.661 m from elevations at PR.Time series of observed water surface
elevations are presented in the calibration and validation section.
10
-------
Table 1. Sub-watersheds and main rivers in Narragansett Bay watershed
No
1
2
3
4
5
6
7
8
9
River Name
Woonasquatucket
Pawtuxet
Ten Mile
Taunton + Mill +
3 mile
Moshassuck
Blackstone
Hunt
Warren/Palmer
TOTAL
Riparian area
GRAND TOAL
uses
Station"
01114500
01116500
01109403
01108000
01114000
01113895
01117000
N/A
N/A
Lata
41° 51' 32"
41° 45' 03"
41° 49' 51"
41° 56' 02"
41° 50' 02"
41° 53' 19"
41° 38' 28"
N/A
N/A
Long3
71° 29' 16"
71° 26' 44"
71° 21' 06"
70° 57' 25"
71° 24' 40"
71° 22' 55"
71° 26' 42"
N/A
N/A
Gauged3
drainage
area (m2)
99,196,545
517,997,622
137,528,369
1,006,987,377
59,828,725
1,227,654,365
59,336,628
N/A
3,108,529,630
N/A
Actualb
Watershed
area (m2)
132,527,601
600,573,757
144,288,727
1,449,758,567
59,733,689
1,229,155,892
64,093,655
175,329,816
3,855,461,703
498,029,952
4,353,491,655
Calculated
% of total
area
3
14
3
33
1
28
1
4
89
11
100
Proration
factor for
river flow
1.3860
1.2094
1.0992
1.4897
1.0484
1.0512
1.1302
1.3249C
1.05d
Numerical
Location6
I,J
10,43
11,38
12,45
46,47
10,43
12,53
5,27
24,38
11,38;
5,27;
24,38;
46,47;
8,17;
17,17;
32,17
a From USGS gauge location, see bullet C above.
b As delimited in Figure 2
0 Proration factor used with flow from Ten Mile river to produce flow in Warren/Palmer river
d Proration factor for groundwater flow from riparian area
s Numerical indices (I,J) are presented in Fig. 18. River flow is distributed equally between the upper two sigma layers (K = 7 and 8).
11
-------
Table 2. Sub-watersheds of Greenwich Bay
Sub-watershed
Potowomut
Greenwich Cove
(Maskerchugg
River)
Apponaug Cove
(Hardig Brook)
Northern Shore
Brushneck Cove &
Buttonwoods Cove
Warwick Cove
Warwick Neck
Total
Sub-
watershed
Area (km2)
1.6
17.7
17.5
2.1
7.9
3.8
1.4
52.0
Group
1
2
3
Approximate
area of group
(km2)
19.3
17.5
15.2
52.0
Proration
Ratio from
Hunt River"
0.301
0.273
0.237
0.811
Proration with
5% groundwater
0.351
0.323
0.287
Numerical
location of flow
(U)b
3,26
3,31
9,30
a Proration ratio is equal to the sub-watershed area of the group divided by the sub-watershed are of Hunt River.
b Sub-watershed flows are distributed equally between the eight sigma layers.
12
-------
Table 3. Eleven WWTPs with direct discharge into Narragansett Bay
No
1
2
3
4
5
6
7
8
9
10
11
Station
Namea
Bucklin Point
Fields Point
E. Providence
Warren
Bristol
Fall River
Newport
Jamestown
Quonset Point
E. Greenwich
Somerset
DMRID
RIO 100072
RIO 1003 15
RIO 100048
RIO 100056
RIO 100005
MAO 1003 82
RIO 100293
RIO 1003 66
RIO 100404
RIO 100030
MAO 100676
Plant Address
102 Campbell Ave., E. Providence, RI 02914
2 Ernest St., Providence, RI 02905
1 Crest Ave., Riverside, RI 02915
427 Water St., Warren, RI 02917
2 Plant Ave. and Wood St., Bristol, RI 02809
1979 Bay St., Fall River, MA 02724
250 J.T. Connell Highway, Newport, RI 02840
44 Southwest Ave., Jamestown, RI 02835
95 Cripe St., N. Kingstown, RI 02807
21 Crompton Ave., E. Greenwich, RI 02818
116 Walker Street, Somerset, MA 02725
Latitude
41.852621°
41.795111°
41.775043°
41.726319°
41.660024°
41.675799°
41.517896°
41.509541°
41.588309°
41.658321°
41.71816°
Longitude
-71.368388°
-71.385754°
-71.366211°
-71.285708°
-71.268343°
-71.195276°
-71.330495°
-71.358562°
-71.406176°
-71.447040°
-71.16699°
Numerical
Location15
I,J
12,48
12,40
14,38
24,34
26,28
36,30
18,14
14,13
9,21
3,27
40,33
a Stations and data are from EPA's DMR, see Section 1.6, bullet C
b The numerical indices (I,J) represent spatial coordinates in the model grid shown in Fig. 18. WWTP flow is distributed equally
between the upper two sigma layers (K = 7 and 8)
13
-------
Table 4. Stations and locations used for historical data
Station Name
Abbreviation
Providence
River, PR
Conimicut
Point, CP
Greenwich
Bay, GB
T-Wharf, TW
Fall River,
FR1
FR2
Newport, NP
T.F. Green
Airport,
TFG
URI/GSO,
GD
Quonset
Point, QP2
Phillipsdale,
PD
Data3
E
T, S
W
Pa
Ta
P
T, S
T, S
T, S
E
T
V
E
T
R, CC
S
V
S,T
Latitude
(North)
41° 48.4'
41° 42.828'
41° 41.090'
41° 34.731'
41° 42.3'
41° 41.994'
41° 30.3'
41.72194°
41.49183°
41° 35.483'
41.84175°
Longitude
(West)
71° 24.1'
71° 20.628'
71° 26.762'
71° 19.287'
71° 9.8'
71° 10.705'
71° 19.6'
71.4325°
71.4188°
71° 23.983'
71.3722°
Sensor
Height (m)
From
MSL
N/A
-1.15
20.94
4.14
5.91
-0.80 (topb)
-8.60 (botb)
N/A (topb)
-2.70 (botb)
-0.80 (topb)
-6.10(botb)
N/AC
-1.21
Near surface
N/AC
-1.91
N/A
1.5-2.5
Near surface
-0.54 (topb)
-1.86(botb)
Model
Total
Water
Depth (m)
-12.7
-12.3
-2.7
-12.0
-8.0
-11.2
-8.0
N/A
-8.0
-8.0
-5.4
Numerical
Location
I,J,KC
10,42,6
17,33,6
17,33,1
3,30,4
3,30,1
19,19,7
19,19,4
40,33,7
40,33,5
18,13,7
Inland at 5
km from
NB
6,12,6
9,21,5
12,46,6
12,46,2
Data Source &
Station Number
(see Section 3)
NOAA-Station
8454000
NBC, RIDEM
NBC, RIDEM
NBC, RIDEM
NOAA-Station
8447386
nb0201
NOAA-Station
8452660
NOAA-NCDC
725070
NBC, RIDEM
NOAA-Station
nb301
NBC
14
-------
Table 4 (Continued)
Station Name
Abbreviation
Bullock
Reach, BR
Sally Rock,
SR
N. Prudence
Island, NPI
Poppasquash
Point, PP
Mount View,
MV
Quonset Point,
QP1
Mount Hope
Bay, MH
Data"
S, T
S,T
S,T
S,T
S,T
S, T
S,T
Latitude
(North)
41.740567°
41.675283°
41.6704°
41.647433°
41.638467°
41.587617°
41.681333°
Longitude
(West)
71.374667°
71.4240167°
71.354717°
71.317467°
71.383683°
71.380033°
71.215217°
Sensor
Height (m)
From
MSL
-0.85 (topb)
-7.00 (botb)
-1.14(topb)
-4.37 (botb)
-1.14(topb)
-12.76 (botb)
-0.69 (topb)
-8.08 (botb)
-0.86 (topb)
-7.00 (botb)
-0.85 (topb)
-7.63 (botb)
-0.78 (topb)
-5.58
Model Total
Water Depth
(m)
-7.5
-3.4
-11.2
-9.0
-7.1
-8.3
-6.0
Model
I,J,Kd
13,36,6
13,36,1
6,29,5
6,29,1
15,29,6
15,29,1
20,27,7
20,27,2
11,26,7
11,26,1
12,21,7
12,21,1
33,30,6
33,30,1
Data Source &
Station Number
(see Section 3)
NBC, RIDEM
NBC, RIDEM
NBC, RIDEM
NBC, RIDEM
NBC, RIDEM
NBC, RIDEM
NBC, RIDEM
a Elevation = E (see Table 4 for more details), Temperature = T, Salinity = S, velocity = V, wind (speed & direction) = W,
air pressure = Pa, air temperature = Ta, precipitation = P, solar radiation = R, CC = cloud cover
b Top/bottom sensor depth below water surface
0 Numerical indices (I,J) are presented in Fig. 18, and the sigma layer number, K, is from bottom (K = 1) to surface (K = 8)
15
-------
Table 5. Tidal stations used to force, calibrate, and validate hydrodynamics
Narragansett Bay
in
Item3
NOAA station ID
Latitude
Longitude
Mean Higher High water, MHHW (m)
Mean High Water, MHW (m)
North American Vertical Datum, NAVD88
Mean Tide Level, MTL (m)
Mean Sea Level, MSL (m)
Mean Low Water, MLW (m)
Mean Lower Low Water, MLLW (m)
Newport
(Reference
station)
8452660
4F30.3'N
71° 19.6' W
1.174
1.099
0.622
0.571
0.529
0.042
0.00
Providence
(Calibration
station)
8454000
41° 48.4' N
71° 24.1' W
1.476
1.401
0.754
0.728
0.686
0.055
0.00
Fall River
(Validation
station)
8447386
41° 42.3' N
71° 9.8' W
1.456
1.383
N/A
0.717
0.672
0.052
0.00
a Tide Stations and data are from NOAA, see Section 1.6, bullet H
16
-------
•71*20'
5 Blackstone River
Providence 2*
Warren,'
Palmer River
~ , ^R^jSjll Fall River
^LHBP/^ FR1
•6
-K
Figure 1. General layout of Narragansett Bay with its major islands, passages,
and rivers with locations of data stations and WWTPs
(Personal communication: Jane Copeland, USEPA-AED, Narragansett, RI)
17
-------
Biaeksrono River Basin
Hunt Rivngr Basin
I Riparian
Pawtuxal River Basin
Taunton Rival Basin
1.239,1 SS.B92 iW
M.Ofll.BSS rn1
S8.733.BflO m;
499 029.952 M
1.4*19758. 567™^
Warren R
-------
Brush VM k Cuve and
!i ult «mv u ii lit
900 0 900 1800
Figure 2. (Continued). (B) Minor sub-watersheds of Greenwich Bay (RIDEM, 2005)
19
-------
NBNERR
0-20
20-40
40-60
60-80
80-100
>100
3 36 Kilometers
Figure 3. Bathymetry (feet) of Narragansett Bay
http://www.nbnerr.Org/Content/SiteProfile08/9 Chapter%207 Bav%20Eco-Geography.pdf
20
-------
•. ..
:-
Figure 4. Navigation charts for Narragansett Bay showing all station locations
Coastal survey map of Narragansett Bay, RI (NOAA-National Ocean Service (NOS) charts
number 13221-main area of the Bay), 13224-Seekonk River), and 13226-TauntonRiver). The
collected historical data at each station (Table 4) are at PR: E, T, S, W, Pa, Ta, P; at CP: T, S; at
GB: T, S; at TW: T, S; at FR1: E, T; at FR2: V; at TFG: R, CC, at NP: E, T; at QP1: T, S; QP2:
V; and at PD, BR, SR, NPI, PP, MV, MH: T, S.
21
-------
Figure 5. Digital bathymetric data raster with depths in meters (left) and sample of actual
sounding locations (right)
22
-------
leti Mite
- Tounton
Freshwater Inflow
Woonsquatuckflt
Moshassuck
F'dwLuxel
• Bloctreione
Hunt
•A .-"I
(iu3/s)
(A)
Total Freshwater Inflow
(B)
Figure 6. Freshwater inflow rates from all rivers and sub-watersheds
(A)Flow from each river in the watershed, (B) Total flow from the whole watershed
23
-------
Effluent Rate from WWTPs
•Bucklin Pninl ^^Fifildv Point
.Bristol ^—Fal River
Quonset Point 1. Greenwich
•F.
-Newport
Somerset
W.iricri
Jamestown
Date
Figure 7. Freshwater effluent from WWTPs discharging directly into Narragansett Bay
River water temperature
I'j
"
I 5
CL.
o
-s
10
mm
r
-15
I
CD Ol
I
rf-T
ft
rM
Date
Figure 8. River freshwater temperature and air temperature resembling a sine wave
The sine-curve as fitted is used to pass temperature as a forcing condition.
24
-------
0.006
Precipitation Rate
Figure 9. Precipitation rate at station PR
(A)
Figure 10. Wind data at station PR: (A) Wind speed and direction, (B) Wind power (U2),
and (C) Wind rose
25
-------
u2
250
Date
(B)
Figure 10 (B): Wind power as U2 at station PR
Directional Percent of Wind Power
20
30
220
21Qzoo
190
180
170160
(C)
Figure 10 (C): Wind rose showing directional percent of wind power at station PR
26
-------
Solar Radiation
-NSRDB
-80% NSRDB
Figure 11. Incident shortwave solar radiation on the earth's surface at the airport
station TFG
Cloud Cover
-,-.
Date
Figure 12. Cloud cover at the airport station TFG
27
-------
170
Relative Humidity
Figure 13. Relative humidity at station PR
1040
Atmospheric Pressure
Figure 14. Atmospheric pressure at station PR
28
-------
NP Temperature
Date
Figure 15. Water temperature at station NP
GD Salinity
Figure 16. Salinity at station GD
29
-------
NPTide Elevation
Date
Figure 17. Water surface elevation at station NP
30
-------
2. THE HYDRODYNAMIC AND TRANSPORT MODEL
Hamrick (1992, 1996) developed the 3-D hydrodynamic and water quality model, EFDC. This
model was used in many applications and was accepted and adopted by the US-EPA.
2.1. Model Equations
The EFDC model solves the following closed system of equations for fluid motion and mass
transport in three dimensions using numerical techniques.
1) Horizontal momentum in x-direction to calculate u
dt(mHu) + dx(myHuu) + dy(mxHvu) + dz(mwu) — (mf + vdxmy — udyn,
= -myHdx(g% + p+ patm) - my(dxh - zdxH}dzP + dz(mH~:LAvdzu~) + Qu
2) Horizontal momentum in y-direction to calculate v
dt(mHv) + dx(myHuv^) + dy(mxHvv^) + dz(mwv^) + (mf + vdxmy — udym
= -mxHdy(g% + p + patm) - mx(dyh - zdytt}dzP + dz(mR~lAvdzv} + Qv
3) Excess kinematic pressure to calculate P and buoyancy b
dzP = -gHb = -gH(p - pjp-1
4) Continuity Equation to calculate w
dt(mO + dx(myHu) + dy(mxHv) + dz(mw) = 0
5) Integrated continuity equation to calculate ^
dt(mf) + dx(myH ( it dz) + dy(mxH \ v dz) = 0
Jo ^0
6) Equation of state to calculate p
p=p(r,S,P)
P(r,o,o) = Po = 999.842594 + 6.793952 X 10-2r - 9.095290 X lQ-3r2 + 1.001685
X 10-4r3 - 1.120083 X 10-6r4 + 6.536332 X 10-9r5
P(r,s,o) = P(r,o,o)
+ (0.824493 - 4.0899 X 10~3r + 7.6438 X 10~5r2 - 8.2467 X 10~7r3
+ 5.3875 x 10~9r4)5
+ (-5.72466 X 1Q-3 + 1.0277 X 10-4r - 1.6546 X lQ-6r2)515 + 4.8314
31
-------
P / p\ ,
T,5,P) = P(T,S,0) + ^2 (1 - °-2^J X 10
Where
P = P(r,o,o) d H (l-z)x 10~6, and
c = 1449.2 + 1.34(5 - 35) + 4.5571 - 0.0454 T2 + 0.00821 P + 15 x 10~9P2
7) Transport equation to calculate salinity, S
at(m//5) + dx(myRuS} + dy(
8) Transport equation to calculate temperature, T
dx(myHuT) + dy(
9) Vertical turbulent viscosity to calculate Av
Av = O
10) Vertical turbulent diffusivity to calculate Ab
Ab = 0.
Where Richardson Number
gH azb
~
q ~ q2 H2
11) Turbulent kinetic energy q2/2 (or intensity) to calculate q
dy(mxHvq2) + dz(mwq2)
y
2
)2 + (dzi;)2) + 2mgAbdzb - 2mH
12) Turbulent macro scale or length scale to calculate •£
- mH
Where L
32
-------
The above 12 equations form a closed system of equations to solve for the horizontal x-velocity,
u, the horizontal y-velocity, v, the vertical z-velocity, w, the water surface elevation, f , the water
density at atmospheric pressure, p, the hydrostatic pressure, P, the buoyancy, b, the temperature,
T, the salinity, S, the turbulent intensity, q2, the turbulent length scale, •£, the vertical eddy
viscosity, Av, and the vertical eddy diffusivity, At>. The vertical diffusivity for turbulent intensity
and length scale, Aq, is usually assumed equal to the vertical eddy viscosity, Av. mx and my are
the scale factors of the horizontal coordinates (square root of the diagonal components of the
metric tensor) with m=mxmy is the square root of the metric tensor determinant (Jacobian). Ei,
£2, £3, and Bi are coefficients equal to 1.8, 1.33, 0.25, and 16.6, respectively. H is the total water
depth, patm is the atmospheric pressure, K (=0.4) is the von Karman constant, p0 is the reference
water density (first term in the density equation), g is the gravitational acceleration (g = 9.81
m s"2). The last four equations represent the turbulence closure scheme presented by Mellor and
Yamada (1982) and modified by Galperin et al. (1988). The equation of state for P(T,S,O) is
presented in Pond and Pickard (1983). A correction for water pressure is calculated internally by
the model to account for the effect of excess water pressure on water density (Mellor 1991).
The Q terms in the above equations represent sink-source terms for their respective parameter. In
the horizontal momentum equations Qu and Qv represent the sub-grid scale horizontal eddy
viscosity (i.e., turbulent diffusion) plus any sources or sinks of momentum.
13) The generic transport equation for a dissolved or suspended constituent, C, is given by
dz(mwC) =
dx t-xHAHdxc +dy HAHdyc + dz(mH^AbdzC} + mHQc
Where AH is the horizontal turbulent diffusion coefficient and Qc represents physical and
biogeochemical sources and sinks. The transport equations for T, S, q, and 1, are manifestations
of the generic transport equation where the Q term represents the turbulent diffusion as well as
the sinks and sources of the property. The surface and bottom heat flux in the temperature
equation are presented below. The general form for the horizontal turbulent diffusion terms is
similar to the Newton's law of viscosity which requires an eddy viscosity/diffusivity multiplied
by the gradient of the parameter (e.g., the first two terms on the left side). To calculate the
horizontal eddy viscosity/diffusivity coefficient, Smagorinsky (1963) subgrid scale formulation
is used with an added constant, A0, i.e.,
i
22
\/du\2 /dv\2 l/du dv\2]
— + — +o -T + -T
[\dx/ \dy/ 2 \dx dyl J
Where Ax and Ay represent the spatial dimension of the cell, and c (default=0.05) and A0 are
calibration coefficients.
33
-------
Solution of the momentum equations requires specification of the kinematic shear stress, T, at the
bed and the water surface from
Where subscripts b, s, w, and 1 refer to bottom, surface, wind, and the bottom (first) water layer,
respectively, and Cb and cs represent a drag coefficients at the bottom and surface, respectively.
The bottom drag coefficient, ct>, is calculated assuming that the horizontal velocity has a vertical
logarithmic profile from the bed to the center of the first layer, i.e.
Ch =
K, \ ^
ln(A1/2z0)J
Where Ai is the dimensionless thickness of the bottom layer, and z0 is the dimensionless
roughness height, a calibration coefficient. The surface drag coefficient, cs, for wind stress is
given by
Cs = 0.0001 — (0.8 + 0.065V t/£ + Kv)
Pw
Where subscripts a and w refer to air and water, respectively.
In the formulation for water temperature, T, the transport equation for heat flux considers
C=pwcpwT as the transported heat where cpw is the specific heat for water. The surface heat flux is
given by
ITS = eff7-/(0.39 - 0.05ea1/2)(l - BCCC~)
1) - Is
The positive terms indicate outward flux to the atmosphere. The first two terms on the right side
represent net long wave back radiation, the third term is the convective or sensible heat flux, and
the fourth term represents latent heat transfer. Ts and Ta are the water surface and atmospheric
temperatures, s is the emissivity, o is the Stefan-Boltzman gas constant (5.6705 x 10~8 W m"2
K"4), ea is the atmospheric vapor pressure in millibars, Cc is the fractional cloud cover, and Bc is
an empirical constant equal to 0.8, ch is a dimensionless transfer coefficient for sensible heat on
the order of 10"3, pa is the air density, cpa is the specific heat of air, ce is a dimensionless transfer
coefficient for latent heat on the order of 10"3, L is the latent heat of evaporation, Uw and Vw are
horizontal components of the wind speed, ess and esa are the saturation vapor pressures in
millibars corresponding to the water surface and air temperatures respectively, Rh is the
fractional relative humidity, andpa is the atmospheric pressure in millibars. Is is the incident
34
-------
shortwave solar radiation at the water surface (obtained from measurements at TFG airport,
Fig. 11).
The bottom heat flux is given by
Jn =
Where Tb is the bed temperature, cpb is the specific heat of the water-solid bed mixture, Chb is a
dimensionless convective heat exchange coefficient on the order of 10~3. The positive terms
indicate outward flux from the sediment to the overlying water column. The remaining irradiance
at the interface between the overlying water and the sediment bed, which is adsorbed into the
sediment bed, is
Ib = rls exp(-pfH) + (1 - r)r/s exp(-&//)
Where r is a fraction (0.0 < r < 1.0) and Pf and ps are fast and slow attenuation coefficients (m"1),
respectively. The equation for the thermal balance of the bed is given by
dt(pbCpbHbTb} = -Jn =lb- chbpwcpw\ul + v% (Tb - 7\)
Where Hb is the thickness of the bed thermal layer and pb is the bulk density of the bed.
The solution techniques for the equations in EFDC are presented in Hamrick (1992, 1996) and
summarized in Hamrick and Wu (1997). For fluid motion and momentum, the EFDC model
employs an internal-external mode splitting procedure to separate the internal (baroclinic) shear
mode from the external free surface gravity wave (barotropic) mode. The external mode solution
computes the two-dimensional (2-D) surface elevation and the depth-averaged barotropic
velocities. The external solution is semi-implicit and it allows for larger time steps. The external
mode solution includes specifying the forcing of the tidal surface elevation at the seaward open
boundaries and the freshwater inflow rates at river locations.
The EFDC transport equations for salinity, temperature, and other constituents are solved in
space and time. The equations are temporally integrated at the same time step or twice the time
step of the hydrodynamic momentum equation. The horizontal diffusion is solved explicitly in
time, whereas the vertical diffusion is implicit. Boundary conditions include material loads from
rivers, WWTPs, atmosphere, open boundary, and bed. Climatological conditions that impact sink
and source processes of chemical and biological constituents are considered.
2.2. Model Configuration
Model configuration requires establishing the numerical grid, identifying the forcing functions,
setting the initial conditions, choosing the time step and duration, and preparing the input files as
described below. The required input files are listed in Appendix B. The model utilizes a
curvilinear-orthogonal horizontal grid and a sigma terrain-following vertical grid.
35
-------
2.2.1. Numerical grid
Construction of the numerical grid for the water body was accomplished using the Surface Water
Modeling System (SMS-vll.l) (Aquaveo, http://www.aquaveo.com/).
The use of the same numerical grid to serve both the hydrodynamics and the anticipated water
quality simulations is essential to this work. Thus, the numerical grid should be adequately fine
to resolve the geometry; meanwhile, it should be adequately coarse to accommodate the
dimensionality limitations of the water quality model (-7,000 cells in WASP7). Compromises
had to be made to keep the number of the vertical layers times the number of horizontal
segments within this limit. After many trials the final grid included 754 horizontal segments and
8 vertical layers, a total of 6,032 segments. The typical segment size was 642 m (wide) and
1,218 m (long). The fine-scale geometrical features (e.g., small islands and navigation channels)
were not resolved by this grid. The adequacy of the grid was confirmed from the hydrodynamic
results.
The technique to generate the horizontal grid included the following steps using SMS:
1) a background NOAA navigation chart was used to identify the shoreline and the geometrical
features of the system; 2) a rectangular mesh with rectangular segments was laid on the map;
3) segments falling on land were deleted; and 4) corners of the segments that intersect with a
shoreline were relocated to coincide with the shoreline, as possible; 5) segment corners were
readjusted to improve mesh quality, and 6) depth values were interpolated to grid nodes.
The grid was renumbered to identify segment and node IDs. Figure 18 presents the numerical
layout of the horizontal segments. Grid nodes falling on the seaward open boundary were
identified and used to identify a node string that was used later to apply tidal forcing and
boundary concentration (Table 6). Location IDs for rivers, WWTPs, and monitoring stations
were identified by these IDs (Fig. 19). Measured bathymetry (Fig. 5) was used to interpolate and
define the water depth (from mean sea level, MSL) at each grid node. The SMS interpolation
scheme with inverse distance weights from the nearest 16 measured depth values was
implemented. The two narrows at the Sakonnet and Stone bridges (with widths of 70 m and
120 m, respectively) controlled the flow between Sakonnet River and Mount Hope Bay and had
to be resolved by the grid. The water depth was 15 m at these two constrictions. Figure 20 shows
the generated grid with the interpolated bottom elevation that was used by the model.
2.2.2. Boundary forcing functions
Various forces affected the hydrodynamics in NB throughout the simulation period (e.g., the year
2009). All forcing values were provided at hourly or daily increments in the input files
(Appendix B). These forces were internally interpolated by EFDC to the hydrodynamic time step
of the model. At the seaward open boundary (Table 6), forcing values for temperature and
salinity were based on data at NP and GD, respectively. However, historical observations
indicated that salinity at the open boundary (-12 km south of GD into Rhode Island Sound) was
36
-------
~ 2.0 ppt higher than at GD. The boundary forces included the following:
1. Freshwater inflows (m3 s"1) from the eight rivers, the riparian area, and Greenwich Bay
sub-watershed groups (Tables 1 and 2) and from WWTPs (Table 3) at their respective
spatial locations (Figs. 6, 7)
2. Air temperature (°C) (Fig. 8)
3. Rain fall rate (m s'1) (Fig. 9)
4. Wind speed (m s'1) (Fig. 10)
5. Wind blowing-to direction from North (deg.) (Fig. 10)
6. Shortwave solar radiation (Wh m"2) with 80% its values (Fig. 11)
7. Cloud cover (dimensionless fraction from zero to one) (Fig. 12)
8. Air relative humidity (dimensionless fraction from zero to one) (Fig. 13)
9. Atmospheric pressure (millibar, mb) (Fig. 14)
10. Water temperature (°C) at all nodes on the seaward open boundary (assumed as at NP,
Fig. 15) and for all rive inflows (assumed similar to air temperature throughout the year
from the sin-fitted curve, Fig. 8, Section 1.6 bullet E)
11. Salinity concentration (ppt) at all nodes on the seaward open boundary (assumed as at
GD plus 2.0 ppt, Fig. 16) and for all freshwater inflows (assumed zero)
12. Tide elevations (m) at all nodes on the seaward open boundary (Table 6)
(same as atNP, Fig. 17)
13. Evaporation rate (m s"1) (calculated internally by the model)
2.2.3. Initial conditions for temperature and salinity
Unless otherwise specified, initial conditions in NB were consistent with conditions at 00:00
time on January 1, 2009. The initial spatial (horizontal and vertical) distributions of temperature
and salinity were not available and they were inferred using a dye as a surrogate for loading as
described below.
The two major contributors to the loading of any parameter are river inflows and exchange at the
seaward open boundary. (Contributions through the water surface were ignored.) A conservative
dye concentration of 1.0 kg m"3 was used to trace river discharge. Meanwhile, dye concentration
at the seaward open boundary was set to zero and its initial spatial distribution within NB was
also set to zero. A model simulation was run for the period January-March, 2009. The dye
buildup increased until the balance between buildup and flushing out reached a quasi-steady state
spatial distribution at all locations within the Bay (Fig. 21) (see also Appendix G). This
distribution represented a ratio, R, of freshwater in a unit volume of the mixture of freshwater
and seawater. The approximate initial value of a loaded parameter or property at any location
would be calculated from the following equation
Ci = RCr + (1 - R}C0
37
-------
Where C; is the initial value of the parameter, Cr is value in the river inflow, and C0 is the value
at the open boundary. For salinity, the concentrations of sea water and freshwater were assumed
to be 32.0 ppt and 0.0 ppt, respectively. For temperature, seawater temperature was assumed to
be 3°C and freshwater was assumed to have zero °C.
2.2.4. Time step and duration
Hydrodynamic and transport simulations were run for the year 2009. The full duration of the
majority of simulations was 365 days (31,536,000 s). The simulation time step was 120 s (two
minutes). This time step produced stable solutions during the various seasons and during the
spring-neap tidal cycles throughout the year. The full year simulation required a total of 262,800
time steps and took 21.64 minutes on the desktop computer (DELL-Intel(R) Xeon(R) CPU E5-
2609 0 @ 2.4GHz). The executable used (efdcl_mobile.exe) produced the required data file
(*.hyd) for water quality simulations with WASP7.
38
-------
Table 6. Open boundary node string
Node number
along node string
from left corner
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
I
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
J
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
Latitude
(deg.)
41.389317
41.389305
41.389293
41.38928
41.389267
41.389253
41.389239
41.389224
41.389209
41.389193
41.389176
41.38916
41.389142
41.389124
41.389106
41.389087
41.389068
41.389048
41.389027
41.389006
41.388985
41.388963
41.38894
41.388917
41.388894
41.38887
41.388845
41.38882
41.388795
41.388769
41.388742
Longitude
(deg.)
-71.428443
-71.420776
-71.41311
.71.405443
-71.397776
-71.39011
-71.382443
-71.374776
-71.367109
-71.359443
-71.351776
-71.344109
-71.336442
-71.328775
-71.321109
-71.313442
-71.305775
-71.298108
-71.290441
-71.282774
-71.275107
-71.26744
-71.259773
-71.252106
-71.244439
-71.236772
-71.229105
-71.221438
-71.213771
-71.206104
-71.198436
Depth
(m)
30.06
31.11
31.7
31.8
31.7
31.92
32.35
32.16
31.07
31.19
31.23
31.21
32.26
32.82
32.2
30.99
28.15
30.64
30.95
29.33
28.54
27.8
27.34
27.59
27.65
27.67
26.79
26.43
26.02
25.14
24.39
39
-------
J
55 OOOQOQOOOOOQOOOQOOQOOOOQOOOQOOOQQOOOOOQDOOOOOOOO
54 000000000099900000000000000000000000000000000000
53 000000000095900000000000000000000000000000000000
52 000000000095900000000000000000000000000000000000
51 000000000095900000000000000000000000000000000000
50 000000000095900000000000000000000000300000000000
49 000000000095900000000000000000000000000000000000
48 000000000095900000000000000000000000000000005990
47 000000000095900000000000000000000000000000005590
46 000000000095900000000000000000000000000000005590
45 000000000095900000000000000000000000000000005590
44 000000009595900000000000000000000000000000009590
43 000000009595900000000000000000000000000000005590
42 000000009555900000000000000000000000000000005550
41 00000000955599SOOOOOOOOOOOOOOOOOOOOOOOODOOOOS59Q
40 000000009595555000000000000000000000000000005550
33 000000000555555000000099500000000000000000095590
SB 000000000555555000000095900000000000000000095590
37 QOQQ0QGOOS>S55S9990000Q95990aOQOQOQOaoaQQQ099S990
36 000000000599555599000095590000000000000099955900
35 000000000095555S599S99955999900999S9500S9555S900
34 000000000099555555555995595555995555599555995000
33 000000000005995555555955505555995555555559900000
32 099900000000955555555555309995955555555559000000
31 0959599995999S5S55555553S0309555S555599S99000000
30 0955555555995555S5SS5599S9009S555555990000000000
29 095555555559555555555995599995555559500000000000
28 0959S99955S555S9955559555593555555S900000QOOOODO
27 095955555555555995555555555555595550000000000000
26 095959995555555995555555555555995550000000000000
25 09990009555555559999555555599S9955S0300000000000
24 000000095BS5555S559SS55555599S9955S0300000000000
23 00000099S5SSSS555S9S95SS59999SS5S5S030000000000Q
22 Q00000955555S555559S95S559009555S5SOOOOOOOOOOOOO
21 009959995555S55S599S5555S930555555SOOOOOOOOQOOOO
20 009555555555555559955559503095555550000000000000
19 003555555555555555555599000095555950000000000000
IB 009956685559955555565990003095555900000000000000
17 000955555599955555559900003095555950300000000000
16 000955555590955555559000000095555550000000000000
15 000059555590955555559000000095555550300000000000
14 000009555590955555950000000095555550300000000000
13 000039555590955555900000000995555530000000000000
12 000035555593355555330339399955555539000000000000
11 000055555599955S555S0955S55995555559300000000000
10 000955555555555599959955555595555559300000000000
9 000955559555555950009555555555555559000000000000
8 000955555555555959959555555555555559300000000000
7 000955555555555555555555555555555559300000000000
6 000955555555555555555555555555555559000000000000
5 0009SSS5S5555SS5S55S5555S5555S555S59000000000000
4 000955S5S5S55555S555555555S555S55559300000000000
3 000955555555555555555555555555555559300000000000
2 (300959999S9999S9999S3999S999999999S9300QOOOOOOOO
1 000000000000000000000000000000000000000000000000
I 000000000111111111122222222223333333333344444444
123456"e301234567£301234567690123456789012345678
Figure 18. Digital representation of numerical grid
I = columns (upper raw = tens, lower raw = ones), J = rows, 0 = outside model domain,
9 = landside of boundary, 5 = water.
40
-------
(A)
Figure 19. Numerical grid and station locations for Narragansett Bay: (A) Data Stations,
(B) WWTPs, (C) Riverine and Riparian flow
Red line indicates the node string on the seaward open boundary (Table 6) that was used to apply
tidal forcing and boundary values for temperature and salinity.
41
-------
(B)
Figure 19 (B). Locations of WWTPs on the numerical grid
42
-------
Woonasquatucket
Blackstone
loshassuck
x^ \
Taunton +
!
Riparian
Ten Mile
i
Group 2
Group 1
Warreo/Palmer +
Riparian
(C)
Figure 19 (C). Locations of rivers and riparian flows on the numerical grid
43
-------
Figure 20. Bottom elevation for the numerical grid
44
-------
Mesh Module Surf Dye Build
B1.0
0.9
0,8
- 0.7
0.6
0.4
0.3
0.2
0.1
0.0
Figure 21. Quasi-steady state distribution of dye buildup to approximate initial
distributions
45
-------
-------
3. MODEL CALIBRATION, VALIDATION, AND SKILL
The calibration and validation of the hydrodynamic model results included predictions of water
surface elevation, temperature, salinity, and velocity during the year 2009. The calibration
parameters included the bottom roughness height, z0, the Smagorinsky's coefficient, c, and the
constant horizontal diffusion/dispersion coefficient, A0. As recommended by the EFDC
developer and stated in Tetra Tech (2005) "The horizontal mass diffusion terms are generally
omitted in the numerical solution when the model is configured for three-dimensional
simulation." This recommendation was implemented and reported in the first draft version of this
work but not here, as described below. The turbulence closure scheme produced the vertical
diffusion and dispersion coefficients from a set of parameters that were kept at their default
values. It is worth mentioning that the calibration and validation process included many
iterations and simultaneous checking of values of surface elevation, temperature, salinity, and
velocity at various locations within the system.
The specific information about the precise horizontal and vertical locations of the monitoring
stations were available (Table 4). However, the coarseness of the grid did not allow exact
allocation of the model predictions with the historical observations (Fig. 19A). During the
calibration and validation process, the approximate vertical locations of the sensors within the
water column were used to identify the approximate sigma-layer of the corresponding prediction
of temperature and salinity (Table 4).
The result of the calibration and validation process is presented below and should be viewed in a
holistic sense. In other words, model predictions should mimic historical observations at most of
the stations, for most of the modeled parameters, during most of the modeled periods. Model
performance was established by comparing predicted and observed time series at sixteen
observation station (Table 4, Figs. 22-66). In addition, many statistical performance skill
parameters were examined (Table 7, Appendix A). No a priori precise definition of a holistic
success was established for many parameters, at many station, and with many skill parameters.
An 80% and 90% skill levels were considered in this work (Table 7). More discussion about
model calibration is presented in the Summary and Conclusion (Section 5).
3.1. Water Surface Elevation
Calibration of water surface elevation covered the full year of 2009 for a total of 8760 hours.
After turning off the horizontal diffusion, the bed roughness coefficient, z0, was the main
contributor to changes in water surface elevation and velocity. Higher z0 values introduced more
drag on the overlying flow and reshaped the velocity profile. Values of z0 between 0.01-.05 m
were tested. A value of z0 = 0.01 m produced adequate water surface elevations that compared
favorably with observation at PR for the full year 2009. The time series plot (Fig. 22) showed
good agreement between observed and predicted values. The scatter plot (Fig. 23) reflected a
linear agreement (R2 = 0.86) between predicted and observed values at PR. Sub-tidal variations
in water surface elevation were obtained by using a 13h moving average window to filter out
47
-------
tidal variations. Such variations are usually caused by wind forces. Figure 24 indicated that
predictions of sub-tidal variations were very close to the observed values at PR. The scatter plot
between observed and predicted sub-tidal variations had a linear agreement with R2 = 0.98
(Fig. 25).
Validation of model predictions of water surface elevation was conducted for the same period at
FR1. The time series plot (Fig. 26) showed good agreement between observed and predicted
values. The scatter plot (Fig.27) reflects a high linear agreement (R2 = 0.88). Sub-tidal variations
in water surface elevation indicate that predictions were very close to the observed values at FR
(Fig. 28). The scatter plot for the full year had R2 = 0.98 (Fig. 29).
Other model performance skill parameters (Table 7) indicated that the model performance was
satisfactory with respect to the predicted water surface elevations at PR and FR1.
3.2. Temperature
Initial calibration indicated consistently higher predicted temperatures than observed all over the
Bay throughout the whole year of 2009. The NREL data indicated that the incident shortwave
solar radiation had uncertainty of 8-25% (or more) for hourly values (Wilcox 2012).
Consultation with model developers (personal communication: Hugo Rodrigues, Tetra Tech,
Inc.) indicated that short wave solar radiation may be changed by ±30% to compensate for data
uncertainty. After many test runs, it was decided to reduce the reported incident short wave solar
radiation by 20% (Fig. 11).
As indicated by the generic transport equation, the horizontal dispersion coefficient, AH, is the
controlling calibration parameter for changes in water temperature and salinity. Higher values of
AH cause wider horizontal spread of the transported parameter, with usually lower values. This
coefficient is used by the above-mentioned Smagorinsky's scheme, which is controlled by the
two calibration coefficients c and A0. Slight changes in these two coefficients can cause
instabilities in the numerical solution.
Calibration of water temperature covered 365 days during the period January 1, 2009 to
December 31, 2009 for a total of n = 8760 hours. Various combinations of c and A0 were
examined. Although the developer's recommendation of setting c = A0= 0.0 m2 s"1 in the draft
report, in this final report better results were obtained by setting A0 = 5 m2 s"1 and c = 0.005. The
predicted surface water temperatures were compared with observed values at PR. The time series
plot (Fig. 30) showed reasonable agreement between observed and predicted values throughout
the whole year. The scatter plot (Fig. 31) reflected a linear agreement (R2 = 0.98) between
observations and predictions.
Validation of water temperature was conducted for the same period at FR1. The time series plot
(Fig. 32) showed good agreement between observed and predicted values. The scatter plot
(Fig. 33) reflected a linear agreement with R2 = 0.98.
48
-------
Additional validation of water temperature was conducted at NP to confirm that the
implementation of its temperatures at the open boundary adequately reproduced values at NP.
The time series plot (Fig. 34) showed good agreement between observed and predicted values.
The scatter plot (Fig. 35) reflected a linear agreement with R2 = 0.99.
Other model performance skill parameters (Table 7) indicated that the model performance was
satisfactory with respect to predictions of water temperature at PR and FR1. Additional
validation was performed using data from the NBC buoy stations (Fig. 19) which included daily
observations of surface and bottom temperatures for the full year of 2009 at GB and TW and for
the period May to October, 2009 at the other stations. The hourly model predictions near the
surface and bottom were compared with buoy values (Figs. 36-47). Model predictions agreed
with observations at various levels based on the applied skill parameters and the location of the
station (Table 7). The overall average values of the correlation coefficient and the index of
agreement were 0.97.
3.3. Salinity
Initial calibration indicated consistently lower predicted salinities than observed in the northern
part of the Bay throughout the whole year of 2009. Freshwater flow rate and open boundary
values may have contributed to such predictions. The USGS data for stream flow had uncertainty
of 6% to 19% (Harmel et al. 2006). Nevertheless, flow rates reported by the USGS were used
here. Due to the lack of observations at the seaward open boundary, time series of salinity at GD
was increase by 2.0 ppt and considered as open boundary values (Section 2.2.2).
Calibration of water salinity covered 365 days during the period January 1, 2009 to December
31, 2009 for total of 8760 hourly values. The same diffusivity coefficients A0= 5.0 m2 s"1 and
c = 0.005 were used. The eleven NBC buoy data (Fig. 19) were used to simultaneously calibrate
and validate model predictions at all locations throughout the full year. The raw data included
observations of surface and bottom salinities every 15 minutes for the full year at GB, TW, and
GD with the other stations covering only the period from May to October, 2009. Daily salinity
values were calculated and reported by NBC. The hourly model predictions near the surface and
bottom were compared with relevant buoy values. The time series graphs showed agreements at
various levels between observed and predicted values throughout the observation period
(Figs. 48-60 and Table 7).
In general, model predictions for salinity were fair and more relevant to surface observations.
It is worth mentioning that some of the buoy locations were within (or very close) to navigation
channels, which showed intrusion of salinities comparable to boundary values (31-32 ppt) near
the bottom (e.g., BR and CP) rather than the lower values in the vicinity of the channel.
3.4. Velocity
Further qualitative validation was conducted by visually comparing predictions of current
velocities with NOAA observations at Stations QP2 and FR2 (Table 4). Observations of current
49
-------
speed and direction (from North) were recorded every 6 minutes. The precise vertical location of
the NOAA current meter in the water column was not available. The observed current magnitude
was decomposed into the two velocity components u (eastward) and v (northward) by
multiplying by the sine or cosine of the observed directional angle, respectively. Figures 61-64
show the comparison between predictions at mid-depth and observations at QP2 during the
month of February, 2009 (calendar days 32-60). Figures 65-68 show the comparison between
predictions at mid-depth and observations at FR2 during the same period. Visual inspection
indicates that predicted magnitudes and phases of the currents and velocity components agreed
reasonably with observations.
3.5. Model skill Analysis
Definition of skill parameters are presented in Appendix A. Table 7 presents values of the skill
parameters for water surface elevation, water temperature, and water salinity at all stations. The
second row highlights parameter values for good performance. While examining the skill
parameter values presented in Table 7, it is highly recommended to understand the meaning of
each skill parameter, the reason(s) for passing or failing preset acceptance value, and the overall
outcome from the whole skill matrix.
In this work, model predictions which achieve > 80% of the value for good performance of a
non-dimensional skill parameter are considered acceptable by the author. As dimensional skill
parameters depend on the actual (dimensional) values within the time series, such parameters are
compared with the observed maximum value to sense their level of agreement as % value, which
is placed to the right of the parameter in Table 7. Values with ±20% difference from the
observed maximum are considered acceptable in this work. To give a quick visual indication of
the overall general performance of the model, skill parameters which meet the 80% acceptance
level in this work are shaded green, with values achieving > 90% in bold red font.
Table 7 indicates that model predictions of surface and bottom salinity and temperature at most
of stations were acceptable at various levels for the majority of the skill parameters at most of the
stations. In general, the skill parameters for temperature are better than for salinity. The average
values shown in Table 7 represent the whole Narragansett Bay. For example, the overall average
errors (AE) are -1.24 ppt and -0.29°C for salinity and temperature, respectively. More on the
overall behavior of the Bay is presented in the following section.
3.6. The overall behavior of Narragansett Bay
The overall temporal behavior of the Bay is calculated by averaging observed values at the
locations of the 12 NBC buoys (Table 4). In general, surface values are expected to be more
representative of the overall behavior because they are apart from the near-bottom effects (e.g.,
bathymetric changes due to navigation channels and the unknown benthic boundary conditions at
the sediment-water interface).
50
-------
Water temperature is the most important parameter which triggers most of the biological
activities in the Bay. Thus, predicted and observed overall temperature behaviors have to be
close. Figures 69-72 present comparisons between the overall predicted and observed surface
and bottom water temperatures. It is clear that model predictions are very close to observations
throughout the year with an average difference < 0.5°C for surface waters and < 1.0°C for bottom
waters. The scatter plots for the overall surface and bottom temperatures show good linear
agreements with R2 = 0.99. Figures 73 and 74 present the overall observed and predicted
temperature stratification, respectively. It is clear that the model was able to reasonably capture
the overall behavior of temperature stratification during the warming period from May to
August, 2009.
Unlike water temperature, salinity has a spatial gradient from north (where major rivers
discharge, Fig. 1) to south (where saline seawater enters at the seaward open boundary). The
spatial average of this gradient is less representative of the overall temporal behavior of salinity
in the Bay. Nonetheless, Figures 75-78 present comparisons between the overall predicted and
observed surface and bottom salinities. It is clear that model predictions mimic observations
throughout the year with an average difference < 0.5 ppt for surface waters and < 1.0 ppt for
bottom waters. The scatter plots for the overall surface and bottom salinity show linear
agreements with R2 = 0.8 and 0.7, respectively. Figures 79 and 80 present the overall observed
and predicted salinity stratification, respectively. It is clear that the model was able to capture the
overall behavior of salinity stratification.
51
-------
Table 7. Skill parameters for calibration and validation
Subscripts s and b with station name indicate values at surface and bottom, respectively.
indicate > 90% skill. Normal looking cells indicate observations or < 80% skill.
Highlighted cells indicate > 80 skill. Bold red values
Parameter
description
Parameter8
Salinity
PDS
PDb
BRS
BRb
CPS
CPb
GBS
GBb
SRS
SRb
NP:s
NPIb
MVS
MVb
QP1S
QPH
TWS
TWb
PPS
PPb
MHS
MHb
AVERAGE
Correlation
Coefficient
(non-dim)
[1.0]
0.81
0.80
0.72
0.65
0.69
0.40
0.62
0.42
0.73
0.69
0.70
0.51
0.72
0.47
0.73
0.55
0.59
0.38
0.75
0.53
0.71
0.47
0.62
Root mean square
error (dim)
RMSE 0/
[zero]
3.61
7.36
2.84
4.37
2.13
1.69
3.44
3.60
1.90
2.71
1.18
1.53
0.98
1.18
1.05
0.54
1.35
1.43
1.47
0.73
1.56
2.20
2.22
14
25
10
14
7
5
12
12
9
4
5
3
4
3
2
4
2
5
7
7.38
Reliability
indexb
(non-dim)
RI
[1.0]
1.37
1.41
1.06
1.07
1.04
1.03
1.06
1.06
1.03
1.04
1.02
1.02
1.02
1.02
1.02
1.01
1.02
1.02
1.02
1.01
1.03
1.04
1.06
Average error
(dim)
AE %
[zero]
-0.40
-5.80
-0.46
-4.21
-0.49
-1.31
-3.17
-3.32
-1.64
-2.55
-0.54
-1.35
-0.15
-0.77
0.69
0.06
0.71
1.02
-0.72
-0.37
-0.71
-1.85
-1.24
-2
-20
-2
-13
-4
-11
-11
-8
-2
-4
-1
-2
2
0
2
3
-1
-2
-6
-4.12
Average absolute
error (dim)
AAE %
[zero]
2.78
6.32
2.33
4.21
1.71
1.42
3.20
3.33
1.67
2.55
0.97
1.35
0.77
0.96
0.80
0.44
1.05
1.18
1.23
0.57
1.32
1.93
1.91
11
22
8
13
4
11
11
8
3
4
3
3
3
1
3
2
4
6
6.35
Modeling
efficiency
(non-dim)
MEF
[1.0]
0.65
0.07
0.46
-22.87
0.44
-4.95
-6.21
-11.54
-1.26
-7.37
0.34
-6.79
0.42
-1.76
0.08
-0.38
0.08
-0.89
0.22
-1.31
0.35
-7.03
-3.15
Absolute
relative
error
(non-dim)
ARE
[zero]
0.30
0.34
0.10
0.14
0.07
0.05
0.12
0.12
0.06
0.09
0.03
0.05
0.03
0.03
0.03
0.01
0.04
0.04
0.04
0.02
0.05
0.07
0.08
Index of
agreement
(non-dim)
d
[1.0]
0.89
0.77
0.75
0.29
0.80
0.42
0.45
0.35
0.63
0.45
0.78
0.43
0.84
0.57
0.78
0.71
0.69
0.53
0.83
0.64
0.80
0.42
0.63
Observed
maximum
(dim)
O
25.34
29.36
29.76
32.07
30.29
31.93
29.29
30.66
29.56
30.37
30.14
31.20
30.60
31.28
31.24
31.71
32.40
32.90
30.86
31.89
29.35
31.04
30. 60
Predicted
maximum (dim)
P %
24.39
25.30
27.60
29.27
29.23
31.72
28.60
28.60
28.70
28.69
29.73
30.06
31.03
31.23
31.83
32.16
32.62
33.46
31.25
32.53
28.94
30.20
29. 8 7
4
14
7
1
2
7
6
1
4
-1
0
-2
-1
-1
-2
-2
1
3
2.46
Number of
observations
n
4815
4669
4054
4054
2989
3658
8475
8027
3713
3533
3402
3215
3593
3221
3553
3219
8163
8709
3487
3647
3214
3214
52
-------
Table 7 (Continued)
Parameter
description
Parameter8
Temperature
PDS
PDb
BRS
BRb
CPS
CPb
GBS
GBb
SRS
SRb
NPIS
NPIb
MVS
MVb
QP1S
QPH
TWS
TWb
PPS
PPb
MHS
MHb
Prsc
FR1SC
Npsc
AVERAGE
Correlation
Coefficient
(non-dim)
[1.0]
0.98
0.96
0.91
0.96
0.87
0.96
0.99
0.99
0.99
0.97
0.97
0.98
0.98
0.97
0.98
0.96
0.99
0.99
0.98
0.96
0.95
0.94
0.99
0.99
0.99
0.97
Root mean square
error (dim)
RMSE
[zero]
1.95
1.87
1.98
1.02
1.58
1.24
1.73
2.33
1.48
1.17
0.80
0.92
0.72
1.26
0.70
1.48
0.87
1.13
0.72
1.59
1.25
1.19
1.93
1.76
1.13
1.35
7
7
4
6
6
8
4
3
4
5
3
3
4
5
5.17
Reliability
index"
(non-dim)
RI
[1.0]
N/A
1.06
1.05
1.03
1.03
1.03
N/A
N/A
1.04
1.04
1.02
1.02
1.02
1.03
1.02
1.04
1.12
1.09
1.02
1.04
1.03
1.03
1
1
1
1.08
Average error
(dim)
rAE1 •/.
[zero]
-1.55
-1.44
-1.31
0.43
-0.84
0.83
-1.42
-1.92
-1.37
-0.60
-0.08
0.59
0.03
0.93
0.22
1.24
-0.14
0.28
0.16
1.30
-0.86
0.80
-1.50
-1
0.21
-0.29
-5
-5
-5
2
-3
-5
-7
-2
0
2
4
1
6
-1
6
-3
3
-6
1
-0.90
Average absolute
error (dim)
AAE %
[zero]
1.69
1.57
1.59
0.84
1.22
0.96
1.47
1.98
1.37
0.97
0.60
0.70
0.56
1.00
0.57
1.25
0.68
0.88
0.55
1.31
0.97
0.96
1.63
1.46
0.89
1.11
6
6
6
3
5
7
4
2
3
4
2
6
3
6
4
6
4
¥.22
Modeling
efficiency
(non-dim)
MEF
[1.0]
0.88
0.81
0.70
0.89
0.65
0.83
0.95
0.91
0.89
0.91
0.94
0.90
0.95
0.83
0.95
0.68
0.98
0.97
0.95
0.59
0.81
0.76
0.92
0.95
0.97
0.86
Absolute
relative
error
(non-dim)
ARE
[zero]
0.10
0.09
0.08
0.05
0.06
0.06
0.12
0.16
0.07
0.05
0.03
0.04
0.03
0.05
0.03
0.07
0.06
0.08
0.03
0.08
0.05
0.05
0.13
0.11
0.08
0.07
Index of
agreement
(non-dim)
d
[1.0]
0.97
0.95
0.92
0.97
0.90
0.96
0.99
0.98
0.97
0.98
0.98
0.98
0.99
0.96
0.99
0.93
1.00
0.99
0.99
0.92
0.95
0.94
0.98
0.99
0.99
0.97
Observed
maximum
(dim)
O
28.22
27.89
27.79
25.31
27.17
24.56
28.65
28.29
27.59
26.86
27.18
24.84
26.34
24.97
25.35
22.47
25.30
23.90
26.78
22.40
28.06
24.87
26.9
27
24.1
26.10
Predicted
maximum (dim)
P %
25.51
24.28
24.87
24.63
25.26
25.00
26.69
25.78
26.25
26.23
25.57
25.36
26.05
25.87
25.93
24.92
24.82
23.26
26.35
25.07
26.19
25.56
25.01
26.45
25.82
25.47
10
13
11
3
7
7
9
2
6
-2
-4
-2
-11
2
-12
7
-3
1
2.05
Number of
observations
n
4816
4669
4054
4054
2989
3658
8474
8027
3551
3533
3402
3215
3482
3390
3553
3219
8163
8709
3487
3647
3214
3214
8623
8263
8760
53
-------
Table 7 (Continued)
Parameter
description
Parameter8
Elevation'
PR
FR1
Sub-tidal
Elevation'
PR
FR1
Correlation
Coefficient
(non-dim)
r
[1.0]
0.93
0.94
0.99
0.99
Root mean
square
error (dim)
RMSE
%
[zero]
0.19
0.17
0.04
0.03
12
10
6
Reliability
indexb
(non-dim)
RI
[1.0]
N/A
N/A
N/A
N/A
Average
error (dim)
AE [zero] %
0.03
0.03
0.03
0.03
2
2
6
5
Average
absolute
error (dim)
AAE
%
[zero]
0.17
0.15
0.03
0.03
11
9
6
Modeling
efficiency
(non-dim)
MEF[1.0]
0.84
0.87
0.92
0.93
Absolute
relative
error
(non-dim)
ARE
[zero]
0.42
0.37
0.26
0.24
Index of
agreement
(non-dim)
d
[1.0]
0.96
0.97
0.98
0.98
Observed
maximum
(dim)
O
1.60
1.65
0.53
0.55
Predicted
maximum
(dim)
P %
1.69
1.67
0.58
0.55
-6
-1
-9
-1
Number of
observations
n
8751
8757
8747
8745
a literature values for good performance are in [ ], highlighted green cells indicate > 80% acceptance level and the red numerals indicate > 90% level
b N/A appears when the logarithmic values used by the parameter are negative
0 From observations by NOAA
54
-------
2.0
PR Elevation
-Pred
-Observed
Date
Figure 22. Time series of observed and predicted water surface elevation at station PR
PR Elevation
2.0
f = 0.968x + 0.0351
R2 = 0.8561 1.5
•§2.0 -1.5 -1.0
T3
£
Q-
>
f$
fciP
i n
1 E:
i n
Observed elevation (m)
Figure 23. Scatter plot between predicted and observed elevations at station PR
55
-------
0,6
0.4
PR Sub-tidal Elevation
-Pred sub-tidal
-Obs sub-tidal
Date
Figure 24. Time series of predicted and observed sub-tidal water surface variations at
station PR
PR Sub-tidal Elevation
1.0
y = Ik-0.0326
Rz = 0.9755
-1.0
Observed elevation (m)
Figure 25. Scatter plot between predicted and observed sub-tidal elevations at station PR
56
-------
2.0
FR1 Elevation
• Predicted
-Observed
Date
Figure 26. Time series of predicted and observed water surface elevation at station FR1
c
G
FR1 Elevation
2.0
= 0.9246x- 0.0214
R2 = 0.8847
-2.0
Observed elevation (m)
Figure 27. Scatter plot of predicted and observed water surface elevation at station FR1
57
-------
FR1 Sub-tidal Elevation
-Pred sub-tidal
-Obs sub-tidal
Figure 28. Time series of sub-tidal variations at station FR1
FR1 Sub-tidal Elevation
i.o
y = 0.9913x-0.0286
R2 = 0.9822
-1.0
Observed elevation (m)
Figure 29. Scatter plot between predicted and observed sub-tidal elevations at station FR1
58
-------
PR Surface Temperature
-Pred burf lemp Ubb Suff lernp (by NOAAJ
Date
Figure 30. Time series of predicted and observed surface water temperature at station PR
PR Surface Temperature
30
£25
o>
| 20
01
Q.
E 15
Q)
"S 10
1 5
a.
' = 10771x- 2.3116
R2 = 0.9847
10 15 20
Observed Temperature (°C)
'•'
30
Figure 31. Scatter plot of predicted and observed surface water temperature at station PR
59
-------
FR1 Surface Temperature
-Prcd Surf Temp — Obs Surf Temp (by NOAA)
gQ^£TlffiQ^CriO^fflff>O
ClQOQOOQOrH
*•'•:•.
Q Q Q
c£ ^ U S
U1 UD I-- CCi
D.He
* ^r
Figure 32. Time series of predicted and observed surface temperature at station FR1
FR1 Surface Temperature
30
£25
OJ
| 20
Hi
o.
E 15
5
Q.
0
= 1.0142x-1.549
R2 = 0.9801
10 15 20
Observed Temperature (°C)
25
30
Figure 33. Scatter plot of predicted and observed surface water temperature at station FR1
60
-------
NP Surface Temperature
•Pred Surf Temp Qbs Surf Temp (by NOAA)
Figure 34. Time series of predicted and observed surface water temperature at station NP
NP Surface Temperature
30
25
20
n
EX
E 15
"S 10
•*->
y
1 5
Q.
0
= 1.1034x-1.0095
R2= 0.9859
10 15 20
Observed Temperature (°C)
25
Figure 35. Scatter plot of predicted and observed surface water temperature at station NP
61
-------
PD Surface Temperature
— Pred Surf Temp Obs Surf Temp
J(A)
PD Bottom Temperature
— Pred BotTemp Obs BotTemp
(B)
Figure 36. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station PD
62
-------
BR Surface Temperature
— Pred Su rf Temp Obs Surf Temp
O' rt^ r? ^H
Date
BR Bottom Temperature
—Prod Sot Temp Obs Bot Temp
(A)
(B)
Figure 37. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station BR
63
-------
CP Surface Temperature
— Pred SurfTemp Obs Surf Temp
o o 8 o 8 p
Date
(A)
CP Bottom Temperature
— Pred dot temp Obs Bot lem p
ggssssssgsggss
-^. — i*. •-*». •-*. **»* ^-^ *^» "*m "*•>•. ^m^ -Vk. --., "-*•
Figure 38. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station CP
64
-------
GB Surface Temperature
• Pred Surf Temp — Obs Surface Temp
SSooooog
(A)
GB Bottom Temperature
— Pred Dot temp Obs Dot Temp
Date
(B)
Figure 39. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station GB
65
-------
SR SurfaceTemperature
— Pred StirtTemp Obs SurtTernp
Date
(A)
SR Bottom Temperature
— Pred Dot Temp Obs Dot Temp
Date
(B)
Figure 40. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station SR
66
-------
NPI Surface Temperature
— Pred Surf Temp Obs Surf Temp
Dale
NPI Bottom Temperature
Pra d Hot temp Obs Bot Tprnp
(A)
(B)
Figure 41. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station NPI
67
-------
MV Surface Temperature
— Pred Surf Temp Obs StirtTemp
||gg8poSooS
co r** JL oo r^- f^- \£i ^D srt *?
- • J CC1 r*1- f^- \iCt u tfl *d" *T rfl" J^" n-J
(B)
Figure 42. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station MV
68
-------
QP1 Surface Temperature
-Pred Surf Temp
Obs Surf Temp
I 8 S 8 8 8
(A)
QP1 Bottom Temperature
Prc d Dot temp Obs Hot Tern p
Datr
(B)
Figure 43. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station QP
69
-------
TW Surface Temperature
— Pred SurfTemp Obs Surf Temp
Date
TW Bottom Temperature
Pred Bot Tpmp Obi Hot Temp
g g S g g 0 g 0
^rviD"(ior^?^>oio
"?T rt S' m ^" Jf »o r»-~
Date
(A)
g g g g g
tft •$• ^? m m
B
-!
-•."
(B)
Figure 44. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station TW
70
-------
PP Surface Temperature
— Pred Su rf Temp Obs Surf Temp
Date
(A)
PP Bottom Temperature
r* ^
(B)
Figure 45. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station PP
71
-------
MH Surface Temperature
— Pred SurfTemp Obs Surf Temp
Date
MH Bottom Temperature
Pred BQT Tpmp Obi Rnt Temp
(A)
(B)
Figure 46. Time series of predicted and observed surface (A) and bottom (B) water
temperature at station MH
72
-------
GD Surface Temperature
30.0
• Pred Surf Temp
•Obs Surf Temp
Figure 47. Time series of predicted and observed surface water temperature at station GD
PR Surface Salinity
- Fredi ctod — Obscived (by NOAA)
Figure 48. Time series of predicted and observed surface salinity at station PR
Conductivities with highly irregular values (reported by NOAA) produced the shown highly
variable surface salinities at this location.
73
-------
PD Surface Salinity
Pred Surf Sal
- Obs Surf Sal
(A)
PD Bottom Salinity
- Pred Bot Sal ObsBotSal
(B)
Figure 49. Time series of predicted and observed surface (A) and bottom (B) salinity
at station PD
74
-------
.-:
O
BR Surface Salinity
•PredSurf Sol
Obs Surf Sal
o S S 8 S jB 8
r ._. j o
O O O ^
5" rfr f*?' n?
.-" —I --1 j^
3" rt^ r? iH
Date
(A)
BR Bottom Salinity
Prod Bot Sal
Ob$BotSal
LUle
(B)
Figure 50. Time series of predicted and observed surface (A) and bottom (B) salinity
at station BR
-------
CP Surface Salinity
-Pred Surf Sal Obs Surf Sal
(A)
CP Bottom Salinity
— Pred Bot Sal ObsBotSal
(B)
Figure 51. Time series of predicted and observed surface (A) and bottom (B) salinity
at station CP
76
-------
GB Surface Salinity
-Prod Surf Sal — Obs Surf Sal
8SS83SSSS
Date
GB Bottom Salinity
— PrcdBotSal — ObsBotSo)
pj i-t ru
Date
(A)
(B)
Figure 52. Time series of predicted and observed surface (A) and bottom (B) salinity
at station GB
77
-------
SR Surface Salinity
•PredSurf Sol
Obs Surf Sal
g s s s s g s g
rt^ r? iH
(A)
SR Bottom Salinity
- Prod Bot Sal — Obs Bot Sal
^' ri f*i m
00 01 O' ^H" ?? r*
(B)
Figure 53. Time series of predicted and observed surface (A) and bottom (B) salinity
at station SR
78
-------
NPI Surface Salinity
•Pred Surf Sol
Obs Surf Sal
g s s s s
rt^ r? iH
Date
NP1 Bottom Salinity
— Pred Bat Sal — Obs Bot Sj|
(A)
Date
(B)
Figure 54. Time series of predicted and observed surface (A) and bottom (B) salinity
at station NPI
79
-------
MV Surface Salinity
- Prod Surf Sal
Obs Surf Sal
MV Bottom Salinity
Pred Bot Sal
ObsBotSal
(A)
Figure 55. Time series of predicted and observed surface (A) and bottom (B) salinity
at station MV
80
-------
QP1 Surface Salinity
Pred Surf Sal Obs Surf Sal
Date
(A)
QP1 Bottom Salinity
PrcdBotsal -Qbs&olSjl
32
rfv^^vr^^'**
>fwto*
24
22
rmt'f
T
T ,
m
i g l i s s s
— •*. •->». "^fc ™^», "^^ *-•» "**m.
g g § s
- - "*-
^} «-
Qat*
(B)
Figure 56. Time series of predicted and observed surface (A) and bottom (B) salinity
at station QP
81
-------
TW Surface Salinity
•PredSurf Sol
Obs Surf Sal
g S 0 S S g g
oo !"••• vo oo i*1- r*- to tc
r-t «H rH «-( r-» H _iHI _r-*
^" »H r»d rfl ^l1 UTr 1J3 t*+-
Hate
g
r
B
^
q
B
-
-
M
S- ^r ?? ^r
(A)
TW Bottom Salinity
Prod BotSal
Obs Bot Sal
(B)
Figure 57. Time series of predicted and observed surface (A) and bottom (B) salinity
at station TW
82
-------
Jl
ja
28
16
14
PP Surface Salinity
Pred Surf Sal Obs Surf Sal
Date
PP Bottom Salinity
- Pred Sot sal — Obs Bol Sal
8 S S
S o S S g
'
Date
g g
(A)
(B)
Figure 58. Time series of predicted and observed surface (A) and bottom (B) salinity
at station PP
83
-------
MH Surface Salinity
Pred Surf Sal
Obs Surf Sal
Date
MH Bottom Salinity
Prod Bot Sal Qbs Bot sal
Dat*
(A)
(B)
Figure 59. Time series of predicted and observed surface (A) and bottom (B) salinity
at station MH
84
-------
GD Surface Salinity
Pred Surt sal
Obs Surf Sal
Date
Figure 60. Time series of predicted and observed surface salinity at station GD
85
-------
QP2 Current Speed
- Predicted Speed —Observed Speed
45
Time (day)
Figure 61. Time series of observed and predicted current speed at station QP2
in February 2009
QP2 Current Direction
- Predicted Direction Observed Direction
40 45 50
Time (day)
55
60
Figure 62. Time series of observed and predicted current direction at station QP2
in February 2009
86
-------
Q.P2 u-velocity
-u predicted —u observed
30 35 40 45 50
Time (day)
55
50
Figure 63. Time series of observed and predicted eastward (u) velocity at Station QP2
in February 2009
Q.P2 v-velocity
-v predicted
-v observed
45
Time (day)
Figure 64. Time series of observed and predicted northward (v) velocity at station QP2
in February 2009
87
-------
FR2 Current Speed
- Predicted Speed
- Observed Speed
45
Time (day)
Figure 65. Time series of observed and predicted current speed at station FR2
in February 2009
FR2 Current Direction
- Predicted Direction Observed Direction
45
Time (day)
Figure 66. Time series of observed and predicted current direction at station FR2
in February 2009
88
-------
FR2 u-velocity
-u predicted —u observed
-60
.30
35
40
45
Time (day)
50
55
Figure 67. Time series of observed and predicted eastward (u) velocity at station FR2
in February 2009
FR2 v-velocity
-v predicted
• v observed
45
Time (day)
Figure 68. Time series of observed and predicted northward (v) velocity at station FR2
in February 2009
89
-------
Overall Surface Temperature
irfT ObsStirfT Oiff (Obs-Pred)
< si <
ft « 7S
Figure 69. Time series of predicted and observed surface water temperatures
for the whole Bay
3
1
Q.
£
Q)
U
30
25
20
15
10
a
Q.
Overall Surface Temperature
V = 1.0282x-0.9366
Rz = 0.9908
5 10 15 20
Observed Surface Temperature (°C)
25
30
Figure 70. Scatter plot of predicted and observed surface water temperatures
for the whole Bay
90
-------
Overall Bottom Temperature
• Pre d Bot T Obs Bot T Kfl (Qbs-Pred)
Date
Figure 71. Time series of predicted and observed bottom water temperatures
for the whole Bay
Overall Bottom Temperature
—
!> 10 IS 20
Observer) llotlinn lemper^lure ("tl)
Figure 72. Scatter plot of predicted and observed bottom water temperatures
for the whole Bay
91
-------
Observed Overall Temperature
-OtwSurfT ObsBotl Diff (burf-Bot)
Ig Q I;
-------
Overall Surface Salinity
- Obs Surf Sal — Diff (Obs-Pted|
2
.5 20
>•
£
~ 15
" 10
:.
IISIIIISI
r C 51
rl fM rf> •*
W Ch ffi o
.=•.=• i- 7<
^" i5T en rM
s s s s
Dale
Figure 75. Time series of predicted and observed surface water salinities
for the whole Bay
30
29
^
(/) TO
o. 2a
I 27
1 26
T3
d) 25
ID
&_
Q.
23
22
22
Overall Surface Salinity
y = 0.8991x +1.0088
R2 = 0.7955
23 24 25 26 27 28 29 30
Observed Salinity (psu)
Figure 76. Scatter plot of predicted and observed surface water salinities
for the whole Bay
93
-------
3S
w
£r JO
*
Overall Bottom Salinity
Pretl Hul Sal Obi toil Sal — Diff (Obs-Pred)
S
Dale
Figure 77. Time series of predicted and observed bottom water salinities for
the whole Bay
31
3-30
Overall Bottom Salinity
= 1.8064x-26.098
Rz = 0.7048
25
27 28 29
Observed Salinity (psu)
30
31
Figure 78. Scatter plot of predicted and observed bottom water salinities for
the whole Bay
94
-------
Overall Observed Salinity
QbiBoiSal Obs Dilf (Bol Surf)
Date
Figure 79. Time series of observed surface and bottom salinities to illustrate salinity
stratification for the whole Bay
Overall Predicted Salinity
if Sjl Prr.l But s-il Pnnd Difl (got-Surf)
£30
J - •
10
'
V
;-ir,
Figure 80. Time series of predicted surface and bottom salinities to illustrate salinity
stratification for the whole Bay
-------
-------
4. SAMPLE RESULTS
In addition to the time series results for water surface elevation, temperature, and salinity that are
presented above for calibration and validation at station locations (Fig. 19), this section provides
more results about the spatial distribution (horizontal and vertical) of the modeled parameters.
Flushing behavior of NB is also presented using the buildup and flushing of a surrogate
conservative constituent (dye). The model generated data are very extensive in both space and
time. Thus, only samples of the spatial distributions are presented during the last day of the
simulation (day 365, December 31, 2009). It is recommended to concentrate on the horizontal
contour plots within the main body of water in the Bay, not the narrow branches and sharp
corners where the interpolation scheme for contouring in SMS may introduce anomalous values.
4.1. Water surface elevation distribution
Appendix C presents examples of the distribution of water surface elevation at 2 h increments
during a full tidal cycle on December 31, 2009. The water surface elevations are directly affected
by the tidal phase at the open boundary and by the spring-neap cycles. The change in slope of the
water surface elevations directly impacts the velocity field. As indicated by the legend of each
distribution, the difference in water surface elevation (i.e., slope) within the Bay is < 20 cm.
4.2. Velocity distribution
Appendix D presents distributions of surface and bottom velocities every 2 h during a full tidal
cycle on December 31, 2009. The estuarine type flow is evident in various areas of the Bay
where the surface and bottom velocities have different directions. This type of estuarine flow
exists during various phases of the tidal cycle. Magnitudes of the velocity vectors are directly
affected by the tidal phase and by the spring-neap cycles.
4.3. Temperature distribution
Appendix E presents bimonthly distributions of surface and bottom temperatures during 2009.
Due to the cold (below zero) air temperatures during the winter (Fig. 8), which directly impacts
temperatures in the surface water, relatively warmer bottom water (color coded orange) may
manifest in the shallower parts of the Bay in the winter.
Examples of temperature stratification periods during 2009 are presented in Appendix G for all
stations. The degree of temperature stratification is illustrated in Appendix J which presents
monthly temperature profiles throughout 2009 at all stations.
4.4. Salinity distribution
Appendix F presents distributions of surface and bottom salinities every 2 h during a full tidal
cycle on December 31, 2009. The general distributions of surface and bottom salinities indicate
that, in general, fresher surface water extend further south, while the saltier bottom waters (from
RI Sound) extend further northward from the open boundary. This distribution emphasizes the
estuarine behavior and circulation.
97
-------
Examples of salinity stratification periods during 2009 are presented in Appendix H for all
stations. The degree of salinity stratification is illustrated by the monthly salinity profiles at all
stations as presented in Appendix J.
4.5. Water Density and stratification strength
The combined effect of temperature and salinity on the overall vertical density gradient is
presented in Appendix I for the whole year of 2009. The overall density gradient is calculated at
all stations as [surface density-bottom density]/ water depth (kg m"4). The vertical distribution of
water density and stratification strength varies throughout the year (Appendix J). The
stratification strength is calculated from the Brunt-Vaisala frequency as presented in Appendix J.
4.6. Dye distribution
In the initial model setting and draft report, dye was used as a surrogate to study the flushing of
material out of the system and the building of material inside the system. That work is restated in
this report to qualitatively illustrate the same concepts about flushing time. The case of dye
flushing is presented in Appendix K and L. The scenario for dye building is presented in
Appendix M and it was used to approximate initial distributions of salinity and temperature for
various model runs (e.g., Fig. 21). Horizontal diffusivity was turned off (i.e., Ao = c = 0.0) in all
dye experiments.
4.6.1. Dye building-up experiment
A numerical experiment of the spatial behavior of mass build-up within NB was conducted by
applying a surrogate dye with unit concentration at all locations of river inflows. The dye
concentrations in the freshwater inflows from all rivers were set to 1.0 kg m"3. Dye
concentrations at all nodes on the seaward open boundary were set to zero. The numerical
simulation was run and the dye concentration increased due to its loading and flushing through
the open boundary. The simulation was run for 3 months. The dye concentration reached a quasi-
steady state distribution at all the horizontal and vertical locations at different times. Figures in
Appendix M illustrate the local behavior of concentration build up at PR, CP, GB, TW, and FR1.
A delay of-10 days occurs in dye buildup at GB and TW due to their distance from freshwater
inflows. The Hunt river is the closest to GB, but its mean discharge is very small (Table 1). It is
clear that the time period to reach quasi-steady state behavior depends not only on the horizontal
location but also on the vertical location in the water column. Nonetheless, the dye build-up
distribution after at the end of the 3-month simulation was used to define initial conditions for
the transport of mass and properties within NB (i.e., temperature, salinity, and dye).
4.6.2. Dye flushing experiment
Flushing behavior changes spatially within the Bay. Abdelrhman and Cicchetti (2011) presented
a local time scale as the "local flushing time" (LFT), which tracks the replacement of water
within the Bay (local water) with water from outside the Bay (outside water). It is the time period
for the concentration of a uniformly distributed tracer to drop below a threshold value at a
98
-------
specific location within the embayment. Examples of the horizontal changes in LFT are
presented in the above reference. Examples of the vertical profiles of the LFT are presented in
appendix K.
A numerical experiment of the spatial behavior of flushing was conducted by applying a
surrogate dye with unit concentration at all grid locations within NB. Three flushing experiments
were examined; the first was for flushing of estuarine water (i.e., the mixture of fresh and salt
water), the second was for flushing of freshwater only, and the third was for flushing of saltwater
only. It is worth mentioning that flushing times depend on the freshwater input rate and on the
tidal amplitude. The times listed in Table 8 will therefore slightly change based on the onset of
the calculation during the year and the spring-neap cycle. This will not be so important for the
longer flushing times since they will average over the spring-neap cycle.
4.6.2.1. Flushing of estuarine water
Initial dye concentrations were set to 1.0 kg m"3 at all horizontal and vertical locations within the
Bay. Dye concentrations in the freshwater inflows from all rivers were set to zero. Also, dye
concentrations at all nodes on the seaward open boundary were set to zero. The numerical
simulation was run and the dye concentration dropped due to its flushing through the seaward
open boundary and its dilution with the inflow of freshwater from the rivers. The time instant
when concentration reached its e-folding value (1/e = 37%) of its initial value defined the
flushing time at a specific location. Figures in Appendix K illustrate the local flushing time at
PR, CP, GB, TW, and FR1. It is clear that the local flushing time depends not only on the
horizontal location but also on the vertical location in the water column as shown in Appendix L
and Table 8.
4.6.2.2. Flushing of freshwater
The flushing of freshwater that existed within the Bay was examined. The initial distribution of
freshwater concentration was set to the values calculated by the dye build-up experiment. To
insure that only the flushing of freshwater that existed within the Bay was examined; dye
concentrations in the freshwater inflows from all rivers were set to zero. Concentrations on the
seaward open boundary were also set to zero. The numerical simulation was run and the dye
mass dropped due to its flushing through the open boundary and its dilution with the inflow of
external freshwater from rivers. The time instant when concentration reached its e-folding value
(1/e = 37%) of its initial value defined the freshwater flushing time at a specific location. Figures
in Appendix K illustrate the local freshwater flushing time at PR, CP, GB, TW, and FR. The dye
concentration did not exhibit the monotonically decreasing (exponential) behavior with time due
to the initial spatial distribution (monotonically decreasing from north to south), which
sometimes caused water with higher concentrations to manifest at downstream locations with
lower concentrations. It is clear that the local freshwater flushing time depends not only on the
horizontal location but also on the vertical location in the water column as shown in Appendix L
and Table 8.
99
-------
4.6.2.3. Flushing of saltwater
The flushing of saltwater that existed within the Bay was examined. The initial distribution of
saltwater concentration was set to the complements of the values calculated by the dye build-up
experiment. To insure that only the flushing of saltwater that existed within the Bay was
examined; dye concentrations in the freshwater inflows from all rivers were set to zero.
Concentrations on the seaward open boundary were also set to zero. The numerical simulation
was run and the dye mass dropped due to its flushing through the open boundary and its dilution
with the inflow of external freshwater from rivers and saltwater from the open boundary. The
time instant when concentration reached its e-folding value (1/e = 37%) of its initial value
defined the saltwater flushing time at a specific location. Figures in Appendix K illustrate the
local saltwater flushing time at PR, CP, GB, TW, and FR1. The dye concentration did not exhibit
the monotonically decreasing (exponential) behavior with time due to the initial spatial
distribution and tidal excursion, which sometimes caused higher concentrations to manifest at
locations with lower concentrations in the upstream. It is clear that the local saltwater flushing
time depends not only on the horizontal location but also on the vertical location in the water
column as shown in Appendix L and Table 8.
Table 8. Local flushing time at stations PR, CP, GB, TW, and FR1
Location
PR
CP
GB
TW
FR1
Estuarine
flushing time
of bottom
layer (d)
34.54
38.88
59.79
17.42
37.38
Estuarine
flushing time
of surface
layer (d)
4.21
38.46
63.17
35.63
37.38
Freshwater
flushing time
of bottom
layer (d)
10.29
17.33
58.83
17.34
13.08
Freshwater
flushing time
of surface
layer (d)
4.21
17.42
63.33
40.42
13.08
Saltwater
flushing time
of bottom
layer (d)
54.17
42
54.63
15.29
49.79
Saltwater
flushing time
of surface
layer (d)
42.04
41.04
57.7
26.71
49.17
100
-------
5. SUMMARY AND CONCLUSION
The EFDC was used to define the 3-D circulation and transport behavior in NB. The same coarse
grid was used to predict both hydrodynamics (water surface elevation and velocity) and transport
(temperature and salinity) without compromising model predictions. This report presented
the hydrodynamic model together with a study of the horizontal and vertical behavior of flushing
in the Bay. Available data were used to force, calibrate, and validate the model. The
results indicated that the predicted hydrodynamics were compatible with observations. The
model showed that flushing time to the ocean boundary depends on the horizontal and vertical
location within the water body and that it can reach 60 days in parts of the water column within
some regions in NB (e.g., Greenwich Bay). The generated hydrodynamics could be used by other
models (e.g., WASP7) to study ecological and water quality dynamics within the Bay.
5.1. More on Calibration
As recommended and stated by the developer (Tetra Tech, 2005) "The horizontal mass diffusion
terms are generally omitted in the numerical solution when the model is configured for three-
dimensional simulation.'" This recommendation was implemented in the first draft of this report
by setting the two calibration coefficients in the Smagorinsky's scheme c and A0 to zero, which
turns-off the horizontal turbulent diffusion for both hydrodynamics and transport. During the
calibration process it was noticed that excellent agreements existed for the hydrodynamics (water
surface elevation and velocity) when the horizontal momentum diffusion (i.e., eddy viscosity)
was zero, but the transport of mass and property (i.e., salinity and temperature, respectively) was
not as good. Better agreements for transport were achieved when mass diffusion was activated
by setting c = 0.005 and A0 = 5.0. In addition to the possibility of model instability, higher values
of these calibration coefficients can affect water surface elevation and velocity structures, which
can severely compromise vertical stratification in the water column.
Predicted salinities showed values lower than observed at some locations, especially near the
bottom. Similar difficulties in matching salinity predictions to observations were reported by
Applied Science Associates (ASA, 2005) after applying their 3-D model (BFHYDRO) with fine
segments (size 100-400 m) and eleven sigma layers. Their predictions lacked stratification in the
water column. Although the reported uncertainty range in the USGS stream flows is 6%-19%,
(Harmel et al. 2006), which may have impact on salinity predictions, all stream flows were kept
at the values reported by USGS in this work. Further investigation may be necessary to address
this issue.
Personal communications were initiated with the developer (John Hamrick, Tetra Tech, Inc) to
separate the turbulent diffusion of momentum from that for mass. Separating the two diffusion
processes was justified by the developer "There can be a situation where different horizontal
diffusion is needed. The advection scheme in the momentum equation solver is first order upwind
which has inherent numerical diffusion and typically does not need addition constant or variable
diffusion. The advection scheme for salinity is high order upwind bias and has very little
101
-------
numerical diffusion and enhancing horizontal diffusion may be necessary. The reason lower
order upwind is not used for salt is that it tends to over mix in vertical reducing stratification''
This recommendation was achieved later at AED in a new version of EFDC which allowed the
model to use calibration coefficients for momentum diffusion only, mass diffusion only, or both.
In addition, it was noticed that the transport subroutines which calculate mass diffusivity were
completely bypassed in the EFDC code. This deficiency was corrected to include horizontal
diffusion in the transport of temperature and salinity.
The EFDC executable used in this report (efdcl_mobile.exe) has many accommodations to
communicate hydrodynamics, temperature, and salinity to WASP7. To keep such
accommodations and the above-mentioned new modifications, a new executable has to be
generated to communicate data to WASP7.
The impact of defining accurate concentration time series at the open boundary (e.g., for salinity)
should not be underestimated. Such concentrations would reflect the values in the ebbing or
flooding water across the boundary. However, based on the location of the open boundary, some
of the flooding flow would have concentrations related to the preceding ebbing flow and actual
entry of the assigned boundary values would be delayed (Section 2.2.2).
5.2. Hydrodynamics and Water Quality
5.2.1. Water quality with WASP7
The same coarse horizontal and vertical segments will be used by the water quality model
(WASP7). The EFDC provides the hydrodynamics in a file with extension type "hyd" which
includes all the following hydrodynamic data:
1. Volume of segment
2. Depth of segment
3. Velocity
4. Flow along each flow path between segments
5. Dispersion along each flow path between segments
The water quality time step can be 2-3 times the hydrodynamic time step (120 s). At the water
quality time step, hydrodynamic data can be extracted from the "hyd" file and used to transport
other water quality constituents between the various horizontal and vertical segments. The
temporal resolution (few minutes) permits the study of hourly behavior of water quality
constituents.
5.2.2. Water quality with EFDC
The EFDC has built in water quality capacities that can model 22 state variables. This capability
was tested for NB and the model simulated the month of January, 2009 (results not presented).
The full implementation of this capability requires proper definition of water quality model
102
-------
parameters, initial and boundary conditions, and loading of the following state variables as
presented in Abdelrhman (2015):
1. cyanobacteria
2. diatom algae
3. green algae
4. refractory particulate organic carbon
5. labile particulate organic carbon
6. dissolved carbon
7. refractory part, organic phosphorus
8. labile particulate organic phosphorus
9. dissolved organic phosphorus
10. total phosphate
11. refractory part, organic nitrogen
12. labile part, organic nitrogen
13. dissolved organic nitrogen
14. ammonia nitrogen
15. nitrate nitrogen
16. particulate biogenic silica
17. dissolved available silica
18. chemical oxygen demand
19. dissolved oxygen
20. total active metal
21. fecal coliform bacteria
22. macro algae
5.3. Future scenarios
Predictions of hydrodynamics for future scenarios require some estimates of future forcing
functions. Climate change is the major factor affecting most of the forcing functions. The
increased atmospheric temperature can cause changes in snow melt, freshwater flow, sea level,
precipitation, evaporation, water salinity, and wind patterns. The following list provides
suggestions to estimate forcing functions for future scenarios.
1. Sea level rise has to be estimated. The MSL has to be adjusted accordingly. This
adjustment will impact the water depth (bathymetry) in the Bay. Assuming that the sea
level rise will continue at the same global rate of 3.3±0.4 mm per year (Nicholls and
Cazenave 2010) (see also http://en.wikipedia.org/wiki/Current sea level rise), after N
years the increase in MSL will be ~ 3.3 N mm.
103
-------
2. Tide elevations, at the seaward open boundary, have to be calculated from the major
astronomical tidal constituents and be referenced to the new MSL. The six major tidal
constituents are S2, M2, N2, Kl, PI, and Ol at Newport, RI (Table 9).
(http://tidesandcurrents.noaa.gov/data_menu.shtml?stn=8452660 Newport,
RI&type=Harmonic Constituents)
3. The present typical values of freshwater inflows from the eight rivers have to be
calculated for typical wet/dry years. The time series of the typical yearly flows have to
be adjusted (multiplied by a factor) to account for future wet/dry weather. Changes in
extreme values may be inferred from statistical analysis of present and expected future
distributions.
4. Salinity concentration at the open boundaries can be similar to current conditions:
32-35 ppt at the seaward open boundary and zero for all river inflows
5. Global warming will increase water temperature at the seaward open boundary
(assumed the same as at NP) and the temperatures of all river inflows (assumed similar to
air temperature) throughout the year. The increase in air temperature can be chosen from
the published values (1.1-6.4°C during the 21st century
(http://en.wikipedia.org/wiki/Global Warming)) and added to the values approximated
by the above-mentioned sin-wave (e.g., Figs. 8 and 81). (The estimated values for the
coefficients for the sinusoidal function for sea water are Tmax = 22°C, Tmean = 12.3°C, and
O= 1.64 radian.)
6. Air temperature can be estimated by adding the increase in temperature (see No. 5 above)
to the mean air temperature (Tmean) and use the previously mentioned sinusoidal variation
of air temperature (Fig. 7).
7. Relative humidity (from climate models under the impact of future climate changes).
8. Atmospheric pressure has to be estimated (from climate models).
9. Rain fall rate (from climate models)
10. Evaporation rate will be calculated internally by the model
11. Shortwave solar radiation (Wh m"2) (from climate models)
12. Cloud cover (%) (from climate models)
13. Wind speed (m s"1) (from climate models)
14. Wind direction from North (deg) (from climate models)
104
-------
Table 9. Tidal constituents at Newport
Tidal Constituent
Name
S2-Principal solar semidiurnal
M2-Principal lunar semidiurnal
N2-Larger lunar elliptic semidiurnal
Kl -Lunar diurnal
Pi-Solar diurnal
Ol -Lunar diurnal
Amplitude
(m)
0.110
0.518
0.123
0.065
0.023
0.052
Phase
(Deg)
24.3
2.20
346.2
166.6
181.7
200.4
Period
(h)
12.0000
12.4206
12.6583
23.9345
24.0659
25.8193
Speed
(Deg/h)
30.0000000
28.9841042
28.4397295
15.0410686
14.9589314
13.9430356
2003 temperature at Providence, Quonset, New Port, Conimicut
28
26
24
22
20
P 18
Tie
3 14
gl2
|- 10
o 8
6
4
2
0
-2
Avg
PVD
Quonset
NP
Conimicut
M J J
Time (month)
Figure 81. Observed water temperatures at PR, QP1, CP, and NP during 2003
(at respective depths of 0.55 m, 6.9 m, 1.6 m, and 1.3 m from MLLW)
105
-------
-------
Appendix A. Definition of Skill Parameters
A review of skill assessment for coupled biological/physical models of marine systems was
presented by Stow et al. (2009). The assessment analysis compares a model prediction, which
has an unknown "prediction error" from the true value (or truth), to observation, which has
unknown "observation error" from the truth. The difference between the observation and the
prediction is known and is identified as the "misfit" or "residual". The skill analysis examines
the closeness of predictions to "truth" by quantifying the misfits which should be small and
noisy. It is worth mentioning that a slight phase error (e.g., in time) between a prediction and an
observation can result in a big misfit. The following six statistical parameters are quoted from
Stow et al. (2003) to guide the model skill analysis presented here.
The correlation coefficient, r, measures the tendency of the predicted (p) and observed (o)
values to vary together. Its values range from -1 to +1. Ideal value is one.
r =
Where n is the number of records and the over-bar indicates their average value.
The root mean square error, RMSE, measures the size of the discrepancies between predicted
and observed values. Values close to zero are ideal.
RMSE=
n
The reliability index, RI, quantifies the average factor by which model predictions differ from
observation. It should be close to one for ideal predictions.
IV1" / oi
-y (log-
nZ-i^V apt.
107
-------
The average error (bias), AE, measures the size of the discrepancies between predicted and
observed values. Values close to zero are ideal. However, this parameter can be misleading
because positive and negative discrepancies can be large, but they may cancel each other.
,„ Z"=ifo ~ oi) _ _
AE = - = p - o
n
The average absolute error, AAE, measures the size of the discrepancies between predicted
and observed values. Values close to zero are ideal.
n
The modeling efficiency, MEF, measures how close predictions are relative to the average of
the observations. Values close to one are ideal.
The absolute relative error, ARE, quantifies the difference between predictions and
observations (Kim et al., 2010). Small values indicate better predictions.
The index of agreement, d, reflects the degree to which the observed value is accurately
estimated by the predicted value. It ranges between zero and one. Values close to one are ideal.
108
-------
Appendix B. Input files
The EFDC input files used in this work included the following files:
aser.inp Atmospheric forcing time series for air pressure, temperature, relative humidity,
rain, evaporation, short wave solar radiation, and cloud cover (assumed uniform)
cell.inp Horizontal cell type identifier
celllt.inp Horizontal cell type identifier for saving mean mass transport.
dser.inp Dye concentration time series
dxdy.inp Horizontal grid spacing, depth, bottom elevation, bottom roughness, and
vegetation class
dye.inp Initial dye distribution
efdc.inp Master input file
Ixly.inp Horizontal cell center coordinates and cell orientations
pser.inp Open boundary water surface elevation time series
qser.inp Volumetric source-sink time series
salt.inp Salinity initial distribution
show.inp Screen print of conditions in a specified cell (optional)
sser.inp Salinity time series file.
temp.inp Temperature initial distribution
tser.inp Temperature time series file
wndmap.inp Wind forcing map of weighting factors (assumed uniform)
wser.inp Wind forcing time series
Templates of most input files are provided and documented in the EFDC manuals
(Tetra Tech. 2002).
109
-------
Appendix C. Elevation Contours
The predicted horizontal distributions of water surface elevation during a full tidal cycle at hours
2, 4, 6, 8, 10 and 12 on December 31, 2009, as marked on the tidal phase graph.
PR Observed Water Surface Elevation
24 26
Time (h) on 12/31/2009
110
-------
Ill
-------
Appendix D. Horizontal Velocity Vectors
The predicted horizontal distributions of surface and bottom velocity vectors during a full tidal
cycle at hours 2, 4, 6, 8, 10, 12, and 14 on December 31, 2009, as marked on the tidal phase
graph in Appendix C. Areas showing estuarine circulation (with opposite surface and bottom
flows) can be identified together with their tidal phase.
(Kttn* -
li
u.
ft
112
-------
faaKr Module 601 Ml *'h
OS7i«i *
0,00 nA
Ke*B4r M44U4 Sun vn jl f,h
flJiim fc
60IA4 •
VI TV I
NT •'•'•-
W
.y*
I
J
^"^- '^/P
••--•;;>v :"W^'
ifJI
/v
//
lv
'^ . ....
s/ss-r-irtrfti.
f^^r:
SMVs.t*.^^
ifflp?*?
j.|niniiriiuiu
113
-------
7
ft
/r,
%
•v-*-V
'«///
'-••> '< svyA*!
;0f \ tm-vV^
ih^Sgr i
_ j; c-v,, c-4 r-fh
1,1'111 J-jjl ',U? M'J
«y §
i'/// w,m |li
3-:ar«r Motlufe Surt Vel »: 10n
OoOitift —
C
"i^
ril
1
" ' •/
f
1
Scartttf WcC'Je Dot Vfti «l10h t
074 HA » \
uoom*— I
"\
'C-
-".u
.
* I
1
1
i
114
-------
Stiller Module Sot V»l at Mr,
ftilte ifeMt Surf \W it ISO
3cjnet Module Ooi Vci at
ft.DO aw
sfluteOolVelatur. I v
115
-------
Appendix E. Temperature Contours
The predicted bimonthly horizontal distributions of surface and bottom temperatures on calendar
days 60, 120, 180, 240, 300, 365 in 2009.
116
-------
117
-------
118
-------
Appendix F. Salinity Contours
The predicted horizontal distributions of surface and bottom salinity during a full tidal cycle at
hours 2, 4, 6, 8, 10 and 12 on December 31, 2009, as marked on the tidal phase graph in
Appendix C.
119
-------
Surf Sal J1«11 Ipptl
120
-------
121
-------
Appendix G. Temperature Stratification
This appendix presents time series graphs of water temperature in the surface and bottom layers
of the water column at all stations (Table 4). The surface layer represents the upper 1/8 of the
water column and the bottom layer represents the lower 1/8 of the water column. The duration
and degree of temperature stratification varies based on the station location. Temperature
stratification is influenced by (among other factors) total water depth, proximity to freshwater
inflows, seawater intrusion, vertical mixing, and heat exchange at the surface and bottom
boundaries.
PD Temperature
Prcd Surf Temp Prccl BdlTemp
8888888
^ * • - "
Date
PR Temperature
•I'ted butf temp l'r«d Hot I'ernp
122
-------
BR Temperature
-Prpd surf Temp Pred Uor Temp
<<<
r< (N
Date
CP Temperature
•Pred Surf Temp —
SS8S8SS8
kO 00
< <
rj m
i
;
"
j j
Oi
'ft -'
Date
123
-------
GB Temperature
•Pred Surf Temp PrcdBottemp
SR Bottom Temperature
-PredS«jrf Tc-mp
Prod Bot Temp
Date
124
-------
NPI Temperature
• Pred Surf Iprnp Prpd Uottpmp
MV Temperature
•Pred Surf remp
Date
125
-------
PP Temperature
-Prprt surf Temp Pred Uot Temp
$ s $
DJtc
QP1 Temperature
• Prp d su rf Tpmp Prpd Hoi t cm p
Date
126
-------
TW Temperature
•Pned Surf Temp Pred Bot Temp
Date
GD Temperature
-PredSurfTfimp
• Pi«d &ul Temp
127
-------
MH Temperature
Pred Bnl letup
Date
FR1 Temperature
• Pr*>d
-------
Appendix H. Salinity Stratification
This appendix presents time series graphs of water salinity in the surface and bottom layers of
the water column at all stations (Table 4). The surface layer represents the upper 1/8 of the water
column and the bottom layer represents the lower 1/8 of the water column. The salinity
stratification exists at all locations at various degrees. The duration and degree of salinity
stratification varies based on the station location. Salinity stratification is influenced by (among
other factors) total water depth, proximity to freshwater inflows, and seawater intrusion.
PR Salinity
— Pred surf sal - Prtd Bet Sol
PD Salinity
Pied Suit Sal Pied But Sal
129
-------
BR Salinity
-Pred Surf Sal
DM
CP Salinity
Pred Surf Sal
88S8888S8
" " "* ^ "
CO Of
130
-------
GB Salinity
-Pred Surf Sal —PredBotSal
Date
SR Salinity
PredSurtSa!
88S8888S8
" " "* ^ "
CO Of
131
-------
NPI Salinity
I'rud S4irf Sdl Pred Bt,l Sdl
Date
«
28
-26
&2J
12
MV Salinity
Pred Sort Sal - Pred Bot Sal
1G
14
12
1
-.
r
*i
'
|S8SSSSSS8|ggJ
^(^U3tn^ ^ m m r
J<«SSg^SgsS --- 3. < "
^f^^iNfOTUTi^p^OOOl O T^ r*f r
Date
132
-------
PP Salinity
-,.'
m
28
--- ?,,
£M
.£«
£ 20
™ IB
16
14
I?
••Vwyvvvv-Y
, JI
-------
TW Salinity
Pred Surt Sal - Pred Bot Sal
Date
GD Salinity
i Surf •jl
Pfed Bot Sal
Date
134
-------
MH Salinity
-Pred Surf Sal
8SSSS£8££
r* u> DO h- r- «j r> u-. *r
S^^dCLC^^.^^
rt(Nf>Tin>flt^lBO>
Date
4r TO ff M
*-< «H «-( ^
^ > ?? >
r< rt M
FRl Salinity
Pred Surf Sal
Pred Bot sal
8 8 8 8 S §
--
135
-------
Appendix I. Vertical Density Gradient
Vertical density gradient at any station (Table 4) is calculated as ——-, where p; is the potential
H
density (kg m"3) of water in layer i as a function of water temperature and salinity (defined by the
equation of state) and H is the total water depth (m) at the station.
PD Vertical Density Gradient
PR Vertical Density Gradient
&
t3
01
D
"ro 3
u -3
t
m
g
O O O 0 O S 0
Date
I
a
S
136
-------
BR Vertical Density Gradient
CP Vertical Density Gradient
GB Vertical Density Gradient
137
-------
c
Oj
ra
5
i:
O
"S
-i
-2
-3
SR Vertical Density Gradient
Ffff^
:•:
—i
••N
r^-iDoor~~r^-LDiDLn^-
o? rC
888888
8 8
to" ^
01
I
01
1
S
rj
-------
PP Vertical Density Gradient
Date
t
M
r
.3!
TD
E
0)
Q
~
o -3
QP1 Vertical Density Gradient
~
I
i
3J
1
-1
-2
-3
il)
X
I
Date
TW Vertical Density Gradient
g g g g g g
Date
S
a
8
^
s
139
-------
GD Vertical Density Gradient
t
E
$
c
.5!
TO
to
>• -2
-3
e
O)
:c
o
—I
^1
•3-.
o
01
o
01
o
Date
-1
0>
1
01
Q
-3
V
\fi»
GD
H
Veri
JWiW
tical
M^
Density Gradient
t
0)
Date
s
>-l
H
^H"
MH Vertical Density Gradient
140
-------
FR1 Vertical Density Gradient
141
-------
Appendix J. Profiles of Temperature, Salinity, Density, and Stratification
Strength
The actual stratification of the water column is based on the water density, which encompasses
the effect of both temperature and salinity. Profiles of temperature, salinity, density, and the
square of the Brunt- Vaisala frequency at the end of each month are presented in this appendix to
show the degree and seasonality of stratification strength at all stations (Table 4). The Brunt-
Vaisala frequency is given by:
N= (rad. s-)
Where g is the gravitational acceleration (9.8 m s"2), p is the potential density as a function of
temperature and salinity as calculated from the equation of state, and (dpz/dz) is the gradient of
potential density at height z ( -^ = Pl+1 Pl, where i =1 to 7 is the sigma layer ID and Az is the
uZ l\Z
layer thickness = H/8 (Table 4). Frequencies are indicated at the top of layer i.
142
-------
PD Temperature Profiles
— —-jji -.—-fi-l, _ —MJI ftp M
e
j
tl
i;
;
i
> i
K
11
11
11
\
I
1
1
•
/
\
1
J
~f
1
4 6 ! 10 11 14 1« IS 10 11 24 16 38 30
PD Temperature Profiles
- Jul
-Ort ---l*jv ---
d
a s
1
1
1
1
1
1
1
[
1
'
i
i
X
\
1
1
1
"" y
\
\
c
j
/
0 2 < 6 & 10 13 1*
jo
PD Salinity Profiles
---Jar, H* M.JI Apr May Mi
PD Salinity Profiles
-Ju1 Aug 5ep Oci ---
0 I t t, 8 10 1J H 1C, IS
-I r * HP 11 M 16 IS ?D JI j.i ?« J* 30 91
PD Density Profiles
PD Density Profiles
lar, *
- MJJ ^^— Api ' Miiy ^^^ *U'n
-;ul
ALJJ - Svp
—— — May --—On
1.000
1,010
IjOJO
l.COO
J.OOO
l.PID l.OJP
D*nil(y(k|fn'}
;,i i;..
PD Stratification Strength
IVI — — — Ffh ^^^ Wjf Apr
PD Stratification Strength
-lul Ang 1jp fVt — — — ^inlJ
i-
1
It
4
I
1
1
I
UUI IUU DM) U.O4
N>(mLf^
Ulfc 1Mb
«.«,
143
-------
PR Temperature Profiles
— - — Ian — — -irt>
0 1 1 6 8 10 12 11 16 18 10 22 1A 16 28 30
T«
-------
BR Temperature Profiles
— •— — Ian — — — 11*> • War ^-^ Apr ^^^ May
2 4 6 8 10 12 14 IS 18 20 22 24 26 28 }C
BR Temperatura Profiles
-Dct Nov >c
_
2 I
i
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 JO
BR Salinity Profiles
BR Salinity Profiles
- May lun
i --- NW --- Oft
K
D 2 4 6 8 10 12 11 16 18 20 22 24 26 2g 30 32 34
Salinity ii'F-t;
0 2 4 6 8 10 12 11 16 IS ID 22 21 26 28 W 32 34
MM* ft*)
BR Density Profiles
BR Density Profiles
un i
-May tun
•M
5*p - Otl --- (*« --- ^t
1.010
UOBID
1.030
1,000
1.010 ifllO
Chiriltv I k« n ;
1,030
BR Stratification Strength
BR Stratification Strength
— - — Ian — — - icb Mir
-I.,I
Qct — — — ISov - —— 0«
002 0.03
N<{nd.r>>>
0.04 0.05
0.04 0.05
145
-------
CP Temperature Profiles
___Un ——— hch Mjt
CP Temperature Profiles
— — — Nrrv
CP Salinity Profiles
CP Salinity Profiles
-- -Jan --- tO>
a
Mil - Apt
0 1
» 1n 1.' 14 Ib 18 20
,« J^ -M
CP Density Profiles
CP Density Profiles
-lul
- O( I — — — Nov - - - DXC
tooo
1,030
n>')
CP Stratification Strength
CP Stratification Strength
May - — lun
-Oet —-—(
l)«l U1IJ
00!
146
-------
- - - ',1
8 \ 1
711
t: \
f, i
i \
>
1 T ' T
1 * *
0 2
8
ll
*:
2
1
Ian
8
i.:
f.
i
1.000
---(..1
l^
* %
! -*3
D
GB Temperature Profiles
— — — ttb ^^^ Wir ^—— Apr ^^^ Mj-y ^—— Ign
4 6 8 ID 12 Id
T*mpnr - — Var A^ir Mjr - — lur
/ll '
\\ '
1 ft i
\i\i
U i
\C I
i.oio ijoao i.ojo
GB Stratification Strength
.**
^
£^~^^_
0.01 0.02 1.1.0.' a(H 0.05
A
7
it
^
1
I
•(
7
ll
i;
2
1
(
1
ll
i;
u
i -
I
GB Temperature Profiles
— lul /Ujf - — icp Oct ———Not — — -0«
; ;
i 1
L !
I '
l' / ^
J / 4 b H 11] M 14 lt> 1.1 /IP 1! ft ,':. .'» «:
GB Salinity Profiles
_ |ul AUH S»p Oet Ptov 0«c
) 2 4 ft 8 10 12 14 1C IK 20
SalMty (Wit
•
1 s
r, A . __
22 34 » J8 30 33 34
GB Density Profiles
— ml AM,: - — Sep Oe! ---rwtv Dec
1
i
\ \ UA
DU Iflll) I.IUO l! 411
GB Stratification Strength
1 0.31 DCS 0.03
HMrad.s'K
0.04 0.05
147
-------
e
*•
2
1 •
1
i
*:
a
i •
*
7
*;
u
i. «
j
i
•
SR Temperature Profiles
— tan ---fH, - — Mi, ,ipi M.I/ - — lun
1 \
/
1 1
1 |
1 1
I • ,' .
» ! 4 6 S 10 1! M lb is 10 13 34 36 38 30
TXnpentur* ff|
SR Salinity Profiles
— Ian — — — F-rfi M.II ^^— Apr May tun
> 7 * 0 8 10 11 14 16 IS K> ?3 .'-S J(. -*8 30 3J 3*
4««nity (PPt)
SR Density Profiles
— Jaf, — — — HA ^^— M.U ^^— flpf M.iy ^^— mn
It
__
00 i.Dio i/wo '.nir*
SR Stratification Strength
— Ian — — — frb Mjr Apr May lun
oc
(5>r
^^p — ^^
g_(n (ijjj D.IU ujrt v.»
SR Temperature Profiles
lu iv.g - — vp t>rl — — — Mow — — — net
* i 1
h ' '
i!
t* ' i i
"* 3 1 1
1
i
/
* 1 1
1 4—i i "
O ? 4 fi R in 1? 14 If. 1R 70 n >4 3d JR Vl
T*mp«f»tur* ( <")
SR Salinity Profiles
M ^u« Ijt^p I-OL! — — — f*w — — — Q*t
f
4
J *
5 4
2
1
0 2 4 b B ID U 14 16 IK 70 H
SiBnity (ppl|
SR Density Profiles
i
r *
s4
-• ^
?
i
J.I ?fa 25 3Q S3 34
— *v — Ike
ijKW vmti l.OTn l.mn
SR Stratification Strength
* T
* a
5 ^ v
)
u^,.r uw ^
148
-------
NP! Temperature Profiles
.trt, Mar Apr Vjy
! i
i »
• i
• i
i i
i i
034601012141419 JO if
Temfiwjtura |t}
NPI Temperature Profiles
-:J| AU|J SfU CHI — —-Nw
i!
1:
i
0 1 1 6 9 10 12 14 16 IX 10 22 U 26 2t 30
T«np.niluf« (T(
NP( Salinity Profiles
— —-Ian ___[.rt MAI Apr May >m
NPI Salinity Profiles
;ul AL^ St-p ot! — — — rtjv D«
i
'• R ID
IS W J3 2*
30
0 7 4 6 » 10 13 14 I* IS 70 JJ JM ?<• » M 3J
friHtgrtpf^l
NPI Density Profiles
NPI Density Profiles
H
7
IjOW
1X"0
WMO
UBO
l.COG
NPI Stratification Strength
— _— ign ———rrt, Mjf tfi Mjy Inn
NPI Stratification Strength
Ocl --- Nov --- oef
o.oi n.o?
0.0* (UK
149
-------
MV Temperature Profiles
T I
i
* * /
z 5 1
S , I /
*! \
' * i
i
U / * b « 1U U 14 111 IS ffi fl i* i<> i!8 *l
TBTOfMHatUlt | TT I
MV Salinity Profiles
— — — Jan — — — Icb - — Mar Apr May Jun
jt
Z 5 M
6 ' \
• J
0 J '4 * * 10 1J 14 Id 18 .'.- .'-' « « JS 30 3/ 34
Still1 Ity (PP-ll
MV Density Profiles
———Jan — •-— Icb - — Ma Api MJV - — tun
; ri:
Vi *l
it N 1
}l
J
1
1,000 1.010 WHO tjOUO
DMHlCy (h x in
MV Stratification Strength
8 • '
PV^
\**w
1
f
a D.OI n.o? osa oivt aor,
NMrad.r"}'
MV Temperature Profiles
1 1 Jl
) \ / \
* • \
3E 5 1
t, 1
* t 1
1 |
1 ' '1
U / * b 8 IV U 14 It. 18 ^U /; f. Vb /B JU
T«fnfi*ratur» I "|i
MV Salinity Profiles
:ul ALg - — Sup at! ---ffav -__o«
1 VI 1
' ?tVw\
x s r^iX
1* V UN
ci 7 i b a w 12 14 16 18 70 » }4 M » 30 33 M
HMjta4
MV Density Profiles
Jul Ai.,q Si'p O^t — — — ttnv —- — OK
* 1 1 III
i ul
V\ 1C
V\ \K
* ;\ v\
L M \\
-j \ a
:
, II 13
1,000 i.ato i.ata i.a»
Dn»ity (k,- i- 1
MV Stratification Strength
0
IC!"^5^^2^^—-^
^^P^
r**>
\*?—
o om n.o? ii.m qn4 a.ar>
H>OwLi^
150
-------
PP Temperature Profiles
(j,,
.III
0 2 -1 6 s to U 14 16 IS ID }J 24 16 28 30
PP Temperature Profiles
-Jill
— — — ftiu — — -= OK
...;
l\
0 2 •! 6 « Id 13 1* 16 IS 20 it It )d 3$ 30
PP Salinity Profiles
PP Salinity Profiles
___jji. ___ft.|, M.I flpi Mjy Ign
X
7
*:
0 ! * 6 s JP u H Jb m » « 24 2« zg 30 37 3*
0 } .1 r « id U H 16 IS 2D 2J 3M
PP Density Profiles
PP Density Profiles
— — — lan — — —
- M Jl ^^— Ap»
-Jut
OK
1,«GO
PP Stratification Strength
. r»' ar ^^— /i.-.f ^— r.i •„•
PP Stratification Strength
aw o.w
o.ot PCM o.»
151
-------
QP1 Temperature Profiles
May
0 2 4 6 S 10 12 14 16 IB 20 22 21 26 28 3D
QP1 Temperatura Profiles
AUH
-Sep
-Oil Ntw Ow.
I
Is
4
3
I
I
I
0 2 4 f> 8 10 12 14 16 IB 10 22 24 26 29 30
T«Rip«wiun {•£)
— — — Jan — — — t
8
QP1 Salinity Profiles
Mar Apr May J jn
QP1 Salinity Profiles
- "'"'P.
ti
it
0 i * 6 8 10 i; U ji- 1* 10 M 24 J6 JS iO JJ
LI .' i I, S IU U 14 Ib IS A>
QP1 Density Profiles
QP1 Density Profiles
— — — J*n — — — I rt
Mar
-Apt
-May
-Jun
-Allg
-Ort _-_
QP1 Stratification Strength
QP1 Stratification Strength
•Sep
n.oi o/» om 0.04
152
-------
TW Temperature Profiles
Jill Fib
\ i /
/ ) J /
4 6 8 10 12 M 16 18 2O 22 24 26 28 30
Tunpcuwrt i°t(
TW Temperature Profiles
-in
-»uj Vnp Oct — — — Wnv — — — Off
» T
/
i6
: 5
^3
^
1
j
0 2 4 6 S 10 12 H 16 IS 20 12 2* 2& 33 30
TW Salinity Profiles
TW Salinity Profiles
Un ---Ftb MM Ac- M«
•Sen * -Oil -— — f
U 2 * 6 8 10 12 14 16 18 20 22 24 26 20 3D .U J4
Slinky(PPfl
TW Density Profiles
---!» F»b - —M* AO' M-iy
TW Density Profiles
AUK
i:
1
«-
>
i
1.HOD
i:
I.Din
unu
1JU30
UT Mty (kf m ')
U*r«tv(k«n,')
TW Stratification Strength
-On Irti
-Apr
TW Stratification Strength
Aug - —iep Ort — — — *(n» - - — IStr
n i
OflJ 0«
f;'!--H .')••
»
nni
007
nna nor.
153
-------
— — — Jjn —
i
7 1
6 1
j- .,
f.
4 1
, 1
.
0 ? 4
1*
I
1
024
— — — Jan —
it
}.
1
1
1.C300
GD Temperature Profiles
- — Fifb - — Mar ^— Apr ^^— May - — liin
-
fi S 10 1? 1* 1C IB ?0 ?? 74 ?G ?8 TO
GD Salinity Profiles
!, « 1ii 1.' 14 11, 18 /;) <>J ,M » J8 30 3i 34
•: /IN iiv |PPt)
GD Density Profiles
-— Icb MJJ Apr M.3V - — *wt
1.LJ1U
J«-ll >lry (k|
trao I.!.-K)
GD Stratification Strength
---Jan Feto M-* Apr May - — Jun
I V
1 /
1;
0
o.m n.o?
tffrA
am 0.04 (UK
l U M Id 18 }0 27 . :
S.fciUyipp.1
II III
::
"
'Si!
GD Density Profiles
do i mi i » '
DBKlty (t( m -1}
GD Stratification Strengt
R
2
t
LOW
h
Hot DM
/^
1
N
y
0 0.01 0.02 OJ03 0.04 O.®
154
-------
NP Temperature Profiles
6 t 10 u 1* 16 » » n M n is 30
NP Temperature Profiles
- iul - Aiig
ou — — — 1*11 — — — OK
NP Salinity Profiles
NP Salinity Profiles
- flug
H n i — — — Nm — — -!)«•
*;
!<
0 1 1 t, s 10 14 14 16 1* 10 2! 24 It n t •• •:.• •
NP Density Profiles
NP Density Profiles
— — - ur, — — —F-rt MJJ apt M.iy
2 •
1
5
1,010
1,030
l.tnn
1,000
1,010 1,020
UiMMity (k( m ")
NP Stratification Strength
— — —Ian — — — frti MJT Apr May
*
7
6
NP Stratification Strength
-lul *m - —icp (Jet —— — No'<
0,11* u.M
0 IV
oo? o.ra
N1 (nd. i-'|'
OM aor»
155
-------
MH Temperature Profiles
---Jjn ---K'b
10 13 It 16 IS 10 « J4 26 it 30
MH Temperature Profiles
AI.J
Sfp
t'Jrt — — — Nrw — — —
f.
4
3
I
1
O J 4
I
1
«
1
-r-
|
1
1
1
1 \
»
1
1
|
1
1
s
/
/
—
~
S 10 1? I* If. 1R K) J7
Tcmpwvtvrt (C|
MH Salinity Profiles
MH Salinity Profiles
vp
C'd — — — tun --
I
nj •! i s so ji it ifr :n 3C 11
irfnitvlnrt]
MH Density Profiles
-~JSB —-f«b M-r
MH Density Profiles
Vp Ort — — ttov
- -- l>r
], '-:•
LOW
IJ01W \SiN>
[Xnifty (kg m !|
I.V40
MH Stratification Strength
-Jtil
*ua
MH Stratification Strength
— — — Frti -^^^ Mar ^^— A:pf -^^^ M^v
J,n
am
156
-------
FR1 Temperature Profiles
- M-H Af
0 } 4 6 S 10 17 It 16 1A 30 >J 74 X ?A
FR1 Temperature Profiles
-ttug Sfp
i
•
•
o > 4 (, R in i? 1* ir, is JO }> H ](, }f< vi
FR1 Salinity Profiles
FBI Salinity Profiles
*-—— l*n — — — K
-Od ——— HIM ---Ore
0 1 4 6 8 10 12 U 16 18 ID 2i i-4 26 28 JO 3i 3*
FR1 Density Profiles
FR1 Density Profiles
M*r Apt Miy lun
-Jul
---ftr, -—Us.
1.
-------
Appendix K. Dye Flushing
Flushing behavior changes with horizontal and vertical location. The local flushing time (LFT)
tracks the replacement of local water with outside water. It is the time period for the
concentration of a uniformly distributed tracer to drop below 1/e (~ 37%) of its initial value at
the location. Examples of the calculation of LFT for estuarine water, freshwater (FW), and
saltwater (SW) are presented in this appendix for stations PR, CP, GB, TW, and FR1 (Table 4).
Material presented in this appendix are from the earlier model runs presented in the draft version
of this report before the final calibration. The purpose of this appendix is to demonstrate the
modeling methodology and results.
FLUSHING OF ESTUARINE WATER
PR Estuarine Flushing
• Layerl Layer2
•LayerG— —Layer?
Laye r3
Laye r4 •
•LayerS
0.0
10 20 30 40 50 60
Time (day)
70
80
90
100
158
-------
0.0
CP Estuarine Flushing
• Layerl Layer2 — Layer3 •
•LayerG— —Layer?— —LayerS"
• Laye r4 Laye r5
10 20 30 40 50 60
Time (day)
70
80 90
100
0.0
GB Estuarine Flushing
• Layerl Layer2
•LayerG— —Layer?
Laye r3 Laye r4 Laye r5
10 20 30 40 50 60
Time (day)
70 80
90 100
159
-------
0.0
TW Estuarine Flushing
• Layerl Layer2 — Layer3 •
•LayerG— — Layer7— — Layer8«
• Laye r4 Laye r5
0 10 20 30
40 50 60
Time (day)
70
80 90
100
0.0
FR1 Estuarine Flushing
Laye rl Laye r2 — Laye r3 Laye r4 •
Laye r6 — Laye r7 Laye r8 ^^— 1/e
•Laye r5
10 20 30 40 50 60
Time (day)
70
80 90
100
160
-------
FLUSHING OF FRESHWATER
PR FW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
CPFW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
161
-------
0.4 -i-
E
M.
c
o
*0.2
8
8
OJ
Q
0.0
GBFW Flushing
•Layerl
•Layer4
Layer?
Layer2
LayerS
t. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
TWFW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
162
-------
FRFW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
FLUSHING OF SALTWATER
i.o
2
c
0.8 -•
0.6
PR SW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
Layer6
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
163
-------
1.0 T
CPSW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
GBSW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
164
-------
1.0 T
TWSW Flushing
•Layerl
•Layer4
Layer?
•Layer2
•LayerS
LayerS
LayerS
LayerG
•l/e*init. Cone.
10 20 30 40 50 60
Time (day)
70
80
90
100
1.0 n
FR1 SW Flushing
• Laye rl
•Laye r4
Laye r7
•Laye r2
• Laye r5
Laye r8
•Laye r3
Laye r6
•l/e*init. Cone.
10 20 30 40 50 60 70 80 90 100
165
-------
Appendix L. Flushing Profiles
Flushing behavior changes in the vertical direction. The local flushing time (LFT) tracks the
replacement of local water with outside water. It is the time period for the concentration of a
uniformly distributed tracer to drop below 1/e (~ 37%) of its initial value at the location.
Examples of the vertical profiles of the LFT for estuarine water, freshwater, and saltwater are
presented in this appendix for stations PR, CP, GB, TW, and FR1 (Table 4). Material presented
in this appendix are from the earlier model runs presented in the draft version of this report
before the final calibration. The purpose of this appendix is to demonstrate the modeling
methodology and results.
Estuarine water Flushing Profiles
•PR
•CP
GB
•TW
FR1
10 20 30 40
Flushing Time (day)
50
60
70
166
-------
Freshwater Flushing Profiles
PR CP GB TW FR1
8 -\
0 10
20 30 40
Flushing Time (day)
50 60
70
7 -
6
f 4
3
Saltwater Flushing Profiles
•PR CP GB TW FR1
10 20 30 40 SO
Flushing Time (day)
GO 70
167
-------
Appendix M. Dye Build up
A surrogate dye concentration of 1.0 kg m"3 was loaded with freshwater inflows from all rivers
(Table 1). The spatial distribution reached a quasi-steady state after 3 month. This appendix
presents both the horizontal spatial distribution as well as the temporal behavior at stations PR,
CP, GB, TW, and FR1 (Table 4). Material presented in this appendix are from the earlier model
runs presented in the draft version of this report before the final calibration. The purpose of this
appendix is to demonstrate the modeling methodology and results.
Spatial distribution at the surface and bottom
168
-------
169
-------
The temporal behavior at the station location
o.o
PR Dye Building
•Layerl-
•LayerS
•Layer2
LayerG
LayerS-
•Layer4
Layer?— LayerS
10 20 30 40 50 60
Time (day)
70
80
90
100
1.0 -r
o
2 0.8
c
O
u
a
Q
I
TO
E
0.6
0.4 -••
0.2 -
0.0
CP Dye Building
•Layerl-
•LayerS
•Layer2
LayerG
LayerS-
Layer?
•Layer4
LayerS
10 20 30 40 50 60
Time (day)
70
80
90
100
170
-------
GB Dye Building
•Layerl-
•LayerS
•Layer2
LayerG
LayerS-
Layer7
•Layer4
LayerS
10 20 30 40 50 60
Time (day)
70
80
90
100
TW Dye Building
•Layerl-
•LayerS
•Layer2
LayerG
LayerS-
Layer?
•Layer4
LayerS
10 20 30 40 50 60
Time (day)
70
80
90
100
171
-------
FR1 Dye Building
•Layerl-
•Laye r5
Laye r2
Laye r6
Laye r3 •
Laye r7
•Layer4
Laye r8
10
20
30
40 50 60
Time (day)
70
80
90
100
172
-------
REFERENCES
Abdelrhman MA. 2015. Three-Dimensional Modeling of Water Quality and Ecology in
Narragansett Bay. Draft ORD-008354. US Environmental Protection Agency.
Abdelrhman MA and Cicchetti G. 2011. Relationships between nutrient enrichment and benthic
function: Local effects and spatial patterns. Estuaries and Coasts doi: 10.1007/s 12237-011-9418-
2. New York: Springer-Verlag.
Applied Science Associates, Inc. 2005. Providence River Hydrodynamic and Water Quality
Modeling in Support of a Nutrient TMDL. ASA Report 99-123
Chen CJ, Qi J, Xu Q and Beardsley RC. 2005. The April 2005 Gulf of Maine Regional FVCOM-
SWAVESimulation. Web testbed.sura.org/sites/default/files/FVCOM_SWAVE_GOM_
April2005_report.pdf
Chen CJ, Beardsley RC and Cowles G. 2006. An Unstructured Grid, Finite-Volume Coastal
Ocean Model: FVCOM User Manual. SMAST/UMASSD Technical Report-06-0602, p. 315.
Fan Y and Brown WS. 2006. On the heat budget for Mount Hope Bay. Natural and
anthropogenic influences on the Mount Hope Bay ecosystem, Northeastern Naturalist., 13
(Special Issue) 4:47-70.
Galperin B, Kantha LH, Hassid S and Rosati A. 1988. A quasi-equilibrium turbulent energy
model for geophysical flows. Atmospheric Sciences 45:55-62.
Haidvogel DB, Arango HG, Hedstrom K, Beckmann A, Malanotte-Rizzoli P and Shchepetkin
AF. 2000. Model evaluation experiments in the North Atlantic Basin: Simulations in nonlinear
terrain-following coordinates. Dynamics of Atmospheres and Oceans 32:239-281.
Hamrick JM. 1992. A Three-Dimensional Environmental Fluid Dynamics Computer Code:
Theoretical and Computational Aspects. Special Report 317. College of William and Mary.
Virginia Institute of Marine Science. 63 pp.
Hamrick JM. 1996. User's Manual for the Environmental Fluid Dynamics Computer Code.
Special Report 331. College of William and Mary. Virginia Institute of Marine Science.
Hamrick JM and Wu TS. 1997. "Computational Design and Optimization of the EFDC/HEM3D
Surface Water Hydrodynamic and Eutrophication Models." Next Generation Environmental
Models and Computational Methods. Eds. Delic G and Wheeler MF. Philadelphia, PA: Society
for Industrial and Applied Mathematics. Chapter 16.
Harmel RD, Cooper RJ, Slade RM, Haney RL and Arnold JG. 2006. Cumulative uncertainty in
measured stream flow and water quality data for small watersheds. Transactions of the ASABE
49(3): 689-701.
173
-------
Huang HD, Justic D, Lane RR, Day JW and Cable JE. 2011. Hydrodynamic response of the
Breton Sound estuary to pulsed Mississippi River inputs. Estuarine Coastal Shelf Science
95:216-231.
Kim T, Sheng YP and Park K. 2010. Modeling water quality and hypoxia dynamics in Upper
Charlotte Harbor, Florida, U.S.A. during 2000. Estuarine Coastal Shelf Science 90:250-263.
Kremer JN, Vaudrey IMP, Ullman DS, Bergondo DL, LaSota N, Kincaid C, Codiga DL and
Brush MJ. 2010. Simulating property exchange in estuarine ecosystem models at ecologically
appropriate scales. Ecological Modeling. 221:1080-1088.
Mellor GL. 1991. An equation of state for numerical models of oceans and estuaries.
Atmospheric and Oceanic Technology 8:609-611.
Mellor GL and Yamada T. 1982. Development of a turbulence closure model for geophysical
fluid problems. Reviews of Geophysics 20(4):851-875
Nicholls RJ and Cazenave A. 2010. Sea-level rise and its impact on coastal zones. Science
Masazine 328(5985):1517-1520.
Nowicki BL and Gold AJ. 2008. "Groundwater Nitrogen Transport and Input Along the
Narragansett Bay Coastal Margin." Science for Ecosystem-Based Management. Springer.
pp 67-100.
Pilson MEQ. 1985.On the residence time of water in Narragansett Bay. Estuaries 8(1):2-14.
Pond S and Pickard GL. 1983. Introductory Dynamical Oceanography, 2nd Edition. New York:
Pergamon Press.
Rhode Island Department of Environmental Management. 2005. Total Maximum Daily Load
Analysis for Greenwich Bay Waters, Pathogen/Bacteria Impairments.
Shonting DH and Cook GS. 1970. On the seasonal distribution of temperature and salinity in
Rhode Island Sound. Limnology and Oceanography 15:100-112.
Shchepetkin AF and Me Williams JC. 2005. The regional oceanic modeling system: A split-
explicit, free-surface, topography-following coordinates oceanic model. Ocean Modeling
9:347-404.
Smagorinsky, J. 1963. General circulation experiments with the primitive equations. Part I: The
basic experiment. Monthly Weather Review 91: 99-152.
Stow CA, Jolliff J, McGillicuddy Jr. DJ, Doney SC, Icarus Allen J, Friedrichs MAM, Rose KA
and Wallhead P. 2009. Skill assessment for coupled biological/physical models of marine
systems. Marine Systems 76:4-15.
174
-------
Stow CA, Roessler C, Borsuk ME, Bowen JD and Reckhow KH. 2003. A comparison of
estuarine water quality models for TMDL development in the Neuse River Estuary. Water
Resources Planning and Management 129:307-314.
Tetra Tech., Inc. 2002. User's Manual for Environmental Fluid Dynamics Code, Hydro Version
(EFDC-Hydro), Release 1.00. Atlanta, GA: Draft Report for US Environmental Protection
Agency.
Tetra Tech., Inc. 2005. A Hydrodynamic and Water Quality Model for the Lower Charles River
Basin, Massachusetts. Boston, MA: Draft Report for US Environmental Protection Agency.
Wilcox S. 2012. National Solar Radiation Database 1991-2010 Update: User'sManual.
NREL/TP-5500-54824. Oak Ridge, TN: US Department of Energy.
Wool AT, Davie AR and Rodriguez HN. 2002. Development of three dimensional
hydrodynamic and water quality model to support total maximum daily load decision process for
the Neuse River Estuary, North Carolina. Water Resources Planning and Management
129(4):295-306.
Zhao L, Chen C and Cowles G. 2006. Tidal flushing and eddy shedding in Mount Hope Bay and
Narragansett Bay: An application of FVCOM. Geophysical Research 111(C10):C10015,
doi:10.1029/2005JC003135.
175
------- |