NOAA Technical Memorandum NOS CS 11

NOS Experimental Galveston Bay 

Nowcast/Forecast System: 

SF6 Tracer Studies 

Silver Spring, Maryland 

June 2007 

noaa National Oceanic and Atmospheric Administration 

U.S. DEPARTMENT OF COMMERCE 

National Ocean Service 

Coast Survey Development Laboratory 
Office of Coast Survey 

National Ocean Service 

National Oceanic and Atmospheric Administration 

U.S. Department of Commerce

The Office of Coast Survey (OCS) is the Nation’s only official chartmaker. As the oldest 

United States scientific organization, dating from 1807, this office has a long history. 

Today it promotes safe navigation by managing the National Oceanic and Atmospheric 

Administration’s (NOAA) nautical chart and oceanographic data collection and 

information programs. 

There are four components of OCS: 

The Coast Survey Development Laboratory develops new and efficient techniques 

to accomplish Coast Survey missions and to produce new and improved products 

and services for the maritime community and other coastal users. 

The Marine Chart Division acquires marine navigational data to construct and 

maintain nautical charts, Coast Pilots, and related marine products for the United 

States. 

The Hydrographic Surveys Division directs programs for ship and shore-based 

hydrographic survey units and conducts general hydrographic survey operations. 

The Navigational Services Division is the focal point for Coast Survey customer 

service activities, concentrating predominately on charting issues, fast-response 

hydrographic surveys, and Coast Pilot updates. 

ii
NOAA Technical Memorandum NOS CS 11 

NOS Experimental Galveston Bay 

Nowcast/Forecast System: 

SF6 Tracer Studies 

Richard A. Schmalz, Jr. 

Office of Coast Survey, Coast Survey Development Laboratory, Silver Spring, 

MD 

 

June 2007

noaa National Oceanic and Atmospheric Administration 

U. S. DEPARTMENT National Oceanic and National Ocean Service 

OF COMMERCE Atmospheric Administration John H. Dunnigan 

Carlos Gutierrez, Secretary Conrad C. Lautenbacher, Jr., Assistant Administrator 

 VADM USN (Ret.), Under Secretary 

 Office of Coast Survey Coast Survey Development Laboratory 

 Captain Steven R. Barnum, NOAA Mary Erickson 

iii
NOTICE 

Mention of a commercial company or product does not constitute an 

endorsement by NOAA. Use for publicity or advertising purposes of 

information from this publication concerning proprietary products or the 

tests of such products is not authorized. 

iv
 TABLE OF CONTENTS 

LIST OF FIGURES ....................................................................................................................... vi 

LIST OF TABLES....................................................................................................................... viii 

BASE MAP.................................................................................................................................... ix 

EXECUTIVE SUMMARY .............................................................................................................x 

1. INTRODUCTION .......................................................................................................................1 

2. LONG WAVE HYDRODYNAMIC MODEL............................................................................3 

2.1. Governing Equations ..........................................................................................................3 

2.2. Model Grids ........................................................................................................................6 

2.3. Air-Water Gas Transfer ......................................................................................................6 

3. SHORT WAVE HYDRODYNAMIC MODEL........................................................................11 

3.1. Governing Equations ........................................................................................................11 

3.2. Initial Validation ...............................................................................................................11 

3.3. Joint Computation.............................................................................................................14 

3.4. Refraction, Diffraction, and Reflection Limitations.........................................................14 

4. WAVE-CURRENT INTERACTION........................................................................................15 

4.1. Surface Drag Coefficient Adjustment...............................................................................15 

4.2. Bottom Friction Adjustment .............................................................................................15 

4.3. Surface Gas Transfer Velocity Adjustment ......................................................................16 

4.4. Wave Breaking Limitations.............................................................................................. 16 

5. NOVEMBER 2004 SF6 TRACER SIMULATION..................................................................19 

5.1. Set-Up ...............................................................................................................................20 

5.2. Long Wave Results...........................................................................................................20 

5.2.1. Water Surface Elevation.........................................................................................20 

5.2.2. Prediction Depth Currents......................................................................................20 

5.2.3. Temperature............................................................................................................21 

5.2.4. Salinity....................................................................................................................21 

5.3. Short Wave Results...........................................................................................................33 

5.3.1. Significant Wave Height ........................................................................................33 

5.3.2. Wave Direction.......................................................................................................33 

5.3.3. Wave Period ...........................................................................................................33

5.4. SF6 Tracer Results.............................................................................................................40 

5.4.1. Concentration Fields...............................................................................................40 

5.4.2. Residence Time ......................................................................................................45 

5.4.3. Exposure Level and Duration.................................................................................45 

5.4.4. Turnover Time........................................................................................................53 

5.4.5. Mass Inventory.......................................................................................................53 

6. CONCLUSIONS AND FUTURE WORK ................................................................................55 

ACKNOWLEDGMENTS .............................................................................................................55 

REFERENCES ..............................................................................................................................55 

v
LIST OF FIGURES 

Figure 1.1. Map showing Galveston Bay and the Houston Ship Channel 

 including PORTS stations...............................................................................................................2 

Figure 2.1.Galveston Bay model grid..............................................................................................8 

Figure 2.2. Galveston Bay model bathymetry, contours in meters..................................................9 

Figure 2.3. Houston Ship Channel model grid embedded in Bay model domain. ........................10 

Figure 3.1. Galveston Bay wave model versus observed wave parameters 

 at Eagle Point for January 25-30, 1997.......................................................................12 

Figure 3.2. Galveston Bay wave model versus observed wave parameters 

 at 42035 for January 25-30, 1997. ...............................................................................13 

Figure 5.1. Patrick Bayou Superfund site. Injection site is at confluence 

 of Patrick Bayou and the Houston Ship Channel. .......................................................19 

Figure 5.2. Freshwater Inflows for the Trinity River, San Jacinto River, and 

 Buffalo Bayou, November 1-12, 2004.........................................................................22 

Figure 5.3. Winds and sea level atmospheric pressure at Bolivar Roads, 

 November 1-12, 2004 ..................................................................................................23 

Figure 5.4. Winds and sea level atmospheric pressure at Eagle Point, 

 November 1-12, 2004.. ................................................................................................24 

Figure 5.5. Winds and sea level atmospheric pressure at Morgans Point, 

 November 1-12, 2004.. ................................................................................................25 

Figure 5.6. GBM simulated and observed water surface elevations 

 at Galveston Pleasure Pier, Bolivar Roads, and Galveston Pier 21, 

 November 1-12, 2004.. ................................................................................................26 

Figure 5.7. HSCM simulated and observed water surface elevations 

 at Eagle Point and Morgans Point and observed GBM water 

 level residuals, November 1-12, 2004.. .......................................................................27 

Figure 5.8. GBM simulated and observed prediction depth currents 

 at Bolivar Roads, November 1-12, 2004.. ..................................................................28 

Figure 5.9. HSCM simulated and observed prediction depth currents 

 at Morgans Point, November 1-12, 2004....................................................................29 

Figure 5.10. GBM simulated and observed temperature and surface salinity 

 at Bolivar Roads, November 1-12, 2004 ...................................................................30 

Figure 5.11. HSCM simulated and observed temperature and surface salinity 

 at Eagle Point, November 1-12, 2004........................................................................31 

Figure 5.12. HSCM simulated and observed temperature and surface salinity 

 at Morgans Point, November 1-12, 2004...................................................................32 

Figure 5.13. GBM and HSCM simulated significant wave height, 

 November 1, 2004......................................................................................................34 

Figure 5.14. GBM and HSCM simulated significant wave height, 

 November 12, 2004....................................................................................................35 

Figure 5.15. GBM and HSCM simulated significant wave direction, 

 November 1, 2004......................................................................................................36 

Figure 5.16. GBM and HSCM simulated significant wave direction, 

 November 12, 2004....................................................................................................37 

vi
vii

LIST OF FIGURES 

Figure 5.17. GBM and HSCM simulated significant wave period, 

 November 1, 2004......................................................................................................38 

Figure 5.18. GBM and HSCM simulated significant wave period, 

 November 12, 2004....................................................................................................39 

Figure 5.19. GBM simulated near surface SF6 concentration 

 at 00 CST on November 12, 2004 (Upper panel surface flux (CS1), 

 Lower Panel no surface flux (CS2)) ..........................................................................41 

Figure 5.20. GBM simulated near bottom SF6 concentration 

 at 00 CST on November 12, 2004 (Upper panel surface flux (CS1), 

 Lower Panel no surface flux (CS2)). .........................................................................42 

Figure 5.21. HSCM simulated near surface SF6 concentration 

 at 00 CST on November 12, 2004 (Upper panel surface flux (CS1), 

 Lower Panel no surface flux (CS2)) ..........................................................................43 

Figure 5.22. HSCM simulated near bottom SF6 concentration 

 at 00 CST on November 12, 2004 (Upper panel surface flux (CS1), 

 Lower Panel no surface flux (CS2)). .........................................................................44 

Figure 5.23. GBM simulated near surface residence time 

 (Upper panel surface flux (CS1), Lower Panel no surface flux (CS2)).....................46 

Figure 5.24. GBM simulated near bottom residence time 

 (Upper panel surface flux (CS1), Lower Panel no surface flux (CS2))......................47 

Figure 5.25. HSCM simulated near surface residence time 

 (Upper panel surface flux (CS1), Lower Panel no surface flux (CS2)).....................48 

Figure 5.26. HSCM simulated near bottom residence time 

 (Upper panel surface flux (CS1), Lower Panel no surface flux (CS2)).....................49 

Figure 5.27. Simulated near surface no surface flux (CS2) average exposure levels (fmol/L) 

 (Upper panel GBM, Lower Panel HSCM) ................................................................50 

Figure 5.28. Simulated near surface no surface flux (CS2) maximum exposure levels (fmol/L) 

 (Upper panel GBM, Lower Panel HSCM) ................................................................51 

Figure 5.29. Simulated near surface no surface flux (CS2) exposure duration (days) 

 (Upper panel GBM, Lower Panel HSCM) ................................................................52 
viii

LIST OF TABLES 

Table 5.1. SF6 Injection Characteristics at 00 CST on November 2, 2004. 

 Note CS1 represents SF6 with gas transfer, while CS2 

 represents SF6 with no gas transfer.............................................................................20 

Table 5.2. SF6 Turnover times (days). Note CS1 represents SF6 with gas transfer, 

 while CS2 represents SF6 with no gas transfer 

 and k=1 (surface),2,3,4, 5(bottom) sigma levels. ........................................................53 

Table 5.3. GBM SF6 Mass inventory (10-3 fmoles). Initial mass at 00 CST 

 on November 2, 2004 and final mass at 00CST on November 12, 2004. 

 Note CS1 represents SF6 with gas transfer, while CS2 

 represents SF6 with no gas transfer.............................................................................54 

Table 5.4. HSCM SF6 Mass inventory (10-3 fmoles). Initial mass at 00 CST 

 on November 2, 2004 and final mass at 00CST on November 12, 2004. 

 Note CS1 represents SF6 with gas transfer, while CS2 

 represents SF6 with no gas transfer.............................................................................54 
BUFFALO BAYOU 

SAN JACINTO RIVER TRINITY RIVER 

42035

Galveston Bay Base Map

ix
x

EXECUTIVE SUMMARY 

The Bay and Channel hydrodynamic models within the Galveston Bay Operational Forecast 

System (GBOFS) have been used to simulate the dispersion of a tracer in the Houston Ship 

Channel and upper Galveston Bay. With two concentration algorithms added to each model, the 

hydrodynamic model simulates the movement of the passive tracer sulfur hexafluoride (SF6) 

with and without surface gas transfer over a ten day period, November 2-11, 2004. To 

accomplish this, GBOFS has been extended to enable a hindcast on demand capability for up to a 

one month duration. 

The two hydrodynamic models have been set-up for simulating the tracer movement for a 

planned release of 1.0 mole of SF6 at the confluence of Patrick Bayou and the Houston Ship 

Channel. The Bay model is forced with USGS river inflows, PORTS observed surface winds and 

PORTS water levels at Galveston Pleasure Pier at the lateral open boundaries, which extend to 

the 20m isobath in the near Gulf. The Bay model is one-way coupled to the Channel model, 

which is also forced with USGS river inflows and PORTS observed surface winds. Discharges at 

the Trinity River, San Jacinto River, and Buffalo Bayou are considered as a one-dimensional 

inflow with a zero input of tracer. The simulated water levels, currents, and density are validated 

with PORTS observations to ensure model accuracy. Model simulated tracer concentration 

distributions and the total tracer mass balance are studied and were used to design the May 17 -

27, 2005 field experiment conducted by Columbia University. The areal extent of the 10 day 

residence time contours for both flux and no flux SF6 conditions were used to determine the 

necessary extent of boat survey coverage. While the actual survey was conducted in May 2005 

and the numerical simulations were for November 2004, the numerical and experimental 

residence time estimates were very similar and demonstrated the utility of the numerical models 

in planning actual tracer release field studies. 

In addition, surface gravity wave algorithms have been incorporated in both hydrodynamic 

models to simulate short period waves with and without wave-current interaction. Results 

presented indicate the feasibility of including the wave algorithms within GBOFS. 

Recommendations for additional SF6 simulation and wave-current interaction experiments are 

presented. 
1

1. INTRODUCTION

Recently attention has been focused on investigating the computation of residence time in 

estuaries and bays via multi-dimensional numerical modeling to aid in NOAA-EPA exposure 

assessment studies. In conjunction with a NOS Partnership Project, a workshop was held in June 

2004 (NOS, 2005) to further study the problem. One area of application was the Houston Ship 

Channel and the potential for toxic releases from the EPA Superfund site at Patrick Bayou. Wei 

(2004) has studied the transport of sulfur hexafluoride (SF6) within New York Harbor using the 

New York Harbor Operational Forecast System (NYOFS) model. Here, we attempt to build on 

this work and investigate the transport of SF6 within the Houston Ship Channel and Galveston 

Bay. 

Galveston Bay (Figure 1.1) has a complex geometry with the Houston Ship Channel bisecting 

the Bay and connecting the Port of Houston to the near Gulf. The Houston Ship Channel is 

important for both safe navigation and hydrodynamics in the Bay. Tidal currents through the 

channel play an important role in determining the dispersion characteristics. Flows from Buffalo 

Bayou, the San Jacinto River (Lake Houston Dam), and the Trinity River provide freshwater to 

Galveston Bay. These river inflows interact with the tidal currents to further complicate the 

circulation and transport in Galveston Bay. 

Columbia University researchers conducted a field experiment in the Houston Ship Channel to 

study the circulation, mixing, and the transport and the fate of solutes using SF6 during the period 

17-27 May 2005. Approximately 1.0 mol of SF6 was injected in the Houston Ship Channel at 

Patrick Bayou (downstream of a EPA superfund site) and the SF6 tracer was observed over the 

next 10 days using a high-resolution measurement system similar to that used in June 2003 in the 

East River, NY. Measured data were processed and compiled for dispersion characteristics 

interpretation (Schlosser et al., 2006). Detailed measurement system description and results from 

the NY study are described in Caplow et al. (2003). 

NOAA’s National Ocean Service (NOS) has developed the Galveston Bay Operational Forecast 

System (GBOFS) to simulate water levels, current velocities, and density for use by mariners 

navigating in Galveston Bay and the Houston Ship Channel (NOS, 1999; Schmalz, 1998b, 

1998c; Schmalz, 2000a, 2000b; Schmalz and Richardson, 2002) based on the Princeton Ocean 

Model (POM, Blumberg and Mellor, 1987). This forecast system (Schmalz, 2004) has been 

running operationally since June 2004 utilizing the near real-time water level and current 

information from NOS’ Physical Oceanographic Real-Time System (PORTS). The 

hydrodynamic models in this system have been used to simulate the SF6 transport in the Houston 

Ship Channel and upper Galveston Bay during the ten day period 2-11 November 2004 to 

examine the dispersion characteristics of the upper Houston Ship Channel. Areal extents of tracer 

concentration above background levels of 2 fmol/L were computed for the case of zero and nonzero SF6 surface fluxes and bracketed the areal extent of the measured SF6 concentrations above 

background level measured. (Refer to Figure 3 in Appendix 3 of Schlosser et al., 2006) and 

demonstrated the ability of numerical models to assist in the planning of the tracer release 

experiments in Galveston Bay. This report documents the modeling work performed to simulate 

the transport of SF6. The model set-up, simulation procedures, and results are described. 

Conclusions and future work based on the simulation results are discussed. 
Figure 1.1. Map showing Galveston Bay and the Houston Ship Channel including PORTS 

stations.

2
2. LONG WAVE HYDRODYNAMIC MODEL

2.1. Governing Equations

A three-dimensional sigma coordinate Galveston Bay and near shelf model (GBM) has been 

developed (Schmalz, 1996) based on a version of the Blumberg and Mellor (1987) model 

extended to orthogonal curvilinear coordinates. GBM is a three-dimensional baroclinic 

circulation model, based on POM, for simulating water levels, current velocities, and density. 

The model is forced with: water levels at the near shore open boundaries based on observed 

levels at Galveston Pleasure Pier; freshwater inflows from the Buffalo Bayou, San Jacinto, and 

Trinity Rivers; and surface winds. In addition, a high resolution Houston Ship Channel model 

(HSCM) has been incorporated to provide finer spatial resolution (Schmalz, 1998a; 2000a; 

2000c). The governing equations in a vertical sigma coordinate are briefly given as follows. 

Detailed formulation is contained in Blumberg and Mellor (1987), Mellor (2003b), and Schmalz 

(2001). 

=I 0)1( , (1) 

x

M Fd

x

gD D d

x

U gD

D

K

x

gDfDVU +

∂

∂

∂

∂

+

∂

∂ −⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = ∂∂ +−I ∫∫ σ

σ

ρ σ

ρ

σ ρ

ρσσ

η

σ σ

0

0

0

0

2

)( , (2) 

y

M Fd

y

gD D d

y

V gD

D

K

y

gDfDUV +

∂

∂

∂

∂

+

∂

∂ −⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = ∂∂ ++I ∫∫ σ

σ

ρ σ

ρ

σ ρ

ρσσ

η

σ σ

0

0

0

0

2

)( , (3) 

T

H F T

D

K

T +⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

σσ

)( , (4) 

S

H F S

D

K S +⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

σσ

)( , (5) 

q

q M H F

lB

VU gK Dq

D

q K

D

K

q +−

∂

∂

+

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡ ⎟

⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎟ + ⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

1

3

0

2 2 2

2 2 2 2 )(

σ

ρ

σσ ρσσ

, (6) 

H l

q M FW

B

Dq K VU gE

D

K lE lq

D

K

lq +−⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

∂

∂

+⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ⎟

⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎟ + ⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I ~

1

3

0

3

2 2

1

2

2 )( )(

σ

ρ

σσ ρσσ

, (7) 

where 

σ

ω

∂

∂

+

∂

∂

+

∂

∂

+

∂

∂ =I****

(*)

y

VD

x

UD

t

D and σ = (z-η)/(H+η). Here (x,y,z,t) are the Cartesian 

spatial and temporal coordinates and D=H+η is the total water depth with H the depth and η the 

water surface elevation with respect to model datum. U and V are horizontal velocities, S and T

are salinity and temperature, KM and KH are the vertical kinematic viscosity and diffusivity, 

respectively, Kq is vertical turbulence mixing coefficient, q

2

 is twice the turbulence kinetic 

energy, and l is the turbulence length scale. 

Note 2

2

~

)(1

L

l EW

κ

+= , κ = 0.4 is the von Karman constant with L-1=(η-z)-1+(H+z)-1, and BB

1, E1, 

E2, and E3 are constants. g is the acceleration due to gravity, f is the Coriolis parameter, ρ=f(S,T)

is the water density, and ω is the transformed vertical velocity normal to a sigma surface. The 

relation of ω ( σ / dtHd ) to the Cartesian vertical velocity w is 

3
tt

D

yy

D V

xx

D Uw

∂

∂

+

∂

∂

+⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

+

∂

∂

+⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

+

∂

∂

+= η σ η σ η σω . (8) 

The horizontal viscosity Fx and Fy and diffusion terms F* are defined as 

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+

∂

∂

∂

∂

+⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = x

V

y

U DA

yx

U DA

x

Fx 2 M M (9) 

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+

∂

∂

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = x

V

y

U DA

xy

V DA

y

Fy 2 M M , (10) 

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

∂

∂

+⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = y

DA

yx

DA

x

F H H

* * * , with *=(T,S,q,l,C) (11) 

where AM = AH are the horizontal eddy viscosity and diffusivity, defined by the Smagorinsky 

 formula (1963) 

 5.0 2 2 2

5.0 ⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+

∂

∂ ⎟ + ⎠

⎞ ⎜

⎝

⎛

∂

∂ ΔΔ=

y

V

x

V

y

U

x

U

NM yxCA , (12) 

with CN, a non-dimensional parameter. 

. 

For the passive tracer SF6, the concentration equation is 

C

H F C

D

K C +⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

σσ

)( (13) 

where C is the concentration of SF6.

Boundary conditions at the free surface (σ=0) are as follows: 

*3/2

1

2 ~~~

,0 ),,(, ,0),,,(,, s

H

sysx

M uBqlCTS CTS

H

VU K

H

K ⎟ = == ⎠

⎞ ⎜

⎝

⎛

∂

∂

∂

∂

∂

∂ ⎟ =

⎠

⎞ ⎜

⎝

⎛

∂

∂

∂

∂ = σσσ

ρττ

σσ

ρω (14) 

Boundary conditions at the bottom (σ=-1) are as follows: 

*3/2

1

2 ,0 ),,(, ,0),0,0,0(,, b

H

bybx

M uBql CTS

H

VU K

H

K ⎟ ===⎠

⎞ ⎜

⎝

⎛

∂

∂

∂

∂

∂

∂ ⎟ =

⎠

⎞ ⎜

⎝

⎛

∂

∂

∂

∂ = σσσ

ρττ

σσ

ρω (15) 

where τs and τb are the wind stress and bottom friction and ),,( are surface fluxes. ~~~

CTS

The above equation set is transformed into orthogonal horizontal coordinates after Blumberg and 

Herring (1987) with 

dt

d hU

dt

d hUdhdhdydxdsdhdydhdx 2

22

1

11

2

2

2

2

2

1

2

1

222

11 22 , , and , ε ε = ε = ε +=+= εε = = . 

=I 0)1( (16) 

12

0

11

0

0 11

1

112

2

~

1 )( Fd

h

D d

h

D U gD K

D

K

hh

gD U DU M

M f ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂

+

∂

∂ −⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = ∂

∂ +−I ∫ ∫ σ

σ

ρ σ

ε

σ

ε

ρ

ρσσε

η

σ σ

(17) 

4
21

0

22

0

0 22

2

221

1

~

2 )( Fd

h

D d

h

D U gD K

D

K

hh

gD U DU M

M f ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂

+

∂

∂ −⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ = ∂

∂ ++I ∫ ∫ σ

σ

ρ σ

ε

σ

ε

ρ

ε ρσσ

η

σ σ

(18) 

T

H F T

D

K

T +⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

σσ

)( (19) 

S

H F S

D

K S +⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

σσ

)( (20) 

q

q M H F

lB

U U gK Dq

D

q K

D

K

q +−

∂

∂

+

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡ ⎟

⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎟ +

⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

1

3

0

2

2

2

1

2

2 2 2 2 )(

σ

ρ

σσ ρσσ

 (21) 

H l

q M FW

B

Dq K U U gE

D

K lE lq

D

K

lq +−⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

∂

∂

+

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ⎟

⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎟ +

⎠

⎞ ⎜

⎝

⎛

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I ~

1

3

3

2

2

2

1

1

2

2 )( )(

σ

ρ

σσ σ ρσ

 (22) 

where 

σ

ω

ε ε ∂

∂

+

∂

∂

+

∂

∂

+

∂

∂ =I *)( *)(*)( *)( (*)

221

21

121

12

hh

DUh

hh

DUh

t

D

2

1

21

1

1

2

21

2

~

∂εε

∂ − ∂∂

+=

h

hh

Uh

hh

U f f

The relationship of ω with the Cartesian vertical velocity w is 

tt

D D

h

D U

h

U

w

∂

∂

+

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

+

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

+

∂

∂

+= η σ

ε

η

ε

σ

ε

η

ε

σω

1 22

2

2 11

1 . (23) 

The horizontal viscosity Fij and diffusion terms F* are defined as 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+

∂

∂

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =

ii

jj

jj

ii

M

ii j

ij

M

iji

ij h

Uh

h

Uh DA

h

Uh

DA

hh

F

ε εε εε

2 1 (24) 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

∂

∂ ⎥ +

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =

22

1

11 2

2

121

*

1 * *

ε εε h ε

h DA

h

h DA

hh

F H H , with *=(T,S,q,l,C) (25) 

where 

5.0 2

2

2

2

1

2

2

1

2

1

1

21 5.0 ⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+

∂

∂

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ == ε εεε

U UU U

NHM hhCAA , (26) 

where CN, a non-dimensional parameter, is set to be 0.005 for both Bay and Channel models. 

. 

For the passive tracer SF6, the concentration equation is 

C

H F C

D

K C +⎥

⎦

⎤ ⎢

⎣

⎡

∂

∂

∂

∂ =I

σσ

)( (27) 

where C is the concentration of SF6.

Boundary conditions at the free surface (σ=0) and at the bottom σ=-1) remain as given above, 

5
while the concentration lateral open boundary condition during the outflow is specified with onedimensional advection, ( ) ( ) = 0

∂

∂

+

∂

∂

ii

i

h

DCU

t

DC

ε

, where Ui is the velocity in the normal direction to 

the boundary while the following gradient condition is used for inflow, ( ) ( ) = 0

∂

∂ = ∂

∂

ii

i

h

DCU

t

DC

ε . 

2.2. Model Grids

The GBM computational grid as shown in Figure 2.1 consists of 181x101 horizontal cells (dx = 

254-2482m, dy= 580-3502m) with 5 levels in the vertical. GBM water depths range from 1 m in 

the shallows to 20m along the shelf boundary. The HSCM grid shown in Figure 2.3 was 

developed in three sections. Each grid section was linked in order to develop the final composite 

channel grid consisting of 71 x 211 horizontal cells (dx=63-1007m, dy=133-1268m) with the 

same 5 sigma levels as in the GBM. Note navigation channel depths are order 14m. The two 

models were then nested in a one-way coupling scheme, wherein GBM water surface elevation, 

salinity, temperature, turbulent kinetic energy, and turbulent length scale time histories were 

saved at 6-minute intervals to provide boundary conditions to drive the HSCM. For salinity, 

temperature, turbulent kinetic energy, and turbulent length scale, a one-dimensional (normal to 

the boundary) advection equation is used. On inflow GBM values are advected into the HSCM 

domain, while on outflow HSCM internal values are advected through the boundary. 

2.3. Air-Water Gas Transfer

After injection into the water column, part of the gaseous SF6 tracer exits from the water into the 

air. The gas transfer velocity, kSF6, which measures the air-water SF6 transfer rate, is determined 

based on k600, the transfer velocity relative to a Schmidt Number of 600, which is determined 

based on a wind speed term (0.296 coefficient after Caplow et al., 2003) plus a temperature 

enhancement term reported by Wanninkhof (1992) as follows: 

)109946.4106256.15246.0(5.2296.0 2 2 24

600 U10N Tx s Tx s k − − = + + + (28) 

 

Next the Schmidt number, Sc, is computed as a function of temperature (0 to 30 o

C) for salinity 

at 0 and 35 psu using the following relations given by Wanninkhof (1992): 

 

2 3

35

2 3

0

090558.02168.740.2316.3531

086070.08370.613.2173.3255

s s s

s s s

Sc T T T

Sc T T T

−= + −

−= + − (29) 

The final transfer velocity, kSF6, is then determined after Ho et al. (2002) as follows: 

n

s

SF

Sc

k k

⎟

⎠

⎞ ⎜

⎝

⎛

=

600

600

6 , where Scs, is the Schmidt number based on the surface salinity as determined 

from a linear interpolation from Sc0 and Sc35, and n is an exponent, which is set to 0.5 for no 

wave current interaction and is made a function of significant wave height, if wave current 

6
interaction is considered (Jahne et al., 1987). The surface flux of SF6, , is then determined by: 

~

C

)( 6

~

sSF −= csolCkC , (30) 

where Cs, is the surface SF6 concentration, and csol is the equilibrium solubility of SF6, here 

taken as 2 fmol/L, as reported by Caplow et al. (2003). 

7
Figure 2.1. Galveston Bay/Houston Ship Channel model grid showing locations of water 

level gages and current meter. 

8
Figure 2.2. Galveston Bay model bathymetry, contours in meter.

9
Figure 2.3. Galveston Bay model grid with Houston Ship Channel grid shown in green. 

10
3. SHORT WAVE HYDRODYNAMIC MODEL 

Schmalz (2003) has compared two USACE parametric wave models (CETN-I-6, 1981; Project 

CW-167, 1955) with the finite difference based Donelan (1977) wave model for wind events 

over Galveston Bay. The two parametric models considered were run on the same grid employed 

for the Bay circulation model (Figure 2.1), while a uniform square grid was used for the Donelan 

model. Further details on the Donelan (1977) wave model may be found in Schwab et al. (1984). 

Details of the initial testing may be found in Schmalz (2003). Best results as given in section 3.2 

below were achieved by the following mixture of the two parametric models as presented in the 

governing equations below. 

3.1. Governing Equations 

CETN-I-6 (1981) significant wave height, Hs (ft): 

]

]]/[530.0tanh[

]/[00565.0 tanh[]]/[530.0tanh[283.0 75.02

5.02

75.02

2

A

A

A

A

S

UgD

UgF UgD U

gH = (31)

CW-167 (Project CW-167, 1955): significant wave period, Ts (sec): 

4507.0

2 )(262.6

A

A

S U

gD

g

U T = (32) 

using the following notation: 

g = Gravitational acceleration (32.2 ft/sec2) 

UA = Windspeed (kts) 

F = Fetch (nm) 

D = Total water depth (ft) 

3.2. Initial Validation 

The January 25-30, 1997 time period was studied due to the availability of USACE wave 

measurements off Eagle Point. Application of the models to this period revealed that despite the 

incorporation of shallow water effects in the finite difference based Donelan model by using a 

linear reduction in transfer of wind to wave momentum as described in Schmalz (2003), best 

results were achieved by using a combination of the two simpler parametric wave models. Total 

significant wave height was computed as the sum of the CETN-I-6 results plus the swell, which 

was determined at the open boundary from the NDBC Buoy 42035 measurements. The total 

significant wave period was taken as equal to the CW-167 result alone. Results for January 25-

30, 1997 are shown in Figure 3.1 at Eagle Point and in Figure 3.2 at NDBC Buoy 42035, 

respectively. Note, the peak at Eagle Point is no longer delayed relative to the observations as 

experienced using the Donelan model with the shallow water adjustment and the peak at NDBC 

Buoy 42035 is well reproduced. 

11
 Figure 3.1. Galveston Bay Wave Model vs Observed Wave Parameters at Eagle Point, 

January 25-30 1997. 

12
 

Figure 3.2. Galveston Bay Wave Model vs Observed Wave Parameters at 42035, 

 January 25-30 1997.

13
3.3. Joint Computation 

The mixed parametric wave model has been included as a subroutine within both the Bay and 

Channel long wave models and uses the same wind field. The waves are computed every sixminutes. In the Bay model, fetch data are specified for each octant of wind direction at 25 

locations surrounding the Bay. An inverse distance squared interpolation is performed to 

determine the fetch distribution over the Bay for each wave computation. The long wave model 

provides the updated total water depth. At each wave computation in the Bay model, the swell 

height and fetch at each boundary point required by the Channel model are written on the 

transfer file in addition to the one-way coupled long wave information. Swell effects in the Bay 

model are input as a boundary condition and are reduced from the offshore boundary by a 

inverse distance squared interpolation of empirical reduction factors supplied at the above 25 

locations around the Bay. 

3.4. Refraction, Diffraction, and Reflection Limitations 

Wave refraction due to changes in bathymetry is not specifically treated. The wave field is 

represented here by a single frequency, height, and direction rather than a continuous spectrum. 

Thus no wave-wave interaction is considered. Wave diffraction around breakwaters and jetties 

and wave reflection are also not included. 

To include these three effects, a much more complicated and computationally intensive shallow 

water wave model must be included such as the Delft Technical University SWAN, which 

simulates the following physical phenomena: 1) wave generation and propagation in time and 

space, 2) shoaling, 3) refraction, 4) frequency shifting, 5) nonlinear wave-wave interactions, 6) 

whitecapping, and 7) blocking of waves by current. Note that diffraction and reflection are not 

explicitly modelled in SWAN but can be handled in the Delft Technical University PHAROS, a 

wave penetration finite element model. For further details refer to Rif (1997) for SWAN and 

Hurdle et al. (1989) for PHAROS.

As pointed out by Mellor and Donelan (2006), to perform shallow water wave computations with 

SWAN on the same grid as the POM circulation model a factor of 100 in computer time is 

needed. In practical applications, these computational requirements for shallow water wave 

computations are so severe that Mellor and Donelan (2006) advance a different initial approach. 

14
 4. WAVE-CURRENT INTERACTION MECHANICS 

We consider a wave field with significant wave height, H s, dominant period, C, and wave 

direction, φ. The wave direction is assumed to be in the direction of the wind. To determine the 

wave age, it is necessary to compute the wave phase speed, k

Cp

ω

= , where C

π

ω

2 = , and 

)tanh(

2

rkg

k ω = , with , where with 

D and g as defined previously as reported at 

⎪⎭

⎪

⎬

⎫

⎪⎩

⎪

⎨

⎧

++ ≤

−+ =

− −

2 031.0169.01(

2 )1221(

5.0 2

2 4

S SS S

eeS S

rk

S S f

/ gDS 2 = ω

http://web.mit.edu/fluidsmodules/www/potential_flows/Lectures. 

The friction velocity, u*, is then computed as = 10* CUu DN 10 , and the wave age is formed as the 

ratio of Cp to u*. 

4.1. Surface Drag Coefficient Adjustment 

In the presence of waves, the surface drag coefficient is increased. Following Drennan et al. 

(2003), we employ their Figure 10 in which the drag coefficient, CD10N, and wind speed, U10N, 

are given versus inverse wave age, u*/Cp. Wind speed ranges are from 5 to 20 m/s. The following 

relationships are used in the algorithm to define the adjustment factor of the surface drag, F, 

based the ratio of measured CD10N with respect to inverse wave age to the Smith (1980) 

relationship. 

N p

N p

N p

N p

F U Cu

F U Cu

F U Cu

F U Cu

/ 0.10 /15)5)(05.0(38.1

,10.0/ 0.08 15/)5)(07.0(27.1

,08.0/ 0.06 15/)5)(13.0(16.1

,06.0/ /15)5)(16.0(.1

10 *

10 *

10 *

10 *

p

p

p

−= −

−= − ≤

−= − ≤

−= − ≤

 (34) 

The effective surface drag coefficient, is given as the product of the adjustment factor, F, and the 

Large and Pond (1981) surface drag relation: 

eff ND

N N

N

ND

FCC

U smU

smU

xC

10

10 10

3 10

10 / 11 )065.049.0(

 2.1 / 11

10

=

⎭

⎬

⎫

⎩

⎨

⎧

+

≤ = f (35) 

Note with no wave current interaction, F=1.0. 

4.2. Bottom Friction Adjustment 

In the presence of waves, the near bottom wave orbital velocity based on linear wave theory, Uo

is first computed as )sinh(

5.0

kD

H

U s

o

ω = after Signell et al. (1990). The near bottom excursion 

amplitude, ab, and effective roughness, kb, are then determined as 0 30 and zk U

a b

o

b = = ω

, with z0

15
the bottom roughness. The wave friction velocity, u*w, is determined based on in the 

relationship reported by Grant and Madsen (1982) in the following manner: 

bbb = / akS

2

*

62.0

40.0

 with , 5.0

 23.0 0.1

23.0 0.08 0.1

13.0 08.0 

w ow

b

b b

b b

w Ufu

S

S S

S S

f =

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

= ≤≤

f

p

 (36) 

Next the current friction velocity, u*c, is determined based on the model horizontal velocity 

components (Ub,Vb), at the level nearest the bed as follows (see Davies and Lawrence, 1995): 

 with ,)]/1.0ln(/[2 )(5.0 22

*

2

c 0 c VUfu bbc f = κ z = + (37) 

The total friction velocity is then determined based on the current friction velocity and the wave 

friction velocity in the direction of the current as follows: 

( )2

*

2

* * )cos( T wc wc += uuu θ , where θ wc , is the angle between the wave and wind directions (refer 

to Grant and Madsen, 1979). A roughness adjustment is determined as T

c

u

u

b

T

p S

U

u

F *

* 1

0

* )( 24 −

= as 

reported by Signell (1990). The bottom friction adjustment factor, Fpb, is given 

by

2

0

0

)/ln(

)/ln(

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡ = zFh

zh F

pb

b

pb . In the present study, Fpb cannot exceed 2. The effective bottom 

roughness is the product of the bottom roughness and this adjustment factor. Note in this 

approach, wave and current effects are considered independently (z0 is never altered) and are 

then combined to determine the adjustment factor. In theory, an iterative approach on z0 is 

desired. 

4.3. Surface Gas Transfer Velocity Adjustment 

Jahne et al. (1987) report that in the presence of waves on the surface, two important 

mechanisms are exhibited. First, the hydrodynamic boundary layer changes such that local 

divergences and convergences occur; this results in the exponent of the Schmidt number, n, 

decrease from 2/3 to 1⁄2. Second, additional energy is cycled with the energy gained by the waves 

being transferred to the near-surface turbulence. In our approach, the following relation is used to 

determine the Schmidt number exponent based on significant wave height. 

s

s

s

s

n H

n H

n H

n H

 5.0 50.0

5.0 4.0 55.0

4.0 3.0 60.0

3.0 667.0

p

p

p

=

= ≤

= ≤

= ≤

 (38) 

Note, for the case of no wave-current interaction, we still set n=0.50. 

4.4. Breaking Wave Limitation 

16
Within the context of the wave-current interaction algorithms considered above, wave breaking 

effects are not explicitly considered. Under the present approach, wave-current interactions act to 

increase surface and bottom stress, which leads to an increase in turbulence and enhanced 

vertical mixing. Mellor and Blumberg (2004) have considered wave breaking effects on the sea 

surface temperature. It would appear that additional research is needed on shallow water wave 

breaking effects and ultimately a shallow water wave and circulation model must be coupled. 

See Mellor and Donelan (2006) for one approach in this regard. 

17
18
5. NOVEMBER 2004 SF6 TRACER SIMULATION 

A simulated release of SF6 was conducted from November 2 to November 11, 2004. One mole 

was injected and assumed to be uniformly dissolved over the vertical at 00 CST on 2 November 

at the confluence of Patrick Bayou and the Houston Ship Channel (Figure 5.1). Note during the 

May 2005 experiment conducted by Columbia University approximately 1 mole of SF6 was also 

injected at the same location as indicated in Appendix 2 of Schlosser et al. (2006). After the 

injection, the SF6 tracer was then simulated over the next 10 days under the assumption of 

surface (CS1) and no surface gas (CS2) transfer. This was accomplished by using two separate 

concentration algorithms for CS1 and CS2. Detailed injection characteristics are given in Table 

5.1. The Bay model long and short wave results were used to provide the boundary conditions 

for the finer resolution Houston Ship Channel simulation. Simulations were performed initially 

with no wave current interaction and are described here. Simulation results with wave-current 

interaction yielded essentially the same results, with a change in water surface elevations of 

order 0.5 cm. 

Figure 5.1. Patrick Bayou Superfund site. Injection site is at the confluence of Patrick Bayou and 

the Houston Ship Channel 

19
Table 5.1. SF6 Injection Characteristics at 00 CST on November 2, 2004. Note CS1 represents 

SF6 with surface gas transfer, while CS2 represents SF6 with no surface gas transfer. 

Parameter Galveston Bay Model Houston Ship Channel Model 

Total Mass (moles) 1.0 1.0 

(I,J) Grid Location (50,88) (31,178) 

(DX,DY) grid spacings (m) (378.0, 1043.9) (213.8, 786.3) 

Water depth with respect to 

model datum (m) 

5.2 11.2 

Total water depth (m) 5.5347 11.5452 

Volume (m3) 2,183,710 1,941,180 

CS1 (k=1,5) (fmol/L) 457937 515151 

CS2 (k=1,5) (fmol/L) 457937 515151 

5.1. Set-Up

Observed river discharges for the Trinity River, San Jacinto, and Buffalo Bayou are shown in 

Figure 5.2 and exhibit relatively low flow conditions with a peak total flow of order 10,000 cfs 

on the 4 and 5 of November. Winds and atmospheric pressure as reported at the PORTS met 

stations shown at Bolivar Roads, Eagle Point, and Morgans Point in Figures 5.3-5.5, 

respectively, are uniform over the Bay with wind speeds order 10 to 12 kts. Water level residuals 

at Galveston Pleasure Pier are less than 25 cm, indicating a rather tranquil period. 

The models were set up to simulate the SF6 tracer concentration and to study the dispersion 

characteristics in the Houston Ship Channel and Upper Galveston Bay. The models simulated the 

water levels, currents, salinity, temperature, and SF6 concentration with and without surface gas 

exchange from November 2 to November 12, 2004. The model was spun-up for one day from 

rest before the tracer injection at 00 UTC, November 2. 

5.2. Long Wave Results 

5.2.1. Water Surface Elevation 

Simulated and observed water levels at Galveston Pleasure Pier, Bolivar Roads, and Galveston 

Pier 21, shown in Figure 5.6, indicate that the Bay model is accurately reproducing the water 

elevations at model interior locations. Water levels at Eagle Point and Morgans Point in Figure 

5.7 indicate that the Channel model is also accurately reproducing the water elevations at model 

interior locations. Note the set-down of order 25cm during the period 3-8 November. In the first 

half of this period, the winds are from the North, and the Bay tends to empty with levels below 

the predicted astronomical tide. 

5.2.2. Prediction Depth Currents 

Simulated prediction depth currents at Bolivar Roads (Figure 5.8) from the Bay model and at 

Morgans Point (Figure 5.9) from the Channel model are in general agreement with the data. Note 

20
the currents are bi-directional with strengths reduced as one proceeds from Bolivar Roads (125 

cm/s) to Morgans Point (25-30 cm/s). Both models tend to slightly underestimate the peak 

currents with errors in direction of order 25 degrees. 

5.2.3. Temperature

Temperature at Bolivar Roads (Figure 5.10) from the Bay model and at Eagle Point (Figure 5.11) 

and Morgans Point (Figure 5.12) from the Channel model are initialized to PORTS data. The 

initial SST field and boundary conditions are persisted over the 11 day period. This results in an 

error of order 5 o

C by the end of the simulation period in sea surface temperature. Stratification 

at Bolivar Roads is underestimated during the 3-5 November period. A refinement to the set-up 

procedure would allow a linear increase or decrease of the initial climatological offshore 

boundary condition to be applied over the simulation. 

5.2.4. Salinity

Surface salinity at Bolivar Roads (Figure 5.10) from the Bay model and at Eagle Point (Figure 

5.11) and Morgans Point (Figure 5.12) from the Channel model are also initialized to PORTS 

data. No daily data assimilation of the PORTS salinity are performed. A persistence of the 

climatological boundary condition is used. Note in Figure 5.10 at Bolivar Roads, the presence of 

many lows in the observed time series, due to the advection of large gradients. Without data 

assimilation the errors in surface salinity are order 5 psu at the end of the 11 day simulation 

period. The above linear adjustment to the climatological offshore boundary condition would 

improve the salinity results as well 

21
Figure 5.2. Freshwater Inflows for the Trinity River, San Jacinto River, and Buffalo Bayou, 

November 1-12, 2004. 

22
Figure 5.3. GBM winds and sea level atmospheric pressure at Bolivar Roads, 

November 1-12, 2004. 

23
Figure 5.4. HSCM winds and sea level atmospheric pressure at Eagle Point, 

November 1-12, 2004. 

24
Figure 5.5. HSCM winds and sea level atmospheric pressure at Morgans Point, 

November 1-12, 2004.

25
Figure 5.6. GBM simulated and observed water surface elevations at Galveston Pleasure Pier, 

Bolivar Roads, and Galveston Pier 21, November 1-12, 2004. 

26
Figure 5.7. HSCM simulated and observed water surface elevations at Eagle Point and Morgans 

Point and observed GBM water level residuals, November 1-12, 2004. 

27
Figure 5.8. GBM simulated and observed prediction depth currents at Bolivar Roads, 

November 1-12, 2004. 

28
Figure 5.9. HSCM simulated and observed prediction depth currents at Morgans Point, 

November 1-12, 2004. 

29
Figure 5.10. GBM simulated and observed temperature and surface salinity at Bolivar Roads, 

November 1-12, 20004. 

30
Figure 5.11. HSCM simulated and observed temperature and surface salinity at Eagle Point, 

November 1-12, 2004. 

31
Figure 5.12. HSCM simulated and observed temperature and surface salinity at Morgans Point, 

November 1-12, 2004. 

32
5.3. Short Wave Results 

5.3.1. Significant Wave Height 

Significant wave heights from the Galveston Bay model shown for 1 November 2004 in the top 

panel of Figure 5.13 are order 0.9m in the near shelf regions and are reduced to order 0.5m in the 

lower and upper Galveston Bay and West Bay regions. In the East Bay and Trinity Bay wave 

heights are order 0.3m. In the lower panel of Figure 5.13 significant wave heights from the 

Houston Ship Channel model agree in magnitude with those of the Bay model. 

Significant wave heights from the Galveston Bay model shown for 12 November 2004 in the top 

panel Figure 5.14 are less than 0.5m throughout all Bay regions and on the near shelf. In the 

lower panel of Figure 5.14 significant wave heights from the Houston Ship Channel model agree 

in magnitude with those of the Bay model. 

5.3.2. Wave Direction 

Wave directions from the Galveston Bay model shown for 1 November 2004 in the top panel of 

Figure 5.15 are to the North consistent with the wind directions shown in Figures 5.3-5.5. Wave 

directions are equal to the wind directions as obtained from a two-step Barnes interpolation from 

the PORTS stations. The Barnes interpolation procedure does not guarantee strict equality of the 

interpolated and injected data point. Wave directions from the Houston Ship Channel model in 

the lower panel of Figure 5.15 are in general agreement with those of the Bay model, but may be 

slightly different due to the fact that a separate Barnes interpolation is performed over the 

Channel model domain. 

Wave directions from the Galveston Bay model shown for 12 November 2004 in the top panel of 

Figure 5.16 are from the North consistent with the wind directions shown in Figures 5.3-5.5. 

Again the same remarks regarding the two-step Barnes interpolation from the PORTS stations 

apply. Wave directions from the Houston Ship Channel model in the lower panel of Figure 5.16 

are in general agreement with those of the Bay model, but may be slightly different due to the 

separate Barnes interpolation. 

5.3.3. Wave Period 

Wave periods from the Galveston Bay model shown for 1 November 2004 in the top panel of 

Figure 5.17 are order 3-4 sec, with 5 sec waves found in the near shelf regions and in the deeper 

navigation channel. In the lower panel of Figure 5.17 wave periods from the Houston Ship 

Channel model agree in magnitude with those of the Bay model with longer period waves being 

again found in the navigation channel. 

Wave periods from the Galveston Bay model shown for 12 November 2004 in the top panel 

Figure 5.18 are of the same order as found on 1 November. In the lower panel of Figure 5.18 

wave periods from the Houston Ship Channel model agree in magnitude with those of the Bay 

model. 

33
Figure 5.13. GBM and HSCM simulated significant wave height, November 1, 2004. 

34
Figure 5.14. GBM and HSCM simulated significant wave height, November 12, 2004 

35
Figure 5.15. GBM and HSCM simulated wave direction, November 1, 2004. 

36
Figure 5.16. GBM and HSCM simulated wave direction, November 12, 2004. 

37
Figure 5.17. GBM and HSCM simulated wave period, November 1, 2004.

38
Figure 5.18. GBM and HSCM simulated wave period, November 12, 2004.

39
5.4. SF6 Tracer Results

5.4.1. Concentration Fields

Near surface Galveston Bay Model SF6 concentration fields 10 days after the injection, on 00 

CST on 12 November, 2004 are shown in Figure 5.19. In the top panel, results are shown for the 

case of surface gas transfer, while in the bottom panel results are shown for no surface gas 

transfer. For the no surface flux case, the surface concentrations are elevated from those for the 

surface flux case and present in the lower portions of the Port of Houston at Morgans Point. 

Near bottom Galveston Bay Model SF6 concentration fields are given in Figure 5.20. In the top 

panel, results are shown for the case of surface gas transfer, while in the bottom panel results are 

shown for no surface gas transfer. Note that in the surface flux CS1 fields, the maximum legend 

value is order 103

 and the concentrations are elevated at the bottom (from the surface, Figure 

5.19 top panel). In the no surface flux CS2 fields, the maximum legend value is 104 and the 

fields are uniformly mixed over the vertical. 

Corresponding results for the Houston Ship Channel Model SF6 concentration fields are shown 

at the near surface in Figure 5.21 and at the near bottom in Figure 5.22. The SF6 tracer reaches 

the boundary and therefore there is some uncertainty as to the influence of the zero flux inflow 

condition on the results. The southern extent of the tracer in both surface and bottom figures is 

reduced from the Bay model results. 

40
Figure 5.19. GBM simulated near surface SF6 Concentrations at 00 CST on 12 November, 2004. 

(Upper panel surface flux (CS1), Lower panel no surface flux (CS2)).

41
 Figure 5.20. GBM simulated near bottom SF6 Concentrations at 00 CST on 12 November, 2004. 

(Upper panel surface flux (CS1), Lower panel no surface flux (CS2)).

42
Figure 5.21. HSCM simulated near surface SF6 Concentrations at 00 CST on 12 November, 

2004. (Upper panel surface flux (CS1), Lower panel no surface flux (CS2)).

43
Figure 5.22. HSCM simulated near bottom SF6 Concentrations at 00 CST on 12 November, 

2004. (Upper panel surface flux (CS1), Lower panel no surface flux (CS2)).

44
5.4.2. Residence Time

The time to exposure for SF6 is defined as the time since the start of the injection when the 

concentration at the given point exceeded a threshold value of 2 fmol/L. The residence time is 

defined as the difference in the time to exposure of the given point and the injection point. Since 

at the injection time, the initial concentration of SF6 equaled 457937 fmol/L in all 5 vertical 

levels in the GBM and equaled 515151 fmol/L in all 5 vertical levels in the HSCM, which 

exceeded the threshold value, the time to exposure and residence times are equivalent. 

 In Figure 5.23, the GBM near surface residence times are given for CS1 surface flux and CS2 no 

surface flux. In the case of no surface flux, the extent of the nonzero residence time is greatly 

enhanced. This may be used by the field research team to plan the area of coverage over a ten 

day experiment. In Figure 5.24, the GBM near bottom residence times are given for CS1 and 

CS2, with the areal coverage of nonzero residence time increased for the no surface flux 

condition. The area of nonzero residence time at the surface for CS1 is similar to the bottom 

area. For CS2 the surface and bottom areas are quite different. 

In Figures 5.25 and 5.26, the corresponding HSCM residence time results are shown. Note that 

the southern extent of the nonzero residence times is reduced from that shown in the GBM. The 

extent of the nonzero residence time, however, does reach the lateral boundaries and this may be 

a source of error, since a zero gradient condition is used on inflow. 

5.4.3. Exposure Level and Duration

To further assess the exposure, the average exposure level and maximum exposure levels are 

determined as well as the duration of the exposure. Note the exposure level is determined as the 

concentration of SF6 exceeding the threshold value of 2 fmol/L. While these results are available 

at each of the 5 vertical layers, we examine here only the first near surface layer. The near 

surface average and maximum exposure levels are shown in Figure 5.27 and Figure 5.28, 

respectively, for both of the models. In Figure 5.29, the duration of exposure is presented. These 

type of computations should enable direct exposure assessments. 

45
Figure 5.23. GBM near surface residence time. (Upper panel surface flux (CS1), Lower panel no 

surface flux (CS2)).

46
Figure 5.24. GBM simulated near bottom residence time. (Upper panel surface flux (CS1), 

Lower panel no surface flux (CS2)).

47
Figure 5.25. HSCM simulated near surface residence time. (Upper panel surface flux (CS1), 

Lower panel no surface flux (CS2)). 

48
Figure 5.26. HSCM simulated near bottom residence time. (Upper panel surface flux (CS1), 

Lower panel no surface flux (CS2)).

49
Figure 5.27. Simulated near surface no surface flux (CS2) average exposure levels (fmol/L). 

(Upper panel GBM, Lower panel HSCM). 

50
Figure 5.28. Simulated near surface no surface flux (CS2) maximum exposure levels (fmol/L). 

(Upper panel GBM, Lower panel HSCM).

51
Figure 5.29. Simulated near surface no surface flux (CS2) exposure duration (days). (Upper 

panel GBM, Lower panel HSCM).

52
5.4.4. Turnover Time 

With respect to the initial injection grid cell, the total turnover time is determined as the time at 

which the total mass is reduced by 1/e of the initial mass. The turnover times are given in Table 

5.2 for both Bay and Channel models for each of the 5 levels as well. Note it is possible, that 

once the mass is reduced to 1/e of its initial mass, the mass may increase again. This is not 

indicated in the turnover time and must be determined by examining the concentration time 

history at the given location. 

Table 5.2. SF6 Turnover times (days). Note CS1 represents SF6 with surface gas transfer, while 

CS2 represents SF6 with no surface gas transfer and k=1 (surface),2,3,4, and 5 (bottom) sigma 

levels. 

Parameter Galveston Bay Model Houston Ship Channel Model 

Total (CS1,CS2) (0.1674, 0.1708) (0.1019,0.1030) 

(CS1,CS2,k=1) (0.1965, 0.2056) (0.2024,0.2111) 

(CS1,CS2,k=2) (0.0701, 0.0750) (0.0378,0.0378) 

(CS1,CS2,k=3) (0.1097,0.1486) (0.0497,0.0496) 

(CS1,CS2,k=4) (0.1458,0.1492) (0.0766,0.0766) 

(CS1,CS2,k=5) (0.1514,0.1555) (0.1297,0.1300) 

5.4.5. Mass Inventory 

To check on the mass consistency of the computations, a mass inventory was determined at the 

end of the simulation. The difference in the initial mass and the sum of the final mass, mass lost 

through the surface, and the mass lost through the lateral boundaries was determined as the mass 

error. It should be noted that while the mass balance computations were carried out using double 

precision, the quantities making up the balance were computed based on single precision 

arithmetic. The mass relative error was determined as the ratio of the mass error to the initial 

mass and was order 10-7, which represent the precision of the single precision arithmetic. Results 

are given in Table 5.3 for the Galveston Bay Model and in Table 5.4 for the Houston Ship 

Channel model computations. It is interesting to observe the difference in the results for the Bay 

and Channel models, particularly with respect to final mass. Note the negative mass losses 

through the surface for CS2 in both model computations are not theoretically possible, but are 

due to the single precision arithmetic. Also note that the initial mass should show no nonzero 

numbers after the leading 1, and again shows the limits of the single precision arithmetic. Decay 

coefficients are computed based on an exponential decay from initial to final mass plus surface 

mass loss for CS1 and from initial to final mass for CS2. The decay coefficients for the 

Galveston Bay Model are nearly equal, while there is an order of 3 difference for the Houston 

Ship Channel Model. 

53
Table 5.3. Galveston Bay Model: SF6 Mass inventory (10-3 fmoles). Initial mass at 00CST on 2 

November 2004 and final mass 00CST on 12 November 2004. Note CS1 represents SF6 with 

surface gas transfer, while CS2 represents SF6 with no surface gas transfers. Note the mass 

balance region considered is I = (15,76), J=(79,95). 

Parameter CS1 CS2 

Initial Mass 1000000034488.555 1000000034488.555 

Final Mass 55303991283.367 994803495802.576 

Mass lost through the surface 941127766951.415 -430414.673 

Mass lost through the lateral 

boundaries 

 3568325477.928 55197231770.618 

Mass error -49224.155 -262669.966 

Mass relative error to the 

initial mass (-) 

-4.922x10-8 -2.627x10-7

Decay coefficient (days-1) -0.00036 -0.00052 

Table 5.4. Houston Ship Channel Model: SF6 Mass inventory (10-3 fmoles). Initial mass at 

00CST on 2 November 2004 and final mass 00CST on 12 November 2004. Note CS1 represents 

SF6 with surface gas transfer, while CS2 represents SF6 with no surface gas transfers. Note the 

mass balance region considered is I = (9,64), J=(174,205). 

Parameter CS1 CS2 

Initial Mass 1000000034488.555 1000000034488.555 

Final Mass 156231572263.406 759593945725.248 

Mass lost through the surface 774541589083.615 -945999.427 

Mass lost through the lateral 

boundaries 

 69227572621.036 240406536065.383 

Mass error -699479.503 498697.350 

Mass relative error to the 

initial mass (-) 

-6.994x10-7 4.987x10-7

Decay coefficient (days-1) -0.0073 -0.0274 

54
6. CONCLUSIONS AND FUTURE WORK

The estuarine dispersion characteristics of the Houston Ship Channel and upper Galveston Bay 

have been investigated by the simulated movement of SF6 tracer at Patrick Bayou in November 

2004. With the addition of two concentration subroutines (with and without surface gas transfer) 

to the Bay and Channel hydrodynamic models and forcing with proper boundary conditions, the 

models have been set up to simulate the tracer movement during 2-11 November 2004. 

Computed areal extents of the 10 day residence time contours for both zero and nonzero surface 

fluxes of SF6 bracketed the measured SF6 contours above background concentration and 

demonstrated the utility of the numerical models in planning field experiments, which will 

determine residence and flushing times directly from measured tracer concentration data. 

The concentration subroutine approach is an effective tool for characterizing the dispersion 

features of a soluble substance in an estuary. For other pollutants particle tracking may be a more 

appropriate approach. This report only documents the model simulation of the SF6 tracer 

dispersion using the concentration approach. Lagrangian trajectory modeling, as reported by Wei 

(1994) and Blumberg et al. (2004), could also be applied to study the SF6 tracer dispersion 

experiments and compared with results from the concentration modeling. In addition, the concept 

of age, introduced by Deleersnijder et al. (2001) as applied by Shen and Hass (2004), could also 

be investigated to further characterize dispersion characteristics. The residence time estimates as 

studied by Wang et al. (2004) might also be tried as alternative measures as well. 

Additional wave-current interaction experiments, should be conducted in which the limit on the 

bottom adjustment factor of 2.0 is altered. Additional work in this area, by Mellor (2003a), 

Johnson (1992), and Monahan (2002) should also be considered in an effort to develop an 

iterative approach for z0. It would also be of interest to compare the SWAN and simpler 

parametric model results to further quantify their limitations. 

ACKNOWLEDGMENTS

The solute model work has been carried out in the Coast Survey Development Laboratory 

(CSDL) of NOS’s Office of Coast Survey with the encouragement from the Chief of CSDL’s 

Marine Modeling and Analysis Programs, Dr. Frank Aikman. Discussions with Dr. Eugene Wei 

regarding his work on the simulation of the SF6 tracer in the New York Harbor have greatly 

assisted this work.

REFERENCES 

Blumberg, A. F., and G. L. Mellor, 1987. A description of a three-dimensional coastal ocean 

circulation model. Three-Dimensional Coastal Ocean Models, (ed. Heaps), American 

Geophysical Union, Washington, DC., 1 - 16. 

Blumberg, A. F., and H.J. Herring, 1987. Circulation modeling using orthogonal curvilinear 

coordinates, in Three-Dimensional Models of Marine and Estuarine Dynamics, (ed. Nihoul and 

Jamart), Elsevier Oceanography Series, 45, 55-88. 

Blumberg, A. F., D.J. Dunning, H. Li, S. Heimbuch, and W.R. Geyer, 2004. Use of a Particletracking Model for Predicting Entrainment at Power Plants on the Hudson River, Estuaries, Vol. 

55
27, No. 3, 515-526. 

Caplow, T., P. Schlosser, D.T. Ho, and N. Santella, 2003. Transport Dynamics in a Sheltered 

Estuary and Connecting Tidal Straits: SF6 Tracer Study in New York Harbor, Environmental 

Science Technology, 37, 5116-5126. 

Coastal Engineering Technical Note, CETN-I-6, 1981. Revised method for wave forecasting in 

shallow water, Coastal Engineering Research Center, Vicksburg, MS. 

Davies, A.M. and J. Lawrence, 1995. Modeling the effect of wave-current interaction on the 

three-dimensional wind-driven circulation of the Eastern Irish Sea, Journal of Physical 

Oceanography, 25, 29-45. 

Deleersnijder, E., J-M Campin, E.J.M. Delhez, 2001. The concept of age in marine modeling I. 

Theory and preliminary model results, Journal of Marine Systems, 28, 229-267. 

Donelan, M.A., 1977: A simple numerical model for wave and wind stress prediction, Report, 

National Water Resources Institute, Burlington, Ontario. 

Drennan, W.M, H.C. Graber, D. Hauser, and C. Quentin, 2003. On the wave age dependence of 

wind stress over pure wind seas, Journal of Geophysical Research, 108, 8062. 

Grant, W.D and O.S. Madsen, 1979. Combined wave and current interaction with a rough 

bottom, Journal of Geophysical Research, 84, 1797-1808. 

Grant, W.D. and O.S. Madsen, 1982. Movable bed roughness in unsteady oscillatory flow, 

Journal of Geophysical Research, 87, 469-481. 

Ho, D.T., P. Schlosser, and T. Caplow, 2002. Determination of Longitudinal Dispersion 

Coefficient and Net Advection in the Tidal Hudson River with a Large-Scale, High Resolution 

SF6 Tracer Release Experiment, Environmental Science Technology, 36, 3232-3241. 

Hurdle, J.P., J.K. Kostense, and P.van den Bosch, 1989. Mild-slope model for the wave 

behaviour in and around harbours and coastal structures in areas of variable depth and flow 

conditions, In: Advances in Water Modelling and Measurement. 

Jahne, B., K.O. Munnich, R. B. Singer, A. Dutzi, W. Huber, and P. Libner, 1987. On the 

Parameters Influencing Air-Water Gas Exchange, Journal of Geophysical Research, Vol. 92, C2, 

1937-1949. 

Johnson, H.K., and H.J. Vested, 1992. Effects of water waves on wind shear stress for current 

modeling, Journal of Atmospheric and Oceanic Technology, 9, 850-861. 

Large, W.G. and Pond, S., 1981. Open ocean momentum flux measurements in moderate to 

strong winds, Journal of Physical Oceanography, 11, 324-336. 

56
Mellor, G.L., 2003a. The three-dimensional, current and surface wave equation, Journal of 

Physical Oceanography, 33, 1978-1989. 

Mellor, G.L, 2003b.Users guide for a three-dimensional, primitive equation, numerical ocean 

model, Progr Atmos and Ocean Sci, Princeton University. (Available at 

www.aos.princeton.edu/WWWPUBLIC/htdocs.pom/publications.htm) 

Mellor, G. L. and A. F. Blumberg, 2004. Wave breaking and ocean surface layer thermal 

response, J. Phys. Oceanogr., 34, 693-698. 

Mellor, G. L. and M. A. Donelan, 2006. A surface wave model for coupling with numerical 

ocean circulation models, submitted. 

Monahan E. C., 2002. Oceanic whitecaps: Sea surface features detectable via satellite that are 

indicators of the magnitude of the air-sea gas transfer coefficient, Proceedings Indian Academy 

of Sciences, Earth and Planetary Science, 111, No 3, 315-319. 

NOS, 1999. NOS procedures for developing and implementing operational nowcast and forecast 

systems for PORTS, NOAA Technical Report NOS CO-OPS 0020, Silver Spring, MD. 

NOS, 2005. Report on the National Ocean Service Workshop on Residence Time in Bays and 

Estuaries, NOS, CSDL, Silver Spring, MD.

Project CW-167, 1955. Waves and wind tides in shallow lakes and reservoirs: summary report, 

U.S. Army Corps of Engineers, Jacksonville, FL. 

Ris, R.C., 1997: Spectral modelling of wind waves in coastal areas, (Ph.D. Dissertation Delft 

University of Technology, Department of Civil Engineering), Communications on Hydraulic and 

Geotechnical Engineering, Report No. 97-4, Delft, The Netherlands. 

Schlosser, P., D. Ho, and P. Schmieder, SF6 Tracer Release Experiment in the Houston Ship 

Channel, Final Report, Lamont-Doherty Earth Observatory of Columbia University, 

New York, NY, June 2006. 

Schmalz, R.A., 1996. National Ocean Service Partnership: DGPS-supported hydrosurvey, water 

level measurement, and modeling of Galveston Bay: development and application of the 

numerical circulation model. NOAA, National Ocean Service, Office of Ocean and Earth 

Sciences, NOAA Technical Report NOS OES 012, Silver Spring, MD. 

Schmalz, R.A., 1998a. Development of a Nowcast/Forecast System for Galveston Bay, 

Proceedings of the 5th Estuarine and Coastal Modeling Conference, ASCE, Alexandria, VA, 

441-455. 

Schmalz, R.A., 1998b. Design of a Nowcast/Forecast System for Galveston Bay, Proceedings of 

the 2nd Conference on Coastal Atmospheric and Oceanic Prediction and Process, AMS, 

Phoenix, AZ, 15-22. 

Schmalz, R.A., 1998c. Initial Evaluation of a Nowcast/Forecast System for Galveston Bay 

57
during September 1997, Proceedings of the Ocean Community Conference, MTS, Baltimore, 

MD, 426-430. 

Schmalz, R.A., 2000a: A Nowcast/Forecast System for Galveston Bay, Proceedings of the 2000 

Joint Conference on Water Resources Engineering and Water Resources Planning and 

Management, ASCE, General Conference, Flood Hazard Modeling Applications, CDROM. 

Schmalz, R.A., 2000b. Demonstration of a Nowcast/Forecast System for Galveston Bay, 

Proceedings of the 6th Estuarine and Coastal Modeling Conference, ASCE, New Orleans, LA, 

868-883. 

Schmalz, R.A., 2000c. High-Resolution Houston Ship Channel ADCP and CTD Survey: 

Model-Data Intercomparisons, NOAA Technical Report NOS CS 6, Silver Spring, MD. 

Schmalz, R.A., 2001. Three-Dimensional Hydrodynamic Model Developments for a Galveston 

Bay Nowcast/Forecast System, NOAA Technical Report NOS CS 9, Silver Spring, MD. 

Schmalz, R.A., 2003. Development of a Real Time Wave Model for Galveston Bay, Proceedings 

of the World Water & Environmental Resources Congress, ASCE, Philadephia, PA, CDROM. 

Schmalz, R.A., 2004. Development of the NOS Galveston Bay Operational Nowcast/Forecast 

System, Proceedings of the 7th International Marine Environmental Modelling Seminar, 

Washington, DC. 

 

Schmalz, R.A. and P. H. Richardson, 2002. Evaluation of the NOS Experimental 

Nowcast/Forecast System for a Galveston Bay, NOAA Technical Report NOS CS 14, Silver 

Spring, MD. 

Schwab, D. R., J.R. Bennett, P.C. Liu, and M.A. Donelan, 1984. Application of a Simple 

Numerical Wave Prediction Model to Lake Erie, Journal of Geophysical Research, Vol. 89, C3, 

3586-3592. 

Shen, J. and L. Haas, 2004. Calculating age and residence time in the tidal York River using 

three-dimensional model experiments, Estuarine, Coastal, and Shelf Science, 61, 449-461. 

Signell, R.P., R.C. Beardsley, H.C. Graber, and A. Capotondi, 1990. Effect of wave-current 

interaction on wind-driven circulation in narrow, shallow embayments. Journal of Geophysical 

Research, 95, 9671-9678. 

Smagorinsky, J. 1963. General circulation experiments with the primitive equations, I. The basic 

experiment. Monthly Weather Review, 91, 99-164. 

Smith, S.D. 1980. Wind stress and heat flux over the ocean in gale force winds, Journal of 

Physical Oceanography,10, 709-726. 

Wang, Chi-Fang, Ming-His Hsu, and Albert Y. Kuo, 2004. Residence time of the Danshuei 

River estuary, Taiwan, Estuarine, Coastal, and Shelf Science 60,381-393. 

58
Wanninkhof, R., 1992. Relationship Between Wind Speed and Gas Exchange Over the Ocean, 

Journal of Geophysical Research, Vol. 97, No. C5, 7373-7382. 

Wei, Eugene J., 1994. Lagrangian Drifter Trajectory Modeling. NOAA Technical Report NOS 

OES 004, NOS, NOAA, Silver Spring, MD. 

Wei, Eugene J., 2004. Solute Dispersion Modeling in New York Harbor. NOAA Technical 

Report NOS CS 18, NOS, NOAA, Silver Spring, MD. 

59