Circulation and Euent Dilution

Modeling in Massachusetts Bay

Model Implementation Verication and

Results

Richard PSignell

USGeological Survey

Woods Hole Massachusetts

Harry LJenter

USGeological Survey

Reston Virginia

Alan FBlumberg

HydroQual Inc

Mahwah New Jersey

January 

USGeological Survey Open File Report 	

This report has not been reviewed for conformity

with USGeological Survey editorial standards
Abstract

A threedimensional hydrodynamic model was developed as part of a coopera

tive US Geological SurveyMassachusetts Water Resources Authority program to

study contaminated sediment accumulation and transport in Massachusetts Bay This report details the development of the model and assesses how well the model

represents observed currents and water properties in the bay It also summarizes

circulation and comparative euent dilution simulations from existing and future

Boston sewage outfalls over a threeyear period from October  	 to December


 

The ECOMsi model a semiimplicit version of the Blumberg and Mellor

	
 Estuarine Coastal and Ocean Model is shown to reproduce many of the

important hydrodynamical features of Massachusetts Bay the seasonal evolution

of the pycnocline the mean ow pattern and the strength of subtidal current

uctuations Throughout the simulation period during both vertically wellmixed

and stratied conditions the seasonal statistics of observed currents are well

represented by the model The model is therefore appropriate for studying the

average dilution of sewage euent and other continuously discharged substances

over seasonal time scales

The ability of the model to reproduce individual ow events varies with season

and location within the bay Flow events during unstratied conditions in western

Massachusetts Bay are particularly wellrepresented indicating that the model is

appropriate for studying processes such as the transport of suspended material

from the future outfall site due to winter storms Individual ow events during

stratied conditions and in the oshore Stellwagen Bank region however are less

wellrepresented due to small length scales caused by upwelling and river discharge

events coupled with insucient data to specify open boundary forcing from the

Gulf of Maine Thus while the model might be used to answer issues such as

the frequency with which Gulf of Maine river plumes visit the new outfall site

attempting to predict whether a particular plume would visit the outfall site could

be problematic

Comparative simulations of euent discharged from the existing and future

Boston outfalls show that the region of relatively high euent concentrations 

part euent to  parts sea water is signicantly smaller with the future outfall

and is limited to Western Massachusetts Bay during both unstratied and stratied

seasons The region of even higher concentration  part euent to  parts sea

water that covers much of Boston Harbor with the existing outfall is nonexistent

in the future outfall simulation Additional simulations of chlorination plant failure


predict that the oshore location of the future outfall will lead to dramatically

lower levels of pathogens at area beaches


Contents

Introduction 

 Rationale  

 Physical Setting  

 Motivation for Realistic Simulations  

 Organization  

 Model Implementation 

 Numerical Scheme  

 Model Geometry and Bathymetry  	

 Vertical Mixing  	

 Horizontal Mixing  


	 Surface and Bottom Stress  

 Surface Heat Flux  

 Freshwater Input  

 Open Boundary Conditions  

 Adjunct Western Gulf of Maine Model  	


 Boundary condition modications for  run  

 Time Step  

 SigmaCoordinate Correction  

 ModelData Comparison 

 Temperature Salinity and Density  

 Velocity Time Series Comparison at the Mooring Locations   

 Statistical Comparison of Currents at the Mooring Locations  

 Euent Modeling 	



	 Simulation of chlorination breakdown at the existing and fu

ture outfall locations 

 Summary 

A ModelData Comparison of Time Series Velocity at the Boston

Buoy 


B Mean Currents by Season 

C Comparative Euent Dilution Simulations for the Existing

and Future Outfalls 	

D References 


List of Figures

 The Massachusetts Bay region 


 Massachusetts Bay model grid 

 Data collection sites in Massachusetts Bay  

 Meteorological Data  


 Freshwater input to the Western Gulf of Maine  

	 Adjunct Western Gulf of Maine grid 

 Modeldata comparison of OBC driven ow  


 Inshore Massachusetts Bay boundary condition from WGOM

model 

 Oshore Massachusetts Bay boundary condition from WGOM

model 


 Observed and modeled vertical temperature dierence 

 Modeldata comparison of nearsurface salinity at Stations

BB BS SC and U 

 Modeldata comparison of lower layer salinity at stations BB

BS SC and U 	

 Modeldata comparison of lower layer salinity at stations U

MN RP and U
 

 Modeldata comparison of nearsurface temperature at sta

tions BB BS SC and U 


	 Modeldata comparison of nearsurface temperature at U

RP MN and U
 

 Modeldata comparison of lower layer temperature at BB BS

SC and U 

 Modeldata comparison of lower layer temperature at U MN

RP and U
 


	
 Modeldata comparison of vertical density dierences at sta

tions BB BS SC and U 

 Coherence between nearsurface salinity time series 


 Coherence between nearsurface temperature time series 

 Boston Buoy velocity comparison Jan    Mar 	 	

 Boston Buoy velocity comparison Jun  
  Aug  


 Boston Buoy ModelData velocity statistics at 	 m depth  


 Boston Buoy ModelData velocity statistics at  m depth  

	 ModelData velocity correlations at selected regional stations	


 ModelData comparison of the st EOF structure of the ve

locity time series  	


 Mean ow and lowfrequency variability ellipses for Winter


  	

 Mean ow and lowfrequency variability ellipses for Spring  	

 Mean ow and lowfrequency variability ellipses for Summer


  		


 Mean ow and lowfrequency variability ellipses for Fall 
  	

 Euent comparison at  m depth during Winter 
 


 Euent comparison at  m depth during Winter  

 Euent comparison at  m depth during Summer 
 

 Euent comparison at  m depth during Summer 
 

	 Residence time for euent in Massachusetts Bay and Boston

Harbor  

A Boston Buoy velocity comparison Jan  
  Feb  



A Boston Buoy velocity comparison Feb  
  Apr  



A Boston Buoy velocity comparison Apr  
  Jun  



A Boston Buoy velocity comparison Jun  
  Aug  

	

A	 Boston Buoy velocity comparison Aug  
  Oct  



A Boston Buoy velocity comparison Oct  
  Dec   



A Boston Buoy velocity comparison Dec  
  Jan  



A Boston Buoy velocity comparison Jan    Mar 	 


A Boston Buoy velocity comparison Mar 	   May 
 


A
 Boston Buoy velocity comparison May 
   Jul 	 

A Boston Buoy velocity comparison Jul 	   Sep  

A Boston Buoy velocity comparison Sep    Nov  


A Boston Buoy velocity comparison Nov    Dec 
 

A Boston Buoy velocity comparison Dec 
   Feb  	

A	 Boston Buoy velocity comparison Feb    Apr 
 

A Boston Buoy velocity comparison Apr 
   Jun 	 


A
 Boston Buoy velocity comparison Jun 	   Aug 
 

A Boston Buoy velocity comparison Aug 
   Oct 	 

A Boston Buoy velocity comparison Oct 	   Nov 
 


B Modeled mean ow at  and 
 m depth for Winter 


B Modeled mean ow at  and 
 m depth for Spring 
 

B Modeled mean ow at  and 
 m depth for Summer 
  

B Modeled mean ow at  and 
 m depth for Fall 
 	

B	 Modeled mean ow at  and 
 m depth for Winter 


B Modeled mean ow at  and 
 m depth for Spring  


B Modeled mean ow at  and 
 m depth for Summer   

B Modeled mean ow at  and 
 m depth for Fall  

B Modeled mean ow at  and 
 m depth for Winter 



B
 Modeled mean ow at  and 
 m depth for Spring  


B Modeled mean ow at  and 
 m depth for Summer   


B Modeled mean ow at  and 
 m depth for Fall  


B Modeled mean ow at  and 
 m depth for Winter  


C Euent comparison during Winter 
 at  m depth  


C Euent comparison during Winter 
 at  m depth  



C Euent comparison during Spring 
 at  m depth 


C Euent comparison during Spring 
 at  m depth 


C	 Euent comparison during Summer 
 at  m depth 


C Euent comparison during Summer 
 at  m depth 

C Euent comparison during Fall 
 at  m depth 

C Euent comparison during Fall 
 at  m depth 

C Euent comparison during Winter 
 at  m depth 

C
 Euent comparison during Winter 
 at  m depth 	

C Euent comparison during Spring  at  m depth 

C Euent comparison during Spring  at  m depth 


C Euent comparison during Summer  at  m depth 

C Euent comparison during Summer  at  m depth 



List of Tables

 Freshwater inputs for the model 

 Modeled and observed mean ow speed on the eastern ank

of Stellwagen Bank  


	 Locations monitored in chlorination failure simulations 


	 Bacteria concentrations for breakdown of existing plant 

	 Bacteria concentrations for breakdown of future plant 


Chapter 

Introduction

Rationale

A broad theme of the USGeological Survey USGS National Marine

and Coastal Geology Program is understanding geologic factors which inu

ence the longterm quality and preservation of environments in their natural

stateToward this goal the USGS is developing an integrated package of

wave circulation and sediment transport models that can be used to study

the movement of polluted material in coastal watersAs a rst step of this de

velopment we have implemented and tested a threedimensional circulation

model in Massachusetts Bay a prototypical region of interest

 Physical Setting

Massachusetts Bay is a semienclosed embayment that opens into the

Gulf of Maine at its eastern boundary Figure It is roughly 

 km

long 	
 km wide and has an average water depth of 	 mThe bay is

bounded on the east by Stellwagen Bank which rises to within 
 m of sea

surfaceBoston Harbor empties into western Massachusetts Bay providing

a signicant source of contaminants to the region

Circulation in Massachusetts Bay is driven by a combination of local and

remote processes that vary with season Geyer et al Throughout

the year water ows southward in the western Gulf of Maine and although

most of this current continues owing southward over the eastern ank of

Stellwagen Bank largely bypassing Massachusetts Bay a small branch ows

into the bay and drives a weak counterclockwise ow that enters at Cape


42o

00’

42o

30’

43o

00’N

70o71 00’W o00’ 70o30’

40

40

40

80

80

80

80

80

0 20 KM

Contours in meters

OUTFALL

New York

Boston

Halifax

Figure  Bathymetric map showing Massachusetts and Cape Cod Bays

present sewage outfalls in Boston Harbor diamonds location of future ocean

outfall for treated Boston sewage in western Massachusetts Bay and the

approximate boundary of the Stellwagen Bank Marine Sanctuary



Ann and exits at Race PointThe magnitude of this ow varies from less

than  cm s o Boston to about  cm s along the western shore of

Massachusetts BayThe remotely driven mean circulation pattern is modied by seasonal

current regimesDuring the winter the bay is wellmixed vertically and

the wind is principally from the northwest reinforcing the counterclockwise

circulationIn the spring surface warming and freshwater intrusions from

Gulf of Maine rivers cause the bay to become stratied and the surface

currents in the bay are dominated by the strong and variable density driven

ow associated with these intrusionsBy summer the stratication has

intensied due to further surface warmingThe dominant wind direction in

summer is from the southwest which drives strong upwelling events along

the western and northern shores of the bayIn the fall rapid cooling on the

shallow western side of the bay results in a density eld that temporarily

reverses the surface mean ow

Tidal currents in Massachusetts Bay are dominated by the semidiurnal

Mconstituent period   hoursThe currents are largely bidirec

tional and the magnitude ranges from about 
 cm sin the interior of

Massachusetts Bay to more than 	
 cm s o the tip of Cape Cod and in

the entrances to Boston Harbor and Plymouth Harbor Blumberg et al 

Irish and Signell Although the transport due to tides is generally os

cillatory and therefore doesnt give rise to net transport of material they are

an important source of bottomgenerated turbulent mixing

 Motivation for Realistic Simulations

A ma jor sewage treatment pro ject is currently underway in Boston prin

cipal components of which include upgrade to secondary treatment and relo

cation of sewage outfall pipes from Boston Harbor to a site  miles oshore

hereafter referred to as the future outfall siteDischarge at the future

outfall site is currently estimated to begin in 
The impact of the new

treatment pro ject on Massachusetts Bay is of critical concern to the Mas

sachusetts Water Resources Authority MWRA which is responsible for

the pro ject and to the citizens of MassachusettsSelection of the future out

fall location was based in part on steady twodimensional depthaveraged

modeling results that showed improved dilution of euent at the oshore

locationLacking in this model however was the ability to simulate ow


under stratied and timevariable conditions although attempts were made

to parameterize these eectsWith the pressing need to better understand

and simulate transport processes in Massachusetts Bay the region was chosen

by the USGS as a testbed for developing a prototype transport model

To address watercolumn eects of the future outfall the MWRA con

tracted HydroQual Inc in the Spring of  to develop a threedimensional

water quality model of the bayTo avoid duplication of eort since three

dimensional circulation is required by both sediment transport and water

quality models the USGS agreed as part of a cooperativelyfunded set of

studies with MWRA to supply HydroQual with the necessary hydrodynamic

calculations for their water quality model

October    July   was selected for testing and calibration of

the circulation and water quality models due to a large data collection eort

during this period by the Massachusetts Bays Program Geyer By

the fall of  the USGS had made signicant progress on hydrodynamic

modelingAs water quality results were not yet available comparative di

lution simulations were made to assess the dierence between the existing

and future outfall locationsUsing a conservative tracer the results sug

gested that the future outfall would result in a signicantly smaller region

aected by high concentrations of euentThese results played a signicant

role in the EPA environmental assessment of the impact of the future outfall

on endangered species on Stellwagen Bank EPA Although the im

pact of the outfall is now more realistically represented by the water quality

model which includes biochemical processes HydroQual and Battelle 	

the comparative dilution simulations still provide important insight into the

transport and dispersion of material in Western Massachusetts Bay and are

thus presented in this report

 Organization

This report is organized as follows The model and its conguration are

described in chapter The model is then compared to data in chapter  to

assess its capability to represent temperature salinity density and circulation

in Massachusetts BayThe modeled euent elds from the existing and

future outfall are presented in chapter Simulated chlorination plant failures

at the existing and future outfall sites are shown in chapter 	


Chapter 

Model Implementation

Numerical Scheme

The model used in this study is called ECOMsi a semiimplicit variant

of the threedimensional Estuary Coastal and Ocean Model ECOM de

scribed by Blumberg and Mellor 
ECOMsi was selected because it

includes a free surface nonlinear advective terms coupled density and ve

locity elds river runo heating and cooling of the sea surface a 	 level

turbulence closure scheme to represent vertical mixing Mellor and Yamada

 and is designed to easily allow realistic simulationsIn addition the

combination of orthogonal curvilinear coordinates in the horizontal plane and

sigmacoordinates in the vertical dimension allows grid renement in regions

of interest without sacricing the wellknown characteristics of Cartesian grid

schemes

The basic equations are expressed in a sigma coordinate system



z  

H  



where Hx y is the bottom topography and x y is the surface elevation

The basic governing equations are presented here in Cartesian coordinates to

facilitate discussionThe equations as expressed in curvilinear coordinates

may be found in Blumberg and Mellor 


The continuity equation is



 t 

uD

x 

vD

y 



 
 


the x momentum equation is

uD

t 

uD

x 

uvD

y 

u

 fvD  gD



x 



KMD

u



gD





x

Z



d

gD



D

x

Z







d Fx 

and the y momentum equation is

uD

 t 

uvD

x 

vD

y 

v

 f uD  gDy 



KMD

v



gD





y

Z



d

gD



D

y

Z







d Fy 

where  is the surface elevation u and v are the x and y components of

velocity D is the total water depth D  H    is the transformed vertical

velocity normal to sigma surfaces KM is the vertical eddy viscosity  is

the water density  is a reference water density and Fx and Fy are the

horizontal viscous terms dened by

Fx 



x 

AMD

u

x 



y AMD



u

y 

v

x 	

and

Fy 



x AMD



u

y 

v

x 



y 

AMD

v

y 

 

where AM is the horizontal viscosityThe parameterizations of KM and

AM are discussed in chapters  and  respectivelyThe model also

solves prognostically for temperature salt concentration used here to model

sewage euent turbulence kinetic energy and turbulence macroscale

ECOMsi diers from the Blumberg and Mellor 
 ECOM model in

that it uses a semiimplicit scheme for calculating the free surface there

fore avoiding the gravity wave CFL condition required by explicit schemes

eg Casulli 
This has the advantage that larger time steps may be

taken on the order of minutes rather than tens of secondsA potential

disadvantage of implicit schemes is that they more readily damp free wave

motions but in strongly forced and damped shallow regions such as Mas

sachusetts Bay the eect is smallThis was determined by halving the time


step and observing negligible dierences in simulation resultsAnother dis

advantage is that because the calculation of surface elevation requires solving

a large matrix equation at each time step and ecient solution of this equa

tion requires positive deniteness boundary conditions for elevation must be

formulated in matrix form and must not destroy the positive deniteness of

the matrixWe use a combination of clamped and gravity wave radiation

conditions on the open boundary made possible by the implementation of

the partiallyclamped formulation of Blumberg and Kantha 	 discussed

in chapter 

 Model Geometry and Bathymetry

The model was congured on a  x  horizontal curvilinear orthogonal

grid Figure The grid spacing ranges from approximately 

 m in

Boston Harbor to about 


 m along the open boundary with the Gulf of

MaineThe grid spacing in the vicinity of the future outfall site is approxi

mately 


 mThe grid extends well oshore of Massachusetts Bay to allow

exchange with the Gulf of Maine inside the model domain and extends far

enough north to include the Merrimack River a large source of fresh water

to the regionInitially 
 evenly spaced sigma levels were used in the vertical

dimension but nal runs utilized  additional nearsurface layers at  and

 of the water depth to better resolve river plumes and the surface mixed

layer

The bathymetry for the model was obtained by interpolating NOAA

sounding data on to the model grid using an inverse distance methodThe

resulting grid was Shapiro ltered Shapiro 
	 to remove  grid length

variability and a maximum depth of 
 m was set to eliminate the deep

complex topography oshore of Massachusetts BayThe minimum depth

was set to  m to avoid ooding and drying

 Vertical Mixing

The vertical mixing parameterization was the level 	 turbulence closure

scheme of Mellor and Yamada  with the extensions by Galperin et

al An important feature of these extensions include a lengthscale

limitation that prevents the mixing length from becoming too large during

strongly stratied conditionsThe minimum background mixing value

below which the eddy viscosity never falls was set to 	 
 msThis

	
 71°00’ 70°30’ 70°00’ 69°30’

 42°00’

 42°30’

 43°00’

10 km

Maximum grid size ~= 6000 m

Minimum grid size ~= 600 m

Figure  Model grid for the threedimensional circulation model of Mas

sachusetts and Cape Cod BaysThe curvilinear orthogonal grid allows the

mesh resolution to vary spatially having a minimum grid spacing of 

 m

and a maximum spacing of 


 mThe grid spacing in the vicinity of the

future outfall is roughly 


 mThe 
 m isobath is represented by the

thick grey line


value is comparable to the value of vertical mixing observed by Geyer and

Ledwell  from a dye experiment near the future outfall site

 Horizontal Mixing

The horizontal mixing parameterization used the Smagorinsky  for

mulation in which the magnitude of the horizontal mixing is proportional to

horizontal current shearFor the water quality runs the Smagorinsky coef

cient was set to 
 which results in typical horizontal eddy viscosities and

diusivities between 	  
 msover much of Massachusetts BayIdeally

the magnitude of this mixing would represent processes occurring at scales

smaller than the mesh can resolveIn western Massachusetts Bay the grid

cell spacing is  km and using the rule of thumb that  grid cells motions

are the minimum length that are realistically represented the model should

well resolve mixing processes occurring at scales of 
 km or soFrom Okubo


 we expect that observed horizontal mixing at a scale of 
 km should

be about 
 ms and at a scale of  km should be about  msSince the

model partially resolves shear and thus produces mixing at scales between

 and 
 km the ideal mixing level should probably be somewhere between

 and 
 msTo test the sensitivity to the mixing level in this range

the model was subsequently rerun for a  month period using a Smagorin

sky coecient of 

 ve times smaller than the water quality run and

only small quantitative dierences were observed in the resulting velocity

temperature salinity and euent concentration eldsThis indicates that

the modeled elds are dominated by resolved mixing and transport processes

rather than the parameterized mixing of the Smagorinsky formulationThis

is consistent with observations made during the outfall siting process where

it was concluded that a value of 	 ms best represented the dispersion of

material over the 
 km scale of western Massachusetts Bay Adams et al


 suggesting that the eective mixing in this region may be a bit higher

than Okubos empirical curveThis would not be surprising considering the

strong current shears in this regionThe level of 	 ms was also used in

the outfall siting model studies described by Walton et al 


In addition to the specied horizontal viscosity the velocity eld was

Shapiro ltered every  hours to remove  grid length energyThis was

necessary to prevent  grid length energy along the open boundary from

growing large enough to violate the advective CFL conditionBy applying a



th order Shapiro lter at  hour intervals the longer wavelength energy is

largely unaectedFor example it can be shown over 
 lter applications

two days of simulation the amplitude of  grid length structure is reduced

less than 	 see appendix of Signell 

 Surface and Bottom Stress

The surface wind stress was assumed to be uniform over the model domain

and was determined from wind measurements at the Boston Buoy Station

BB in Figure  using the Large and Pond  formulation

Although the scale of the bay is relatively small compared to most syn

optic scale wind features there is certainly some variability over the bay

caused in part by local eects such as the seabreezeSome of the discrepan

cies between model and data especially during the summer are undoubtably

the result of the idealized wind stress eldObservations were available at 

hourly intervals from the Boston Buoy anemometer height  m except

for the period 
  

 when wind data from Logan Airport was

used insteadA comparison between Boston Buoy and Logan Airport winds

Figure  Geyer et al  suggests that the winds are comparable dur

ing this time period although the wind at Boston Buoy is weaker during the

summer and stronger during the winter than winds at LoganThe hourly

averaged winds were converted to wind stress then averaged over  hourly

intervals before input to the model to reduce the size of the input le

The bottom stress is determined from the quadratic drag law applied at

the point closest to the bed where there is a velocity estimateThis is the

middle of the lowest grid cellBecause the distance between the the lowest

velocity point and the bed varies with location the bottom roughness length

zo is specied at each location and then a constant stress layer is assumed

so that an eective drag coecient can be computed at the velocity point

For this study zo was assumed to be uniform throughout the domain with a

value of 


 m which is equivalent to a drag coecient of 	
at 
 m

above the bottom which is in turn a typical value used for depthaveraged

models of coastal seasIn these runs the bottom sigma layer thickness is


 of the water depth so the lowest velocity point is 	 of the water depth

above the bedWithin Massachusetts Bay the maximum distance between

the lowest velocity point and the bed is  m in Stellwagen Basin where the

model water depth reaches  m


42o

00’N

42o

30’N

43o

00’N

70o71 00’W o00’ 70o30’

40

40

40

80

80

80

80

80

0 20 KM

Contours in meters

BB

BS

MN

SC

U7

RP

U6

U3

SB

U2

Figure  Data collection sites in Massachusetts BayTime series of cur

rents salinity and temperature were collected at the 
 locations indicated

Meteorlogical information was obtained from the Boston Buoy Station BB


Oct Feb Jun Oct Feb Jun Oct Feb Jun Oct

0

10

20

m/yr

1990 1991 1992

Precipation − Evaporation

−10

0

10

20

C

Air Temperature

0

50

100

%

Relative Humidity

950

1000

1050

mb

Barometric Pressure

0

500

1000

watts/m2

Insolation

−5

0

5

Surface Meteorological Boundary Conditions

dynes/cm2

Wind Stress

Figure  Meteorological data used to specify wind stress buoyancy ux

and heat ux boundary conditions at the sea surface



 Surface Heat Flux

The surface heat ux was calculated using techniques described by Weller

et al 	 and the ux was introduced into the top sigma layer at each

grid locationThese techniques produce bulk estimates of latent sensible

and longwave radiation given sea surface temperature air temperature in

solation relative humidity barometric pressure and windThe limited avail

ability of these data led us to use air temperature barometric pressure and

wind data at the Boston Buoy relative humidity at Logan Airport and inso

lation at Woods Hole to represent baywide conditions Figure The sea

surface temperature used in the heat ux estimates was initially determined

from observations at the Boston BuoyThis proved problematic because

Cape Cod Bay is signicantly warmer than the Boston Buoy site during

most of the summer which should result in reduced heat ux all else being

equal but there was no feedback mechanism to allow thisThe result was

that model temperatures in Cape Cod Bay became too warmTo address

this problem later runs used modelgenerated sea surface temperature elds

to calculate surface heat uxUsing the modeled sea surface temperature

eld allows spatial variability in seasurface temperatures to feed back into

the heat ux routinesUnder upwelling conditions for example there can be


or more dierence between surface water temperatures in the Bay which

can in turn result in heat ux dierences of 

 W mor moreAllowing

for spatial dierences in the surface heat ux due to varying water temper

atures led to a signicantly improved correspondence between modeled and

observed water temperatures in Cape Cod BayHourly estimates of heat

ux obtained from the bulk formulae were averaged over  hours to reduce

the size of the model input le

 Freshwater Input

Fresh water was introduced directly into the model at grid cell locations

representing the mouths of the Merrimack Charles Neponset and Mystic

rivers as well as at the existing sewage euent discharges at Nut and Deer

IslandsIn addition input from the Penobscot Androscoggin Kennebec

and Saco Rivers were used to drive an adjunct model of the western Gulf of

Maine that was used to help improve salinity open boundary conditions for

the Massachusetts Bay model described hereThe western Gulf of Maine

adjunct model is described in the open boundary condition chapter 



 Discharge Agency Location Scale

 m s Factor

Merrimack  USGS Lowell MA 


Charles 	 USGS Waltham MA 


Deer Island  MWRA Deer Island 



Nut Island 
 MWRA Nut Island 



Neponset 	  Inferred from Charles 



Mystic 
  Inferred from Charles 



Penobscot 		 USGS Eddington ME 


AndroKenn  USGS NSidneyAubern ME 	

Saco  USGS Cornish ME 

Table  Average freshwater inputs for the model during 
The scale

factor is the amount that the gauged discharge was multiplied by to account

for drainage area downstream of the gauge

To indicate the relative importance of these sources the mean ow during


 is shown in Table  and indicates that the Penobscot Androscog

ginKennebec and Merrimack River systems dominate the freshwater input

to this regionThe Androscoggin and Kennebec Rivers are combined because

they merge before discharge into the Gulf of MaineThe daily ows from

the Neponset and Mystic Rivers were inferred by multiplying the daily ow

from the Charles River by 
 and 
	 respectively factors determined

by the ratio of the annual average

Examination of the daily discharge over the three year simulation period

shows that although the maximum discharge occurs during the spring of

each year there are secondary maxima in the fall Figure Fall 
 in

particular has unusually high discharge almost as large as the preceding and

following spring discharges

	 Open Boundary Conditions

Open boundary conditions are required for elevation temperature salin

ity and tracer concentration used here to represent euentBecause the

elevation eld at each time step is solved for with a fully implicit scheme the

boundary condition must be specied as part of the matrix solution for the


Gulf of Maine inputs

Boston inputs 

Merrimack 

Saco 

Andro−Kenn

Penobscot 

Rivers: Charles+Neponset+Mystic

Effluent: Deer + Nut 

0

500

1000

1500

2000

2500

3000

m3/s

Major Rivers

0

1000

2000

3000

4000

5000

6000

m3/s

Freshwater Input

Gulf of Maine & Boston inputs

Oct Apr Oct Apr Oct Apr Oct

0

20

40

60

80

m3/s

1990 1991 1992

Boston Sources

Figure  Freshwater input to the Western Gulf of MaineThe top panel

shows the combined discharge due to the Kennebec Androscoggin Penobscot

and Merrimack Rivers compared to the Boston dischargeThe middle panel

shows the individual contributions of the ma jor rivers and the lower panel

show the relative contributions of rivers and euent to discharge in Boston

Harbor


elevation eldIn addition the matrix solver used in this code requires the

matrix to be symmetric and positive denite which puts further restraints

on the type of boundary condition imposedThe partially clamped bound

ary condition of Blumberg and Kantha 	 can be shown to satify these

conditions and is expressed as



t  q

gD 

x 

  data

Tlag



where data is specied from dataIf Tlag is very large a pure radiation

condition is approached and if Tlag is very small a purely clamped speci

ed condition is approachedIn the Massachusetts Bay model the oshore

boundary elevation was clamped while the southern boundary radiated grav

ity wavesThe values for elevation used along the oshore boundary were a

combination of M N and Stidal elevations derived from the tidal model

of Lynch and Naimie  and lowfrequency uctuations obtained from

the adjunct Gulf of Maine model described in chapter The Mtidal

elevations obtained from Lynch personal communication  were used

to synthesize Sand Nconstituents by scaling the Mamplitudes by the

observed SMand NMratios at the Boston Buoy and shifting the S

and Nphases by the amount observed at the Boston BuoyThis is a rea

sonable approximation since the M Sand Namplitude and phase vary

in much the same manner in the Gulf of Maine egMoody et al 

Temperature salinity and euent are also specied along the open bound

ariesIf the ow is directed out of the domain the interior values are simply

advected out of the domainWhen outow turns to inow the water prop

erty values slowly move toward specied values over a dened relaxation time

to avoid articial fronts from developingIn the Massachusetts Bay model

the relaxation time was specied to vary linearly from  days at the most

northern boundary cell to 
 days at the southern boundary cell

Initially the temperature and salinity values along the boundary were

specied by climatology interpolated from the Bedford Institute of Oceanog

raphy AFAP database Drinkwater  while the euent concentration

eld was set to zeroAlthough it was hoped that most local river eects

would be due to the Merrimack and the Charles early runs showed that the

gross salinity picture for Massachusetts Bay could not be adequately specied

by parameterizing the other Gulf of Maine rivers by seasonal climatologyTo


address this problem an adjunct model of the western Gulf of Maine was de

veloped to provide better salinity boundary conditions for the Massachusetts

Bay model

 Adjunct Western Gulf of Maine Model

The western Gulf of Maine model grid is 
 x 	 cells and has a typical grid

spacing of  km Figure 	The bathymetry was generated in identical

fashion to the Massachusetts Bay model and the model was driven by the

same wind and heat ux forcingThe elevation boundary conditions were

specied as M Sand Ntides derived from Lynch and Naimie 

There were no subtidal uctuations specied so the lowfrequency elevation

was eectively clamped at zero along the oshore boundaryAt the northern

boundary a mean sea surface slope is imposed to drive a coastal current

Along this boundary the lowfrequency elevation varied linearly from 



m at the coast to 


 m at the open boundary 	 km oshoreThis was

picked to drive a mean alongshore ow of 	
 cm s at locations where

mean currents were observed by Vermersch et al 
The mean ow in

the Gulf of Maine model is held constant throughout the yearRecent data

collected on the eastern ank of Stellwagen Bank at a location determined to

be a good measure of the largescale circulation in the Gulf of Maine Bogden

et al 	 suggests encouragingly that there might be little variability in

the seasonal mean owA comparison of the mean ow observed in spring

summer and fall of  at the Stellwagen Bank mooring with seasonally

averaged mean ow obtained from the model at this location over 


shows remarkable agreement Figure Outside of the buoyant surface

layer the observed mean ow due to large scale circulation in the Gulf of

Maine changes very little with season Table The mean ow at 
 m

for example appears to be about 
 cm s throughout the year and the

mean ow at 	
 m is nearly constant at  cm sThe temperature and salinity along the open boundary were specied by

monthly climatology determined from the Bedford Institute of Oceanography

AFAP databaseModel output averaged over the tidal cycle was output

along the cells that constituted the open boundary for the Massachusetts

Bay model the converted into lowfrequency elevation boundary conditions

and temperature and salinity boundary conditions at standard levels  	

	
 
	 
 	 
 
 
 	
 
	 and 
 mThe resulting boundary

	
 71°00’ 70°30’ 70°00’ 69°30’ 69°00’ 68°30’

 42°00’

 42°30’

 43°00’

 43°30’

 44°00’

 44°30’

Androscoggin/Kennebec Rivers −>

Penobscot River −>

Saco River −>

River −>

Merrimack

<−− Western Gulf of Maine Grid

<−− Massachusetts Bay Grid

Figure 	 Western Gulf of Maine grid showing ma jor riversGrid spacing

is  kmThe isobaths shown are 
 

 and 

 m


10 cm/s

Win Spr Sum Fa ll

5 m

27 m

50 m

Figure  Modeldata comparison of mean ow on the east ank of Stell

wagen BankThe thick grey arrows are observed means from  while

the thin black arrows are modeled means from 
  and The

ow at this location is strongly controlled by largescale circulation in the

Gulf of Maine that is parameterized in the model through the open boundary

conditions

MODEL 
 DATA 

cm s cm s

	 m 
 m 	
 m 	 m 
 m 	
 m

Winter 
 
  nd nd nd

Spring 
    

 

Summer 	 
   
 

Fall    	 
 

Table  Modeled and observed mean ow speed on the eastern ank of

Stellwagen Bank at 	 
 and 	
 m depthsBelow the river plume inuenced

surface layer the mean ow at 
 m is relatively constant at 
 cm s



conditions for Massachusetts Bay as intended represent more variability

and Gulf of Maine dischargedriven events than would have been possible

with open boundary conditions based on climatology Figures  

Interannual variability at frequencies less than the synoptic scale weather

band 
 days is not dramatic although certain events stand out such as

the sustained decrease in surface salinities near the coast during the fall of


 consistent with the abnormally large runo during this time period

 Boundary condition modifications for 

 run

Although the temperature and salinity boundary conditions obtained

from the climatology and the western Gulf of Maine adjunct model result

in a quite realistic representation of the observed temperature and salinity

structure in Massachusetts Bay actual data obtained near the boundary dur

ing the simulation time period can further improve the representationTo

improve the  water quality verication run several sets of CTD data

were obtained from the National Marine Fisheries Service HolzwarthDavis

and Taylor  and used to make very lowfrequency weeks to months

corrections to the temperature and salinity conditions

Time Step

To achieve optimal eciency for these multiyear runs the time step was

varied with season and forcing conditions as determined by trial and error

Although ECOMsi is not restricted by the gravity wave CFL condition it

is limited by the advective CFL conditionDuring the spring and summer

months in particular the time step needed to be decreasedDuring the

winter the time step was usually  s 
 steps per Mtidal cycle except

for a few short periods generally coinciding with winter stormsDuring the

summer the time step was usually 

 s  steps per Mtidal cycle

 SigmaCoordinate Correction

In regions of steep topography sigmacoordinate models are prone to er

rors in the calculation of pressure gradient and horizontal diusion terms

Beckman and Haidvogel In the modeling here initial stratication

results were not encouraging in that the modeled temperature stratication

during the summer was 	
 weaker than observed at the future outfall site

Figure  and in Stellwagen BasinAlthough vertical mixing was initially


1 m 

20 m

Oct Apr Oct Apr Oct Apr Oct

24

26

28

30

32

34

psu

1990 1991 1992

Salinity

1 m 

20 m

0

5

10

15

20

25

C

Temperature

Open Boundary Condition 10 km off New Hampshire Coast (i=67,j=59)

Figure  Inshore Massachusetts Bay boundary condition from WGOM

model


1 m 

20 m 

140 m

1 m 

20 m 

140 m

0

5

10

15

20

25

C

Temperature

Open Boundary Condition 50 km off Cape Ann (i=67,j=39)

Oct Apr Oct Apr Oct Apr Oct

24

26

28

30

32

34

psu

1990 1991 1992

Salinity

Figure  Oshore Massachusetts Bay boundary condition from WGOM

model



15

10

5

0

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

Celcius

1990

Vertical temperature

difference (5 m - 23 m)

at the new outfall

location

model without

sigma correction

model with

sigma correction

observed data

Figure  Observed and modeled temperature dierence 	 m   m at

Station BB see Figure Before the correction the model greatly un

derpredicted the degree of temperature stratication in the summer months

Observed data is missing from late August to November


suspected the problem turned out to be caused by spurious upslope mix

ing caused by coordinate errorsThese errors were signicantly reduced

by removing the domainaveraged vertical salinity temperature and density

proles before calculating horizontal dierences egBeckmann and Haid

vogel The eectiveness of this technique depends upon the degree

to which the vertical structure varies over the domainIn domains such

as Massachusetts Bay in which the stratication is nearly uniform over the

region this technique greatly improves the model prediction in comparison

to observed data Figure The more the vertical structure varies over

the domain however the less useful the technique in reducing coordinate

errorsIn domains where both stratied and wellmixed regions exist er

rors would actually increase in the wellmixed region since subtracting the

domainaveraged vertical prole would result in more vertical structure lo

cally


Chapter 

ModelData Comparison

Temperature
 Salinity and Density

The comparison between modeled and observed temperature salinity and

density at most sites shows that the model generally does a good job at re

producing seasonal cycles and top to bottom density dierences Figures 

There are some clear dierences howeverWhile most of surface salinity

variations during the spring and summer are well represented Figure 

the salinity drop associated with the large fall 
 discharge appears in the

model about three weeks too late and during spring of  the surface

salinities in the model are 
	
 psu too saltyThe lower layer salinities

Figures  and  during this time period are also 
	
 psu too salty

so the vertical density structure is minimally aectedIt does suggest how

ever that if a better match of the salinity magnitude was desired it would be

necessary to improve upon the climatological salinity boundary conditions

for the western Gulf of Maine modelThese boundary conditions specied

a xed seasonal cycle assuming that all interannual variability was due to

local river discharge in the western Gulf of MaineIt is reasonable to ex

pect that interannual variability in discharge from the StLawrence river

and variations in the Gulf of Maine circulation would also be responsible for

interannual salinity uctuations in the water entering Massachusetts BayThe modeled temperatures and salinities appear to track some of the

observed shorter period uctuations in addition to the general seasonal cycle

The coherence between modeled and observed water properties is generally

signicant at 	 condence over periods of 
 days Figures  and 



data 

model

30

31

32

psu

Broad Sound (BS)

30

31

32

psu

Scituate (SC)

data 

model

30

31

32

Near−Surface Salinity (4−5 m)

psu

Boston Buoy (BB)

Jun Aug Oct Dec Feb Apr Jun

30

31

32

psu

1990 1991

Stellwagen Basin (U6)

Figure  Modeldata comparison of nearsurface salinity at Stations BB

BS SC and U


data 

model

31

32

33

psu

Broad Sound (BS) at 19 m

31

32

33

psu

Scituate (SC) at 24 m

data 

model

31

32

33

Lower−Layer Salinity

psu

Boston Buoy (BB) at 23 m

Jun Aug Oct Dec Feb Apr Jun

31

32

33

psu

1990 1991

Stellwagen Basin (U6) at 45 m

Figure  Modeldata comparison of lower layer salinity at stations BB

BS SC and U

	
data 

model

31

32

33

psu

Manomet (MN) at 25 m

31

32

33

psu

Race Point (RP) at 55 m

data 

model

Jun Aug Oct Dec Feb Apr Jun

31

32

33

psu

1990 1991

Cape Cod Bay (U7) at 25 m

31

32

33

Lower−Layer Salinity

psu

Cape Ann (U2) at 60 m

Figure  Modeldata comparison of lower layer salinity at stations U

MN RP and U



data 

model

5

10

15

20

Near−Surface Temperature (4−5 m)

C

Boston Buoy (BB)

5

10

15

20

C

Broad Sound (BS)

5

10

15

20

C

Scituate (SC)

data 

model

Jun Aug Oct Dec Feb Apr Jun

5

10

15

20

C

1990 1991

Stell Basin (U6)

Figure  Modeldata comparison of nearsurface temperature at stations

BB BS SC and U



data 

model

5

10

15

20

Near−Surface Temperature (4−5 m)

C

Cape Ann (U2)

5

10

15

20

C

Race Point (RP)

5

10

15

20

C

Manomet (MN)

data 

model

Jun Aug Oct Dec Feb Apr Jun

5

10

15

20

C

1990 1991

Cape Cod Bay (U7)

Figure 	 Modeldata comparison of nearsurface temperature at U RP MN and U



data 

model

5

10

15

20

Lower−Layer Temperature

C

Boston Buoy (BB) at 23 m

5

10

15

20

C

Broad Sound (BS) at 18 m

5

10

15

20

C

Scituate (SC) at 23 m

data 

model

Jun Aug Oct Dec Feb Apr Jun

5

10

15

20

C

1990 1991

Stellwagen Basin (U6) at 45 m

Figure  Modeldata comparison of lower layer temperature at BB BS

SC and U


data 

model

5

10

15

20

Lower−Layer Temperature

C

Cape Ann (U2) at 60 m

5

10

15

20

C

Manomet (MN) at 25 m

5

10

15

20

C

Race Point (RP) at 55 m

data 

model

Jun Aug Oct Dec Feb Apr Jun

5

10

15

20

C

1990 1991

Cape Cod Bay (U7) at 25 m

Figure  Modeldata comparison of lower layer temperature at U MN

RP and U




data 

model

0

1

2

3

4

Vertical Density Difference

sig−t

Boston Buoy (BB)

[23 m − 5 m]

0

1

2

3

4

sig−t

Broad Sound (BS)

[18 m − 5 m]

0

1

2

3

4

sig−t

Scituate (SC)

[23 m − 5 m]

data 

model

Jun Aug Oct Dec Feb Apr Jun

0

1

2

3

4

sig−t

1990 1991

Stell Basin (U6)

[45 m − 4 m]

Figure  Modeldata comparison of vertical density dierences at stations

BB BS SC and U


The temperature coherence decreases somewhat with increasing frequency

but the salinity coherence doesnt show any obvious dependence on frequency

Model uctuations in salinity are comparable in strength to observations

Figure  except over Stellwagen Basin U where the modeled salinity

uctuations are less that half the strength of the observationsThe strong

observed uctuations of nearsurface salinity at U Figure  are associated

with intrusions of lowsalinity water from the Gulf of Maine Geyer et al

 and it appears that these intrusions are either underrepresented or

are somewhat smeared out by the model

Model uctuations in temperature are comparable in strength to observa

tions Figure 
 except at Scituate SC and Broad Sound BS where the

modeled temperature uctuations are about 	
 stronger than the observa

tionsThe strong modeled uctuations of nearsurface temperature at these

nearshore locations are associated with upwelling events so one explanation

is that the model surface layer responds a bit too strongly to upwelling winds

 Velocity Time Series Comparison at the Mooring

Locations

The quality of the comparison between modeled and observed currents

varies substantially with space and timeAt the Boston Buoy for example

events are frequently wellrepresented during the winter Figure  but

are only occasionally wellrepresented during the summer Figure As

explained in Signell et al  this is because the ow regime in winter is

largely winddriven with large spatial scalesThus small deviations in the

structure or timing of the circulation between between the model and reality

will not result in large errors at the mooring locationsDuring the stratied

seasons however the spatial scales are considerably smaller due to upwelling

intrusions of Gulf of Maine river water etcSmall deviations in the structure

or timing during the stratied seasons will therefore result in larger errors

at the mooring locationsFor example if the model representation of an

intrusion is o by 	 km the whole signal could be missed at a particular site

like the future outfall location

The degree of agreement between model and data can be quantied by

the complex correlation coecient which is high during the winter and con

siderably lower during the summer Figures  and Although the


0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0

0.2

0.4

0.6

0.8

1

Frequency (cycles per day)

Coherence: Modeled/Observed Near−Surface Salinity (4−5 m)

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0

0.5

1

1.5

2

Frequency (cycles per day)

Amplitude Ratio: Modeled/Observed Near−Surface Salinity (4−5 m)

BS

U6

SC

BB

BB

BS

SC

U6

Figure  Top Panel Coherence between modeled and observed near

surface salinities at Broad Sound BS Boston Buoy BB Scituate SC

and Stellwagen Basin UThe 	 condence level is 
Lower Panel

The amplitude ratio of model salinity uctuations to observed salinity uc

tuations


0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0

0.2

0.4

0.6

0.8

1

Frequency (cycles per day)

Coherence: Modeled/Observed Near−Surface Temperature (4−5 m)

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0

0.5

1

1.5

2

Frequency (cycles per day)

Amplitude Ratio: Modeled/Observed Near−Surface Temperature (4−5 m)

U6

BS

BB

SC

SC

BB

U6

BS

Figure 
 Top Panel Coherence between modeled and observed near

surface temperatures at Broad Sound BS Boston Buoy BB Scituate

SC and Stellwagen Basin UThe 	 condence level is 
Lower

Panel The amplitude ratio of model temperature uctuations to observed

temperature uctuations


−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

28 4 11 18 25 4 11 18 25

−10

0

10

cm/s

Jan Feb Mar

1991

MODEL: 23 m

Figure  Boston Buoy velocity comparison Jan    Mar 	 

	
−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

18 25 2 9 16 23 30 6 13

−10

0

10

cm/s

Jun Jul Aug

1990

MODEL: 23 m

Figure  Boston Buoy velocity comparison Jun  
  Aug  



0

0.2

0.4

0.6

0.8

1

Statistical Comparison of Model/Data Velocity at BB_05

Complex Correlation

Data 

Model

0

5

10

15

cm/s

Magnitude of mean current

Data 

Model

Dec Jun Dec Jun Dec Jun Dec

0

5

10

15

cm/s

1990 1991 1992

Magnitude of standard deviation

Figure  Boston Buoy ModelData velocity statistics at 	 m depth



Data 

Model

0

5

10

15

cm/s

Magnitude of mean current

Data 

Model

Dec Jun Dec Jun Dec Jun Dec

0

5

10

15

cm/s

1990 1991 1992

Magnitude of standard deviation

0

0.2

0.4

0.6

0.8

1

Statistical Comparison of Model/Data Velocity at BB_23

Complex Correlation

Figure  Boston Buoy ModelData velocity statistics at  m depth


correlations are considerably lower in the summer the modeled mean and

standard deviation of the lowfrequency currents are comparable to observed

values and do not show any clear bias

The relatively good correlation between observed and modeled currents at

the Boston Buoy is also seen at the other western Massachusetts Bay stations

Scituate and Manomet Figure 	The correlations are signicantly lower

in Cape Cod Bay and eastern Massachusetts BayAnother way of assessing the degree to which the model compares with

data is to look at the empirical orthogonal functions EOFs of the model and

dataAt mooring locations where at least 
 of the time period May 


 June  contained good data time domain EOFs were computed from

the eastwest and northsouth components of the subtidal velocity time

seriesThe rst modes of the observed and modeled velocities explain 	


and  of the variance respectively and have reassuringly similar structure

Figure This indicates that surface currents in western Massachusetts

Bay constitute the largest amount of the energy in both the model and the

data

 Statistical Comparison of Currents at the Mooring

Locations

Reproducing the day to day variations of ocean weather at specic

locations may be required to warrant the use of the model to address some

issues but for other issues it may only be necessary to correctly reproduce the

statistics or climate of the currents over a few weeks or monthsTo assess

the capability of the model to represent the climate the mean current and

lowfrequency ellipses from the data and model were compared at the moored

locations for each seasonThe seasons are dened by distinct dynamical

regimes instead of the conventional dates Winter November  to March 

represents a winddominated wellmixed regime Spring March  to June

 represents a runoinuenced transitional stratication regime Summer

June  to September  represent an upwelling strongly stratied regime

and Fall Sepember  to November  represents a transitional overturning

regime

During the winter as expected from the relatively high correlations at

the moorings the mean and lowfrequency ellipses compare quite well Fig

ure 
The ellipses are shown centered around the tip of the mean ow


BB_05

MN_05

SC_05

BB_23

MN_29

SC_16

U7_04

U7_25

0

0.5

1

W. Mass Bay Upper

0

0.5

1

Cape Cod Bay

0

0.5

1

W. Mass Bay Lower

0

0.5

1

E. Mass Bay Upper

Jul Sep Nov Jan Mar May

0

0.5

1

1990 1991

E. Mass Bay Lower

RP_05

U6_08

U3_04

U2_04

RP_23

U6_25

U3_25

U2_25

Figure 	 ModelData velocity correlations at selected regional stations

	

220 240 260 280 300

90

100

110

120

130

140

150

160

170

DATA: 1990/04/24−1991/06/15

220 240 260 280 300

90

100

110

120

130

140

150

160

170

MODEL: 1990/04/24−1990/04/24

Figure  ModelData comparison of the st EOF structure of the velocity

time series

	
220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

Winter 1990−1991: 1990/10/31−1991/02/28

Wind Stress

117

100

94

117

65

73

73

120

Upper Layer (4−8 m)

Data

Model

10 cm/s

1 dyne/cm2

220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

km

Wind Stress

88

74

94

65 101

120

Lower Layer (18−25 m)

1 dyne/cm2

Figure 
 Mean ow and lowfrequency variability ellipses for Winter 




	
and the arrows and ellipses are scaled to correspond to the distance that par

ticles traveling at the observed velocity would move in  dayThis means that

about twothirds of the time the percentage of time particles move less than

one standard deviation particles released at the moored locations would

be found somewhere inside the ellipse after  dayModeled displacements

are very close to predicted displacements with the exception of Race Point

where the observed mean ow is close to northward while the modeled ow

is close to eastwardAnalysis indicates that this discrepancy is most likely

due to a shift in the modeled pattern of tideinduced residual currents rather

than a ma jor error in representing ow out of Massachusetts BayThere

is strong 
 cm s northward mean ow in the model just  km to the

southwest the result of a small 	 km tideinduced residual eddy o the tip

of the Cape

During the spring the surface current variability increases due to intru

sions of low salinity water from Gulf of Maine rivers and weaker frictional

resistance to wind driving since the surface layer becomes somewhat decou

pled from the lower layerThe model reects the increased intensity of surface

currents and the general pattern of the current climate doing particularly

well at the future outfall site Figure The modeled surface currents

are slightly stronger than the observations near Stellwagen Bank and slightly

weaker than the observations at Scituate and Manomet howeverAt Race

Point there is a discrepancy between the modeled and observed mean ows

very similar to winter consistent with the interpretation that this is due to

slight spatial error in representation of the tidal residual eddy o the tip of

the CapeAs with the surface currents the lower layer currents are also

stronger than observed near Stellwagen Bank but are close to the observed

currents at the other locations except Scituate where the modeled currents

are extremely weak

In the summer strong stratication results in much stronger currents

in the upper layer than in the bottom layerThe model represents this

behavior quite well and again is a close match to the current climate at

the future outfall site Figure The model currents are again stronger

than observed near Stellwagen Bank and Cape Cod Bay surface currents are

also stronger than observed

In the fall the surface mean ow reverses and ows northward at Scituate

and Manomet in both the observations and the model while the mean ow

near Stellwagen Bank still ows to the south Figure 
The ellipses are

	
220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

Spring 1991: 1991/02/28−1991/05/31

Wind Stress

92 92

76

92

86

92 92

92

Upper Layer (4−8 m)

Data

Model

10 cm/s

1 dyne/cm2

220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

km

Wind Stress

92

92

35

92 72

92

92

92

Lower Layer (18−25 m)

1 dyne/cm2

Figure  Mean ow and lowfrequency variability ellipses for Spring 

	
220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

Summer 1990: 1990/05/31−1990/08/31

Wind Stress

76

38

92

35

35

86

92

Upper Layer (4−8 m)

Data

Model

10 cm/s

1 dyne/cm2

220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

km

Wind Stress

89

89

61

92

42

80

92

Lower Layer (18−25 m)

1 dyne/cm2

Figure  Mean ow and lowfrequency variability ellipses for Summer




		
220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

Fall 1990: 1990/08/31−1990/10/31

Wind Stress

48

48

58

57

59

61

61

61

Upper Layer (4−8 m)

Data

Model

10 cm/s

1 dyne/cm2

220 240 260 280 300 320

90

100

110

120

130

140

150

160

170

180

km

km

Wind Stress

58

51

48

48

57 61

61

Lower Layer (18−25 m)

1 dyne/cm2

Figure 
 Mean ow and lowfrequency variability ellipses for Fall 


	
generally quite similar except at Race Point where there appears to be some

ow into Massachusetts Bay that is not apparent in the data

The degree to which the model compares with data can be used to iden

tify appropriate uses of the modeled hydrodynamicsFor example since the

current climate is well represented in western Massachusetts Bay through

out the years it is appropriate to use the model to study the nearcontinuous

release of euent in this region by the existing and future outfallsLikewise

since individual events are wellrepresented in western Massachusetts Bay

during unstratied conditions the model is well suited to studying the trans

port of suspended material from this region during individual winter storms

On the other hand predicting the detailed evolution of a lowsalinity plume

from a spring runo event is not possible at least partially due to inadequate

information about the large scale forcing from the Gulf of Maine

	

Chapter 

Euent Modeling

Since the hydrodynamic model represents the stratication and statistics of

lowfrequency currents in the western part of Massachusetts Bay rather well

it is appropriate to use the model to investigate the transport of euent

from the existing and future outfalls over timescales of weeks to monthsTo

model the euent a conservative tracer was added to the code modeled the

same as temperature and salinity but with no dynamical eectThe euent

concentration was set to 

 at the Nut and Deer Island diusers for the

existing outfall runFor the future outfall run the discharge from Nut and

Deer Islands was added together and introduced at the future outfall again

with a concentration of 

Along the open boundary far from the sources

the euent concentration was set to zeroSince the loading concentration is

arbitrary the euent concentrations were converted to euent dilutions for

the purposes of presentation

Since ECOMsi is a hydrostatic model with grid spacing of about  km in

the vicinity of the future outfall it cannot represent the physics of the tur

bulent entrainment process that occurs when the euent is discharged from

the diusersBy designing the grid so that the size of the diuser grid cell is

comparable to the area of zone of initial dilution however it is possible to

dilute the euent by approximately the same amount as predicted by near

eld models such as EPAs ULINEThis allows ECOMsi to eectively model

the height to which the plume rises Blumberg et al PhD work by

XueYong Zhang in preparation at MIT has shown that with the relatively

strong pycnocline observed in Massachusetts Bay there is considerable lati

tude in choosing the size of the diuser grid cell to produce approximately

	
the right trap height

Euent was discharged and tracked from the existing and future outfalls

from October  to September  and the results have been presented in

a variety of formsSome animation clips comparing the two outfall locations

during winter and summer periods can be seen on the World Wide Web

at httpcrustyerusgsgovThe 

 dilution level is about the level at

which the nutrient signal from the outfall should become lost in background

variability and thus is highlighted in presentations

During the winter highest euent concentrations from both the existing

and future outfalls occur at the surface since the initially diluted euent

is lighter than the wellmixed dense seawater Figure The existing

outfall due to the shallow depth and connes of Boston Harbor results in a

much larger region of poorly diluted euentThere is very little dierence

in the size of the 

 dilution region between 
 and  although the

size of the 

 dilution region changes considerably Figures  and 

During the summer the highest euent concentrations from the existing

outfall are again found at the surface since the euent is eectively discharge

directly into the surface layer due to strong tidal mixing in Boston Harbor

The 

 dilution region is very similar to the winter region Figure 

At the future outfall however the euent plume becomes trapped below

the pynocline since the seawater near the surface has warmed and freshened

to the point where it is lighter than the initially diluted euentAt  m

the depth of maximum areal extent of the 

 euent surface the region

aected by the future outfall is still considerably smaller than the surface

extent of the 

 region from the existing outfall Figure Average

dilution comparisons for other seasons are presented in Appendix C

While the dispersion of euent can be characterized by the plots of the

average euent dilution a residence time for euent in Massachusetts Bay

and Boston Harbor can be determined by dividing the total amount of euent

in the bay or harbor by the rate of loadingThis calculation of residence time

assumes that both the input euent and export from the system are relatively

steadyThe boundary of Massachusetts Bay was specied as a straight line

between Cape Ann and Provincetown and the boundary of Boston Harbor

was specied as a straight line across the harbor mouthThe residence times

for Massachusetts Bay are very similar for the existing and future outfall

scenarios and show a distinct seasonal cycle Figure 	

Overall the residence time for euent in Massachusetts Bay is from 


	
220 240 260 280 300 320

80

90

100

110

120

130

140

150

160

170

180

km

Mean Effluent Dilution: 1990/01/02−1990/03/05

200

400

400

800

800

800

EXISTING OUTFALL

DEPTH = 2 m

220 240 260 280 300 320

80

90

100

110

120

130

140

150

160

170

180

km

km

400

400

800

800

800

800

FUTURE OUTFALL

DEPTH = 2 m

Figure  Euent comparison at  m depth during Winter 




220 240 260 280 300 320

80

90

100

110

120

130

140

150

160

170

180

km

Mean Effluent Dilution: 1990/11/08−1991/03/12

200

400

800

800

EXISTING OUTFALL

DEPTH = 2 m

220 240 260 280 300 320

80

90

100

110

120

130

140

150

160

170

180

km

km

400

800

800

FUTURE OUTFALL

DEPTH = 2 m

Figure  Euent comparison at  m depth during Winter 


220 240 260 280 300 320

80

90

100

110

120

130

140

150

160

170

180

km

Mean Effluent Dilution: 1990/06/06−1990/09/07

400

800

800

800

800

EXISTING OUTFALL

DEPTH = 2 m

220 240 260 280 300 320

80

90

100

110

120

130

140

150

160

170

180

km

km

FUTURE OUTFALL

DEPTH = 2 m

Figure  Euent comparison at  m depth during Summer 



Figure  Euent comparison at  m depth during Summer 



Existing Outfall

Future Outfall 

20

40

60

80

100

120

Residence Time of Effluent in Massachusetts Bay and Boston Harbor

Days

Massachusetts Bay

5

6

7

8

9

10

Days

Boston Harbor

Jan Jun Nov Apr Sep Feb Jul

10

15

20

25

30

m3/s

1990 1991 1992

Effluent Discharge Rate

Figure 	 Residence time for euent determined by dividing the total

amount of euent in a the bay or harbor by the euent loading rate




 days while the residence time in Boston Harbor is 	 daysThe euent

has longer residence times in the winter and shorter residence times in the

summer in both the harbor and the baySince the euent discharge does not

have a clear seasonal cycle the change in residence time is presumably due to

increased ushingFor the bay the residence time calculated via this method

changes over time scales that are comparable to or less than the residence

time itself indicating that an interpretation as total amount of euent in the

bay is more appropriate since the loading is relatively constantEuent

levels in the bay build during the fall and winter reach a peak in March

then fall precipitously in April as the bay is ushed by energetic springtime

currentsThe lowest amounts of total euent are found in summer and as

the water column destraties the euent levels begin to climb again

	
Chapter 

Simulation of chlorination

breakdown at the existing and

future outfall locations

One of the design considerations for the future outfall site was that it should

be more than one tidal excursion oshore so that in the event of a chlorination

failure untreated euent would not be swept into Boston area beaches on

the incoming tideSignell et al showed that even at spring tide

the future outfall location is well oshore of the region from which Boston

Harbor draws water on ood tideTo quantify the expected benet of the

oshore locations however chlorination failures were simulated in ECOM

si and the bacteria concentrations were monitored at beaches throughout

Massachusetts BayChlorination failure was modeled as a breakdown that

releases concentrations of 	





 cells

 ml for a period of  dayA

efolding decay timescale of  hours was used to simulate dieo of the

bacteria

Since currents are variable chlorination failures were simulated every 


days over a one year period from October  to October 
 and the

observed concentration levels were monitored at 	 dierent beaches and

harbors along the coast Table 	

Two runs were conducted one for the existing outfalls at Nut and Deer

Islands and one for the future outfall near the Boston BuoyFor each loca

tion the number of events where the concentration counts

 ml exceeded

various values were calculated and are shown in Tables 	 and 	Events


Number Longitude Latitude Location

 


 	 Glouchester Harbor

 

 
 Marblehead Neck

 

  Nahant Beach

 


  Winthrop Beach

	 




 
	
 Wollaston Beach

 

 


 Nantasket Beach


 

	

  NScituate Beach

 




  Humarock

 

 

 Duxbury Beach


 

 	

 Plymouth Beach

 

	
 


 Sagamore Beach

 

 
 Sandy Neck

 



 


 Great Island

 

 

 Race Point

Table 	 Locations monitored in chlorination failure simulations

were only counted if they occurred at least 	 days after the previous event

The future outfall location obviously has a large benecial eect as the

the total number of events with concentrations greater than 	 counts

 ml

dropping from  to  and the number of events greater than 
 counts



ml dropping from 
 to 
The existing outfall caused  events greater

than 


 counts

 ml whereas the highest concentration for the future

outfall was less than 
 counts

 ml



Counts

 ml Nahant Winthrop Wollaston Nantasket NScituate

	 	  	 
  


	 
 	  
  

	 	 
  
  

	 	
  
 
  

	 	 
  
  


	 	
 
  
 
 


	 	

 
    


	 


 
 
  
 


	 	

 
 
  
 


	 	


 
 
 
 
 


Table 	 Bacteria concentrations for breakdown of existing plantNote

Glouchester Harbor Marblehead Neck Humarock Duxbury Beach Ply

mouth Beach Sagamore Beach Sandy Neck Great Island and Race Point

had no events that exceeded 	 counts

 ml

Counts

 ml Nahant Winthrop Wollaston Nantasket NScituate

	 	  
 
  

	 
 
 
 
 
 


	 	 
 
 
 
 


	 	
 
 
 
 
 


	 	 
 
 
 
 


	 	
 
 
 
 
 


	 	

 
 
 
 
 


	 


 
 
 
 
 


	 	

 
 
 
 
 


	 	


 
 
 
 
 


Table 	 Bacteria concentrations for breakdown of future plant


Chapter 	

Summary

Comparison of model results with observations indicated that ECOMsi sim

ulations are appropriate for studying the continuous release of euent in

western Massachusetts Bay as they reproduce the most important features

of the observed stratication and current regime in this regionThe seasonal

development and breakdown of the thermocline occur within a few weeks of

the correct time the vertical structure has the right toptobottom density

dierence and about the right vertical prole though the modeled pycnocline

is not quite as sharp as reality and the currents responsible for transporting

material have about the right mean ow and level of variabilityThe fact

that the best comparisons are found on the western side of Massachusetts

Bay reects the strong eect of local wind driving in this region an eect

that can be directly modeled rather than parameterized by boundary condi

tionsFurther oshore near Stellwagen Bank for example the currents and

water properties are largely determined by Gulf of Maine conditions and

the model does not compare as wellAlthough signicant improvements to

salinity boundary conditions were obtained by including an adjunct model

to model the western Gulf of Maine boundary conditions remain one of the

largest roadblocks to obtaining more realistic simulationsIn the future

improvement of boundary condition information may be obtained from a

timevarying coarser scale model of the entire Gulf of Maine or possibly from

data collected from sensors at strategic locations and assimilated into model

runs

The euent simulations show that with either the existing or future out

fall the region of Massachusetts Bay that experiences relatively high euent


concentration  part euent 

 parts sea water is relatively small and

conned to western Massachusetts BayThe region of high concentration

	
 that covers much of Boston Harbor with the existing outfall is com

pletely eliminated with the future outfall scenarioThis is consistent with

the results of the twodimensional modeling conducted as part of the outfall

siting process MWRA The subsurface areal extent of the bottom

trapped summertime plume is about the same size as the winter plume for

the future outfall scenarioAlthough the initial vertical mixing of the euent

is conned to the lower half of column during the summer this eect is o

set by stronger currents and shears during the summer that act to increase

dispersionThe total amount of euent in the bay obtained from the fu

ture outfall is within 
 of that obtained from the existing outfall and the

time series track each other closelyEuent levels build up over the winter

reaching their peak in March before being greatly reduced by the increased

ushing of the system in April




Appendix A

ModelData Comparison of

Time Series Velocity at the

Boston Buoy



−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

 1 8 15 22 29 5 12 19 26

−10

0

10

cm/s

Jan Feb

1990

MODEL: 23 m

Figure A Boston Buoy velocity comparison Jan  
  Feb  




−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

26 5 12 19 26 2 9 16 23

−10

0

10

cm/s

Feb Mar Apr

1990

MODEL: 23 m

Figure A Boston Buoy velocity comparison Feb  
  Apr  




−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

23 30 7 14 21 28 4 11 18

−10

0

10

cm/s

Apr May Jun

1990

MODEL: 23 m

Figure A Boston Buoy velocity comparison Apr  
  Jun  




−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

18 25 2 9 16 23 30 6 13

−10

0

10

cm/s

Jun Jul Aug

1990

MODEL: 23 m

Figure A Boston Buoy velocity comparison Jun  
  Aug  



	
−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

13 20 27 3 10 17 24 1 8

−10

0

10

cm/s

Aug Sep Oct

1990

MODEL: 23 m

Figure A	 Boston Buoy velocity comparison Aug  
  Oct  




−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

 8 15 22 29 5 12 19 26 3

−10

0

10

cm/s

Oct Nov Dec

1990

MODEL: 23 m

Figure A Boston Buoy velocity comparison Oct  
  Dec  




−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

0

0.5

1

cm/s

DATA: 23 m

 3 10 17 24 31 7 14 21 28

−10

0

10

cm/s

Dec Jan

1990 1991

MODEL: 23 m

Figure A Boston Buoy velocity comparison Dec  
  Jan  




−4

−2

0

2

dynes/cm2

Wind Stress at BB

−10

0

10

cm/s

DATA: 5 m

−10

0

10

cm/s

MODEL: 5 m

−10

0

10

cm/s

DATA: 23 m

28 4 11 18 25 4 11 18 25

−10

0

10

cm/s

Jan Feb Mar

1991

MODEL: 23 m

Figure A Boston Buoy velocity comparison Jan    Mar