Development, Testing, and Application of 

a Coupled Hydrodynamic Surface-Water/

Groundwater Model (FTLOADDS) with Heat 

and Salinity Transport in the Ten Thousand 

Islands/Picayune Strand Restoration 

Project Area, Florida

By Eric D. Swain and Jeremy D. Decker

Prepared as part of the 

U.S. Geological Survey Priority Ecosystems Science Initiative

Scientific Investigations Report 2009-5146

U.S. Department of the Interior

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

KEN SALAZAR, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

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

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

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

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

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

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

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

U.S. Government.

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

reproduce any copyrighted materials contained within this report.

Suggested citation:

Swain, E.D., and Decker, J.D., 2009, Development, Testing, and Application of a Coupled Hydrodynamic SurfaceWater/Groundwater Model (FTLOADDS) with Heat and Salinity Transport in the Ten Thousand Islands/Picayune 

Strand Restoration Project Area, Florida: U.S. Geological Survey Scientific Investigations Report 2009-5146, 42 p.

ISBN 978 1 4113 2597 5
iii

Acknowledgments

Assistance in the development of model parameters and priorities was obtained through 

consultation with Kim Dryden and Terry Doyle at the U.S. Fish and Wildlife Service, Dewey 

Worth and Janet Starnes at the South Florida Water Management District, and David Bauman 

at the U.S. Army Corps of Engineers. 

Modeling pre- and post-processing support was provided by USGS colleague Melinda Lohmann, 

and biological studies support was provided by USGS colleagues Catherine Langtimm and 

Brad Stith.

v

Contents

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

1–Introduction ................................................................................................................................................2

1.1–Purpose and Scope.......................................................................................................................2

1.2–Description of Study Area and Water-Management History.................................................2

1.3–Proposed Restoration Efforts ......................................................................................................4

2–Model Code Enhancements.....................................................................................................................5

2.1–Integration of Surface-Water and Groundwater Flow and Transport Models....................5

2.2–Incorporation of Heat Transport into the Surface-Water Model...........................................7

2.3–Representing Surface-Water Flow over Hydraulic Barriers..................................................8

3–Model Input Development......................................................................................................................10

3.1–Surface-Water Model.................................................................................................................10

3.1.1–Topography .......................................................................................................................10

3.1.2–Frictional Resistance ......................................................................................................10

3.1.3–Representation of Weirs and Levees ...........................................................................13

3.1.4–Heat Transport Parameters and Atmospheric Data ..................................................13

3.1.5–Boundary Conditions.......................................................................................................15

3.2–Groundwater Model....................................................................................................................15

3.3–Restoration Changes...................................................................................................................19

3.3.1–Canal Plugs.......................................................................................................................19

3.3.2–Pump Stations and Spreader Channels.......................................................................19

3.3.3–Levees................................................................................................................................19

3.3.4–Road Removal and Culverts...........................................................................................19

4–Testing and Application of Coupled Model with Heat and Salinity Transport...............................20

4.1–Existing Conditions Calibration and Results............................................................................20

4.1.1–Calibration Parameters...................................................................................................20

4.1.2–Existing Conditions Simulation Results........................................................................22

4.2–Sensitivity Analysis......................................................................................................................30

4.3–Simulation of Picayune Strand Restoration Project Scenario.............................................33

4.4–Ecological Comparison: Manatee Locations ..........................................................................34

4.5–Model Limitations and Capability to Predict Restoration Changes.....................................39

5–Summary ...................................................................................................................................................40

6–References Cited .....................................................................................................................................41

Figures

1. Map showing Ten Thousand Islands model domain, Picayune Strand 

Restoration Project area, and water-management features ................................................3

2. Map showing major features of management or restoration Alternative 3D 

proposed for the Picayune Strand Restoration Project ........................................................4

3. Diagram showing surface-water flow over a weir.................................................................8
vi

4-9. Maps showing—

4. Land-surface elevation in the Ten Thousand Islands area ........................................11

5. Land-surface elevations in Ten Thousand Islands model area model with the

Picayune Strand Restoration Project represented........................................................12

6. Frictional resistance in the Ten Thousand Islands model area ..................................14

7. Location of Florida Automated Weather Network weather station, National 

Oceanic and Atmospheric Administration tidal stations, stations used to 

define groundwater boundaries, and stations for groundwater comparison 

in the Ten Thousand Islands model domain..................................................................16

8. Groundwater transmissivity values applied to the Ten Thousand Islands 

model area ..........................................................................................................................18

9. Initial groundwater salinity conditions in the Ten Thousand Islands 

model area ..........................................................................................................................21

10-13. Graphs showing—

10. Comparison of simulated existing conditions flows at the I-75 northern 

canal boundaries and measured flows at structure FU1............................................22

11. Comparison of measured and simulated flows at FU1................................................22

12. Comparison of measured and simulated stage at FU1................................................24

13. Comparison of measured and simulated water levels for selected sites in 

the Ten Thousand Islands area .......................................................................................25

14. Map showing surface-water salinity and temperature stations in the Ten 

Thousand Islands area ..............................................................................................................26

15-17. Graphs showing—

15. Comparison of measured and simulated salinity for selected sites in the 

Ten Thousand Islands area..............................................................................................27

16. Comparison of measured and simulated water temperature for selected 

sites in the Ten Thousand Islands area .........................................................................31

17. Comparison of measured and simulated water temperature for the Port of 

the Islands site ...................................................................................................................33

18-20. Maps showing—

18. Simulated surface-water temperature for December 21, 2003, for existing 

conditions, and Picayune Strand Restoration Project-implemented 

conditions in the Ten Thousand Islands model area ..................................................35

19. Manatee locations and measured surface-water salinity and temperature 

on November 16 and 22, 2002, in the Ten Thousand Islands model area.................36

20. Salinity of temperature nodes for manatee travel corridors in the model area.........37

Tables

1. Flow simulation statistics at structure FU1 for existing conditions and sensitivity 

runs ...............................................................................................................................................23

2. Transport simulation statistics at field sites for existing conditions and sensitivity runs ........29

3. Simulation statistics at field sites comparing Picayune Strand Restoration 

Project-implemented conditions with existing conditions ..................................................34

4. Simulation statistics for selected nodes in manatee travel corridors comparing 

Picayune Strand Restoration Project-implemented conditions with existing 

conditions.....................................................................................................................................38
vii

Conversion Factors, Acronyms, Abbreviations, and Datums

Multiply By To obtain

centimeter (cm) 0.3937 inch (in.)

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

square kilometer (km2) 0.3861 square mile (mi2)

cubic meter per second (m3/s) 35.31 cubic foot per second (ft3/s)

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

meter squared per second (m2/s) 10.76 foot squared per second (ft2/s) 

meter squared per day (m2/d) 10.76 foot squared per day (ft2/d) 

CERP Comprehensive Everglades Restoration Plan

DEM Digital elevation model

FAWN Florida Automated Weather Network

FTLOADDS Flow and Transport in a Linked Overland/Aquifer Density-Dependent System

GIS Geographic information system

LABINS LAnd Boundary INformation System

mg/L milligrams per liter

NOAA National Oceanic and Atmospheric Administration

PEV Percent explained variance

PSRP Picayune Strand Restoration Project

PSU Practical salinity units

RMSD Root mean square difference

RMSE Root mean square error 

SGGE Southern Golden Gate Estates

SFWMD South Florida Water Management District

SFWMM South Florida Water Management Model

s/m seconds per meter

SWIFT2D Surface-Water Integrated Flow and Transport in Two Dimensions

TIME Tides and Inflows in the Mangrove Everglades

TTI Ten Thousand Islands

USGS United States Geological Survey
viii

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Vertical coordinate information is referenced to the North American Vertical Datum of 1988 

(NAVD 88)

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83)

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

*Transmissivity: The standard unit for transmissivity is cubic foot per day per square foot times 

foot of aquifer thickness [(ft3/d)/ft2]ft. In this report, the mathematically reduced form, foot 

squared per day (ft2/d), is used for convenience.

Specific conductance is given in microsiemens per centimeter at 25 degrees Celsius (μS/cm 

at 25°C).
Abstract 

A numerical model application was developed for the coastal area inland of the Ten Thousand Islands (TTI) in southwestern 

Florida using the Flow and Transport in a Linked Overland/Aquifer Density-Dependent System (FTLOADDS) model. This 

model couples a two-dimensional dynamic surface-water model with a three-dimensional groundwater model, and has been 

applied to several locations in southern Florida. The model application solves equations for salt transport in groundwater and 

surface water, and also simulates surface-water temperature using a newly enhanced heat transport algorithm. One of the purposes of the TTI application is to simulate hydrologic factors that relate to habitat suitability for the West Indian Manatee. Both 

salinity and temperature have been shown to be important factors for manatee survival. The inland area of the TTI domain is the 

location of the Picayune Strand Restoration Project, which is designed to restore predevelopment hydrology through the filling 

and plugging of canals, construction of spreader channels, and the construction of levees and pump stations. The effects of these 

changes are simulated to determine their effects on manatee habitat.

The TTI application utilizes a large amount of input data for both surface-water and groundwater flow simulations. These 

data include topography, frictional resistance, atmospheric data including rainfall and air temperature, aquifer properties, and 

boundary conditions for tidal levels, inflows, groundwater heads, and salinities. Calibration was achieved by adjusting the 

parameters having the largest uncertainty: surface-water inflows, the surface-water transport dispersion coefficient, and evapotranspiration. A sensitivity analysis did not indicate that further parameter changes would yield an overall improvement in 

simulation results. The agreement between field data from GPS-tracked manatees and TTI application results demonstrates that 

the model can predict the salinity and temperature fluctuations which affect manatee behavior. Comparison of the existing conditions simulation with the simulation incorporating restoration changes indicated that the restoration would increase the period of 

inundation for most of the coastal wetlands. Generally, surface-water salinity was lowered by restoration changes in most of the 

wetlands areas, especially during the early dry season. However, the opposite pattern was observed in the primary canal habitat 

for manatees, namely, the Port of the Islands. Salinities at this location tended to be moderately elevated during the dry season, 

and unchanged during the wet season. Water temperatures were in close agreement between the existing conditions and restoration simulations, although minimum temperatures at the Port of the Islands were slightly higher in the restoration simulation as a 

result of the additional surface-water ponding and warming that occurs in adjacent wetlands.

The TTI application output was used to generate salinity and temperature time series for comparison to manatee field tracking data and an individually-based manatee-behavior model. Overlaying field data with salinity and temperature results from the 

TTI application reflects the effect of warm water availability and the periodic need for low-salinity drinking water on manatee 

movements. The manatee-behavior model uses the TTI application data at specific model nodes along the main manatee travel 

corridors to determine manatee migration patterns. The differences between the existing conditions and restoration scenarios can 

then be compared for manatee refugia. The TTI application can be used to test a variety of hydrologic conditions and their effect 

on important criteria.

Development, Testing, and Application of a Coupled 

Hydrodynamic Surface-Water/Groundwater Model 

(FTLOADDS) with Heat and Salinity Transport in the Ten 

Thousand Islands/Picayune Strand Restoration Project 

Area, Florida

By Eric D. Swain and Jeremy D. Decker
2 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

1–Introduction

The Picayune Strand Hydrologic Restoration Project (PSRP) is part of a much larger effort known to as the Comprehensive Everglades Restoration Plan (CERP), which provides a framework and guide to restore, protect, and preserve the water 

resources in central and southern Florida. CERP projects cover 16 counties and affect an area over 46,600 km2, with the overall 

goal of capturing unused freshwater flow and redirecting it to areas in need. The majority of the water will be used for environmental restoration, and any remaining water will be used to benefit agriculture and urban supply.

The purpose of the PSRP is to restore the natural hydrology of a 220 km2 rural area in western Collier County in southwestern 

Florida. This area was overdrained during the 1960s as part of a failed housing development (Southern Golden Gate Estates), 

and still retains the associated canals, structures and roadways. The proposed restoration process involves the removal or modification of these features with the goal of restoring the predevelopment hydrology and improving downstream coastal areas by 

reducing freshwater drainage and elevating groundwater levels. The effects of the restoration are difficult to predict, due to the 

complicating factors of groundwater/surface-water interactions, flat hydraulic gradients and topography, salinity exchange at the 

coast, and the complex coastal geometry. A numerical model that properly simulates these factors is needed to determine how 

PSRP water redistribution affects the frequency and duration of inundation (hydroperiod), salinity, and coastal discharge.

The Ten Thousand Islands (TTI) are a chain of islands and mangrove islets off the coast of southwestern Florida, just 

offshore from the PSRP area (fig. 1), and are an important coastal ecosystem. The Faka Union Canal is the major hydraulic connection between the PSRP and TTI areas (fig. 1). The Port of the Islands is a resort and marina on the Faka Union Canal and an 

important habitat for manatees and other aquatic life. Numerous field studies incorporating aerial surveys and satellite telemetry 

have been implemented to track manatee migration and determine the factors that affect their behavior (Langtimm and others, 

2006). Transmitters attached to the manatees allow them to act as mobile temperature and salinity probes (Stith and others, 2007). 

The Port of the Islands basin is used as a passive thermal refuge by manatees during periods of cold weather. Consequently, any 

changes to the hydrologic system must be considered in relation to their effects on these refugia if adverse changes to manatee 

habitat are to be avoided. 

An individually based manatee-behavior model is being developed that uses environmental data, including salinity and 

temperature, on a network structure using a Markov-Chain approach to predict manatee movement (Stith and others, 2006). For 

a hydrologic model to supply information to the manatee-behavior model, the ability to simulate heat and salinity transport in 

the PSRP and TTI areas is essential (Swain and Stith, 2006). In 2005, the U.S. Geological Survey (USGS) initiated a study to 

develop a coupled surface-water and groundwater model that incorporates salt transport and surface-water heat transport. The 

effort is conducted as part of the USGS Priority Ecosystems Science (PES) Initiative, with the purpose of understanding and 

predicting the results of the proposed PSRP changes. Model results will help quantify the effects to the ecosystem of the PSRP 

and supply input for the individual-based manatee-behavior model. Consequently, the model is an essential tool for providing 

reliable scientific information for the management of the water and biological resources of southern Florida.

1.1–Purpose and Scope

The purpose of this report is to describe the development, testing, and application of a flow and transport simulation model 

for the TTI area in support of ongoing Everglades restoration projects. The simulation code is summarized, including a description of the code structure and input requirements. The characteristics of the study area are presented, in addition to physical and 

hydrologic data used to develop a simulation of the period 1998–2005, and the calibration and sensitivity analysis. Measured 

and simulated values of two ecologically important parameters (salinity and temperature) are presented for current conditions 

and compared with corresponding values from simulations of proposed restoration changes.

1.2–Description of Study Area and Water-Management History

The TTI area is located within western Collier County and includes the TTI model domain (fig. 1). This application domain 

includes the Southern Golden Gate Estates (SGGE) area, South Belle Meade, Rookery Bay National Estuarine Research Reserve, 

Collier Seminole State Park, Ten Thousand Islands National Wildlife Refuge, and Fakahatchee Strand Preserve State Park. The 

SGGE area is of particular interest and encompasses about 243 km2, bordered by I–75 to the north and U.S. Highway 41 to the south. 

During the early 1960s, the Gulf American Corporation initiated work on the SGGE development, eventually dredging 

about 77 km of canals and constructing 467 km of shell-rock roads. Construction of the canal system began in 1968, and the 

main canal (Faka Union Canal) and the western feeder canal (Miller Canal) were completed by 1970. The two eastern feeder 

canals (Merritt Canal and Prairie Canal) were excavated in 1970–71, and the expansion of Prairie Canal was finished in 1973. 

Seven weir control structures were constructed to limit drainage from the area (fig. 1); however, annual surface-water runoff is 

typically more than double preconstruction amounts (Swayze and McPherson, 1977). 
1–Introduction 3

Figure 1. Ten Thousand Islands model domain, Picayune Strand Restoration Project area, and water-management features.

STUDY

AREA

FLORIDA

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

Naples

41

75

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

Barron River

Fakahatchee

Strand

Faka Union Canal

Port of the Islands

FU-1

South Belle Meade

Miller Canal

Merritt Canal

Prairie Canal

Okaloacoochee Slough

Southern Golden

Gate Estates

Faka Union Canal

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2855000

2865000

2875000

2885000

2895000

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

CONTROL STRUCTURE

SURFACE-WATER TIDAL

BOUNDARY

SURFACE-WATER NO-FLOW

BOUNDARY

SURFACE-WATER DISCHARGE

BOUNDARY TEN THOUSAND ISLANDS NATIONAL

WILDLIFE REFUGE

EVERGLADES NATIONAL PARK

FAKAHATCHEE STRAND PRESERVE

STATE PARK

PICAYUNE STRAND RESTORATION

PROJECT AREA

ROOKERY BAY NATIONAL ESTUARINE

RESEARCH RESERVE

COLLIER SEMINOLE STATE PARK

Note: Some shaded areas overlap

GULF OF MEXICO
4 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Canal construction efforts for the SGGE development predate State and Federal laws concerning drainage standards or the 

regulation of wetlands, and led to the over-drainage of adjacent lands, including the Fakahatchee Strand Preserve State Park 

(U.S. Army Corps of Engineers, 2004). Groundwater levels were lowered as much as 1.2 m near Prairie Canal and the low-water 

level was reduced by 0.6 m in the center of Fakahatchee Strand. Before channelization, groundwater gradients in the study area 

were oriented to the southwest toward Okaloacoochee Slough. These gradients shifted westward following canal construction, 

indicating that channelization affected groundwater movement in the area (Swayze and McPherson, 1977). Channelization also 

decreased inundation areas and periods within the study area (Klein, 1980). Examination of pollen and spore content in core 

samples from Fakahatchee Strand Preserve State Park indicate a gradual shift from a wet to relatively dry pollen and spore 

assemblage in the area during the past hundred years. This drying began after the construction of the first canals in the area in 

1928, and accelerated after logging burrows and tramways were built during the late 1940s and early 1950s. Completion of the 

SGGE canal system in 1973 further affected the drying process. This drying trend can be explained by comparing the timing of 

these construction efforts with observed changes in the gathered pollen samples (Donders and others, 2005).

1.3–Proposed Restoration Efforts

The changes implemented within the TTI application for the restoration scenario are intended to represent the recommended 

plan for the restoration of Picayune Strand as described by the U.S. Army Corps of Engineers (2004). This project is a part of 

the Comprehensive Everglades Restoration Plan (CERP) and proposes hydrologic changes to produce the following benefits: 

the restoration of historic wetland ecosystems, improved sheetflow and overland flow toward coastal areas, reduced surge flows 

through existing canals, improved aquifer recharge and reduced surface drainage, improved habitat for fish and wildlife, and a 

reduction in exotic species. Within the documentation, the final plan is referred to as Alternative 3D, shown in figure 2. This plan 

involves installing 83 canal plugs to the four canal reaches to restrict canal flow, three pump stations to supply spreader channels connected to Miller, Faka Union, and Merritt Canals, creating levees for lands adjacent to the Port of the Islands and two 

other residential areas, and removing roads in the SGGE area. The South Florida Water Management District (SFWMD) plans 

to maintain existing inflow quantities and patterns into Miller, Faka Union, and Merritt Canals at I–75 along the northern model 

boundary (fig. 1).

Figure 2. Major features of management or restoration Alternative 3D proposed for the Picayune Strand 

Restoration Project (PSRP). Modified from U.S. Army Corps of Engineers (2004).

SPREADER CHANNEL

PUMP STATION

CANAL PLUG

ROAD NOT REMOVED

AS PART OF PSRP

LEVEE (PROPOSED)

0 1 2 MILES

0 1 2 KILOMETERS

CONTROL STRUCTURE

81°40 ́ 81°30 ́

26°10 ́

26°00 ́

EXPLANATION

Private Lands

levee

POI SW levee POI N levee

41

75

Miller Canal

Faka Union Canal

Merritt Canal

Prairie Canal

35.4 m /s 3

74.5 m /s 3

22.7 m /s 3

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

FU1

Belle Meade
2–Model Code Enhancements 5

2–Model Code Enhancements

Several model code enhancements were necessary to meet the objectives of the study. These enhancements included integrating surface-water and groundwater flow and transport models, incorporating heat transport into the surface-water model, and 

representing surface-water flow over hydraulic barriers.

2.1–Integration of Surface-Water and Groundwater Flow and Transport Models

The numerical model applied to the TTI area is referred to as Flow and Transport in a Linked Overland/Aquifer DensityDependent System (FTLOADDS). Two preexisting models were combined to form FTLOADDS: the two-dimensional hydrodynamic surface-water model SWIFT2D (Schaffranek, 2004), and the three-dimensional groundwater model SEAWAT (Guo and 

Langevin, 2002). Both models can simulate variable-density salinity transport, and the coupling involves representing leakage 

and salt flux between surface water and groundwater. SWIFT2D computes vertically integrated two-dimensional forms of the 

equations of surface-water mass and momentum conservation, and transport equations for salt, heat, and other constituents. The 

equations for mass and momentum conservation in two dimensions are as follows (Leenderste, 1987):

∂

∂

+

∂( )

∂

+

∂( )

∂ = t

HU

x

HV

y

0

, (1)

∂

∂

+

∂

∂

+

∂

∂

− +

∂

∂

+

∂

∂

+

+ − U

t

U U

x

V U

y

fV g

x

gH

x

gU U V

C H

ζ aW

ρ

ρ θρ

2

2 2 1 2

2

2 ( ) si / nψ

ρH

k U

x

U

y − ∂∂

+

∂

∂











 =

2

2

2

2 0

, (2)

∂

∂

+

∂

∂

+

∂

∂

+ +

∂

∂

+

∂

∂

+

+ − V

t

U V

x

V V

y

fU g

x

gH

x

gV U V

C H

ζ aW

ρ

ρ θρ

2

2 2 1 2

2

2 ( ) co / sψ

ρH

k V

x

V

y − ∂∂

+

∂

∂











 =

2

2

2

2 0

, (3)

where: 

C is the Chezy resistance coefficient,

f is the Coriolis parameter,

g is acceleration due to gravity, 

H is temporal depth (h + ζ),

h is the distance from a horizontal reference plane to the channel bottom,

k is the horizontal exchange coefficient,

U is the vertically averaged velocity component in the x direction,

V is the vertically averaged velocity component in the y direction,

W is wind speed,

ζ is water-surface elevation relative to horizontal reference plane,

	 θ is the wind stress coefficient,

	 ρ is water density,

	 ρa is air density, and

ψ is the angle between wind direction and the positive y direction.

In these equations, x and y are Cartesian coordinates in the horizontal plane, and t represents time. Equation 1 expresses mass 

conservation, and equations 2 and 3 express momentum conservation in the x- and y-coordinate directions.
6 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

The equation for constituent transport and conservation in two-dimensions is as follows (Leenderste, 1987):

∂

∂

+

∂

∂

+

∂

∂ − ∂ ∂ ∂∂ − ∂ ∂ ∂∂

+ = ( ) HP ( ) ( ) ( / ) ( / )

t

HUP

x

HVP

y

HD P x

x

HD P y

y

HS x y 0

 , (4)

where Dx and Dy are the diffusion coefficients of dissolved substances, P is the vector of vertically averaged dissolved constituent 

concentrations, and S is the source of fluid with dissolved substances.

The code was modified for application to coastal wetlands, such as the Everglades, to allow the input of spatially variable 

rainfall, computation of evapotranspiration, variation of frictional resistance with depth, and other necessary features (Swain and 

others, 2004, Swain, 2005). SEAWAT combines the three-dimensional groundwater flow model MODFLOW with the solutetransport code MT3DMS to incorporate the effect of salt transport. The three-dimensional variable-density groundwater flow 

equation is as follows (Guo and Langevin, 2002):

0 0 0

0 ,0

0

' s

s s s

s

h C h z S q t C t       

 

  −  ∂ ∂ ∂ ∇⋅ ∇ + ∇ = + −       ∂ ∂ ∂  

K0

, (5)

where:

∇ is three-dimensional divergence; ∂/∂x + ∂/∂y + ∂/∂z ,

μ is dynamic viscosity,

μ0 is the dynamic viscosity of freshwater,

K0 is the hydraulic conductivity tensor, 

h0 is hydraulic head [L] measured in terms of equivalent freshwater,

Ss,0 is specific storage [L-1], 

	 θ is porosity,

Cs is salt concentration [ML-3], and 

q ́s is a source or sink [T-1] of fluid with density ρs. 

The solute transport equation is as follows (Guo and Langevin, 2002):

1+ 



 





∂( )

∂ = ∇⋅ ⋅ ( ) ∇ − ∇⋅( ) − +∑ ρ

θ

θ

θ b d

k

k k

s s

k

n

K C

t

D C Cq q C' R

, (6)

where: 

ρb is the bulk density (mass of the solids divided by the total volume), 

Kd is the distribution coefficient, 

Ck is the concentration of species k,

D is the hydrodynamic dispersion coefficient tensor, 

q is specific discharge [LT-1], 

Cs

k is the source or sink concentration of species k, and 

Rn is a reaction term.

The version of SEAWAT used in FTLOADDS does not incorporate heat transport and temperature simulations for groundwater, 

although a newer version of SEAWAT does have this capability (Thorne and others, 2006; Langevin and others, 2008) and could 

be incorporated at a later date. Combining SWIFT2D and SEAWAT to account for leakage and salt flux between surface water 

and groundwater is accomplished by programming FTLOADDS to call both models, calculate leakage and mass exchange, and 

pass the necessary information between them.

The TTI application uses the FTLOADDS formulation with modifications for wetland applications as described in Swain 

(2005) and Wang and others (2007). In addition to the features described in these studies, some additional code modifications were 

made for the TTI application, namely, incorporating heat transport for surface-water temperature simulations and representing flow 

over hydraulic barriers such as weirs and levees.
2–Model Code Enhancements 7

2.2–Incorporation of Heat Transport into the Surface-Water Model

To simulate ecologically important parameters in the TTI area, FTLOADDS must be able to compute temperature and heat 

transport in the surface-water system. The original version of SWIFT2D , called SIMSYS2D (Leenderste, 1987), had a limited 

ability to calculate radiation, conduction, and aerodynamic convection heat-flux terms; little associated literature; and lacked 

references and/or derivations for the equations and constants implemented. The reformulated version used in the current study 

allows the use of clearly defined parameters that can be derived from field measurements.

As part of the current study, modifications to the heat-transport formulation in SWIFT2D were tested in the Tides and 

Inflows in the Mangroves of the Everglades (TIME) simulation area (Wang and others, 2007) in Everglades National Park. In the 

test, latent heat flux was computed from precalculated potential evapotranspiration and was found to have insufficient accuracy 

(Swain and Decker, 2007). In contrast, a Dalton formulation that relies on wind speed, humidity, and temperature was used for 

computing both latent and sensible heat flux (Brutsaert 1982):

H

LCU e e

P LH

a w a = ρ β ( ) −

, (7)

H C SH = Ucp a  ( ) T T w a − , (8)

where 

HLH is latent heat flux, 

HSH is sensible heat flux, 

L is latent heat of vaporization,

C is the mass-transfer coefficient determined from field data,

U is wind velocity,

ew is saturation vapor pressure,

ea is ambient vapor pressure,

P is atmospheric pressure,

β is the ratio of the molecular weight of water to dry air,

cp is the constant pressure specific heat of air,

Tw is water temperature, and

Ta is air temperature.

As an alternative approach, equation 7 was replaced with the following to implement a Penman formulation for latent heat 

computations:

H

H H c r e e

r r LH

ASR R a p a s a

s a

= ( ) − + −

+ + ( )

∆

∆

( / )( )

/

ρ

γ 1 , (9)

where 

HASR is absorbed solar radiation,

HR is long-wave radiation exchange flux,

∆ is the slope of the saturation vapor-pressure curve = (ew-es)/(Tw-Ta),

γ is the psychometric constant = cpP/βL,

rs is bulk stomatal resistance, and

ra is the aerodynamic resistance term. 

Equation 8 is still used for computing sensible heat when the Penman method is used.

In this study, the Dalton formulation (incorporating equations 7 and 8) was compared to the Penman formulation (incorporating 

equations 8 and 9). The Dalton and Penman formulations produced similar daily-average temperatures with slightly different 

daily fluctuations (Decker and Swain, 2008). The Penman equation is considered more physically based and, for that reason, is 

used in the TTI application.
8 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Radiation and sensible heat-flux terms were computed using solar radiation, net radiation, wind, humidity, and temperature 

data collected at field sites. The field data in the TIME domain were used in the TTI application to compute albedo, the sensibleheat mass-transfer coefficient, and water and air emissivity coefficients. To compute latent heat using the Penman method, the 

aerodynamic resistance term was computed using the formulation in Wang and others (2007, p. 6) and the bulk stomatal resistance was calibrated in the TIME simulation. The model formulation was also modified to account for the additional heat capacity of the soil underlying the wetland flow. The simulated surface-water temperatures in the TIME domain closely matches daily 

average temperatures measured inland and at the coast, indicating the utility of the formulation and coefficient development 

(Decker and Swain, 2008).

2.3–Representing Surface-Water Flow over Hydraulic Barriers

SWIFT2D originally required four separate subroutines to represent flow over hydraulic structures. The formulation defined 

boundary conditions on model grid rows or columns that are defined by the structure flow equations (Schaffranek, 2004). 

Locations at which hydraulic structures are defined must be specified as boundaries, and the flows computed by the structure 

equations used as boundary conditions. For the current study, it was decided to replace this approach with a simpler one that 

computes an effective Chézy coefficient for the model-cell face that represents the same resistance as the hydraulic structure. 

Because the structure flow equation and the Chézy formula do not have the same proportionality between discharge and water 

levels, the effective Chézy coefficient is a function of the water level and hydraulic gradient. The weir equation for the situation 

shown in figure 3 can be stated as follows:

Q C h h h h = − w u ( ) p w up − 1 , (10)

where 

Q is discharge, 

Cw is the weir coefficient (which can vary depending on flow conditions), 

hup is the headwater stage, 

hw is the elevation of the weir-crest, and 

h1 is an elevation that depends on submergence as: h1 = hw if hdwn ≤ hw (free flow condition) where hdwn is the 

tailwater stage. h1 = hdwn if hdwn > hw (submerged flow condition).

Figure 3. Surface-water flow over a weir.

hup hdwn

hw

hup

hdwn

hw

DOWNSTREAM STAGE

UPSTREAM STAGE

ELEVATION OF WEIR CREST

EXPLANATION

Water surface

Land surface

Elevation datum
2–Model Code Enhancements 9

The Chézy equation for a cell face can be stated by the equation:

Q C w h z h h

l z up

up dwn = − ( )1 5. −

, (11)

where

Cz is the Chézy friction coefficient, 

z is the land elevation, 

l is the length between headwater and tailwater stages, and 

w is the flow width.

The SWIFT2D model grid consists of square cells, meaning l = w; therefore, equations 10 and 11 can be combined to yield 

the following equation:

C C

w

h h

h z

h h

h h z

w up w

up

up

up dwn

= ( ) −

( ) −

−

− 1 5

1

.

. (12)

Substituting in the definitions of h1, the following are obtained:

C

C

w h h

h h

h z z

wf

up dwn

up w

up

=

( ) −

−

−













1 5.

 if hdwn ≤ hw, and (13)

C C

w

h h

h z z

ws up w

up

= ( ) −

( ) − 1 5.

 if hdwn > hw, (14)

where Cwf and Cws are values of Cw that differ for the two flow conditions, namely, free flow and submerged flow. These coefficients 

have the dimensions of length1.5/time, and are related to two dimensionless flow coefficients Cdf (for free flow) Cds (for submerged 

flow) by the following two equations:

C C w g wf = df 2 , (15)

C C w g ws = ds 2 , (16)

where g is gravitational acceleration. Cdf, Cds, w, and g are constants, and therefore, have been combined in Cwf and Cws. Because 

flow can occur in either direction over the weir, hup and hdwn are the higher and lower stages, respectively, regardless of location.

Equations 13 and 14 are used in SWIFT2D to set the Chézy friction coefficient for the cell faces where weirs are defined. 

This representation can be applied at other types of surface-water barriers such as levees and roads; in this case, the flow coefficients must be chosen to represent the effective resistance of the barrier. For levees and roads that are rarely overtopped, the 

accuracy of the coefficients is less crucial.
10 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

3–Model Input Development

A large variety of data input is needed for FTLOADDS, given the information needed for the formulation and the inclusion of 

both surface-water and groundwater regimes. The 8-year simulation period used in this study is from 1998 through 2005. Model 

data input includes topography and bathymetry, frictional resistance, hydraulic structure coefficients, parameters controlling 

heat fluxes between the atmosphere and land (including evapotranspiration), wind speeds and directions, boundary levels, flows, 

salinities, aquifer depths, and hydraulic conductivities. Less data are available in the TTI domain than in the TIME domain in 

Everglades National Park (Wang and others, 2007), which has been studied extensively.

3.1–Surface-Water Model

The surface-water simulation was executed with a 2-minute timestep and a uniform finite-difference grid with 500-m spacing. 

The maximum model dimensions are 94 cells east-to-west and 85 cells north-to-south. Two constituents were simulated, salt and 

heat, and the parameter that defines the minimum depth at which a cell is considered wet was set to 0.1 m above land surface. 

Stage was assigned a spatially uniform value of 0.01 m NAVD 88.

3.1.1–Topography

Topographic data for the TTI application included land-surface elevations and oceanic bathymetry. Data used to assign land 

surface elevations were compiled from the State of Florida Land Boundary Information System (LABINS) (Florida Department 

of Environmental Protection, 2008). Digital elevation models (DEMs) provided from LABINS have a resolution of 30 m and 

were downloaded based on the Florida Department of Environmental Protection quadrangles. These models were then resampled using ArcMAP software (ESRI, 2008) for 50-m grid cells. Combined land elevation and bathymetry data were obtained 

from the National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center. As with the DEM 

data, the bathymetry data were resampled at a grid-cell spacing of 50 m. The DEM data from LABINS were used for inland cell 

elevations with coastal and offshore elevations specified from NOAA bathymetry data. The decision to use DEM or bathymetry 

values for cells along the coast was made based on elevation values, average tidal values, and geographic location. The 50-m 

DEM/bathymetry values were then averaged into larger 500-m cells and used for the model simulations.

Major canals in the model area range from 12 to 60 m in width, which is substantially smaller than the model grid size. It 

was decided to represent the canals by changing the properties of the corresponding grid cells―a method commonly used to 

represent smaller scale features in a numerical model without reducing the grid size (Neelz and Pender, 2007; Wang and others, 

2007). To provide a model cell volume equivalent to the canal and the land within the cell, the canal bottom elevation and the 

adjacent land elevation are weighted by the ratio of the canal width and the cell width. Figure 4 shows the land-surface elevations in each model grid cell for the TTI domain.

The topographic grid for the PSRP scenario represents the plugging of the canals shown in figure 2. Multiple canal plugs 

within a single model cell are represented as the removal of the canal, because the plugs are designed to fill the canal to land 

surface. This representation is necessary to remove the storage of the canal reach. Figure 5 shows the land-surface elevations 

in each model grid cell for the TTI application PSRP scenario. The most noticeable change is to the areas where short spreader 

channels have been added (fig. 2).

3.1.2–Frictional Resistance

The process of determining the frictional resistance of model cells began with assigning onshore/offshore/coastal cell types. 

Cells determined to be offshore based on cell topography and mean tidal water levels were assigned Manning’s n values based 

on a standard open-water value of 0.02 s/m1/3 (French, 1985, p. 130). Coastal cells were those only partially inundated at mean 

tidal levels, and were assigned weighted values of Manning’s n proportional to the fraction of each cell offshore and in the 

coastal wetlands. The values used for this weighting were again 0.02 s/m1/3 for offshore portions and a nominal wetland value of 

0.4 s/m1/3 (Schaffranek and others, 1999). 

The overland frictional resistance for onshore model cells was derived from several sources. First, existing land-use 

maps and aerial photography were used to determine vegetation types present within model cells. The primary source for 

this information was the Existing Condition Land Use Map (U.S. Army Corps of Engineers, 2004), which was derived from 

SFWMD 2000 land-use Geographic Information System (GIS) coverage. Within this map, the hundreds of land-use categories 

originally used by SFWMD were reduced into 21 different vegetation types. Further simplifications were made in assigning 

model cells a vegetation type based on this coverage to obtain seven primary vegetation/land-use types within the model area. 
3–Model Input Development 11

Figure 4. Land-surface elevation in the Ten Thousand Islands area. Control structure FU1 location is shown for 

reference.

FU1

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

ELEVATION, IN METERS NAVD 88

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

-12 -10 -8 -6 -4 -2 0 2 4

GULF OF MEXICO

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

CONTROL STRUCTURE
12 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 5. Land-surface elevations in Ten Thousand Islands model area model with the Picayune Strand Restoration 

Project represented. Surface-water barriers are shown in white. Control structure FU1 location is shown for reference.

-12 -10 -8 -6 -4 -2 0 2 4

GULF OF MEXICO

ELEVATION, IN METERS NAVD 88

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

FU1

EXPLANATION

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

0 5 10 MILES

0 5 10 KILOMETERS

CONTROL STRUCTURE
3–Model Input Development 13

Manning’s n values, which ranged from 0.25 to 0.45 s/m1/3, were assigned to each vegetation type based on research conducted 

within Everglades National Park (Lee and Carter, 1999), the Project Implementation Report (U.S. Army Corps of Engineers, 

2004), and standard published values (French, 1985).

The frictional resistance of the canal banks is also represented by a cell-face Manning’s n value. Because the flow depths 

along the canal bank tend to be very low, especially when a berm exists, the effective frictional resistance at the overland/canal 

interface tends to be higher than in other cells. The South Florida Water Management Model (SFWMM) represents regional 

overland and canal flow on a larger scale grid. The interaction of the overland flow and the canals is represented by a resistance 

term of the form (South Florida Water Management District, 2005):

n adb = , (17)

where n is Manning’s friction, d is the water depth, and a and b are constants. For the connection of overland flow and canals, 

the values in the SFWMM varied from a = 0.5 to 2.0, and b = -0.77 to 0 (South Florida Water Management District, 2005). An n

value of 2.0 s/m1/3 was used for the canal banks in the TTI application and seems consistent with those used in the SFWMM.

In addition to using vegetation maps, small-scale topographical features were represented by manipulating the frictional 

resistance across overland model cell faces. Major canals and waterways were mapped within model cells, and the Manning’s n

value of the cell faces are calculated as the average of the canal and adjacent wetland resistances weighted by the ratio of the 

canal width to the remainder of the cell face. The preceding methods were used to obtain the Manning’s n frictional resistance 

map shown in figure 6. The color defining the frictional resistance of a given cell face is centered on the model cell face as 

shown in the figure inset. The anisotropic effects of incorporating the canals are evidenced by the low resistance along the canals 

and the high resistance along the canal banks (fig. 6). As part of the sensitivity analysis presented later, the Manning’s n values 

are adjusted to examine their effect on model results.

3.1.3–Representation of Weirs and Levees

The stage and flow data collected at FU1 were used to estimate values of Cwf and Cws for the four weirs in the model area 

(fig. 1). Using the 1998–2005 period of record for water levels and flows at FU1 and equation 10, mean values of Cwf = 

104.3 m1.5/s and Cws = 202.6 m1.5/s were calculated. At FU1, w = 57 m, which yields corresponding dimensionless flow coefficients, Cdf = 0.413 and Cds = 0.802. These values are similar to values developed by SFWMD for weirs of this design: 0.36 

and 0.90 for Cdf and Cds, respectively (Otero, 1994). The calculated values for Cwf and Cws are used in the model for the weir at 

FU1. The coefficients Cwf and Cws for the other weirs in the model area were calculated by multiplying the computed Cdf and Cds

values by w g2 at each structure, yielding the following values:

Weir name Width w

(m)

Cwf, 

using Cdf = 0.413

(m1.5/s)

Cws, 

using Cds = 0.802

(m1.5/s)

FU1 57 104.3 202.6

Miller-1 15 27.45 53.32

FU2 22 40.26 78.20

Merrit-1 6.5 11.89 23.10

For the PSRP scenario shown in figure 2, hydraulic barriers are used to represent the proposed levees, as well as canal plugs 

separated by a distance greater than one cell width. The canal plugs are specified to be at land surface, and the elevations of the 

levees are specified from data in U.S. Army Corps of Engineers (2004). The flow coefficient for levee overtopping is estimated 

by using equation 10 for the case when submerged with a Chezy value equivalent to the average Manning’s n of 0.4 s/m1/3 and 

using a nominal depth of water of 0.6 m; yielding a Cws value of 22.4 m1.5/s.

3.1.4–Heat Transport Parameters and Atmospheric Data

The Dalton formula was applied initially using equations 7 and 8, which use the same mass-transfer coefficient C for sensible heat and latent heat. Mean simulated temperature is inversely proportional to the mass-transfer coefficient, and calibrating 

the coefficient to match mean temperatures is straightforward. When the Penman equation is used to compute latent heat using 

equation 9, the Dalton equation 7 is still used for sensible heat computations and incorporates the same mass transfer coefficient 

as when both equations 7 and 8 are used.
14 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 6. Frictional resistance in the Ten Thousand Islands model area. Control structure FU1 location is shown 

for reference.

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

FU1

GULF OF MEXICO

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Manning’s n

0.1 0.2 0.3 0.4 0.5

Orientation of model grid and colors

indicating cell-face friction values
3–Model Input Development 15

Based on parameters developed for the TIME model area (Decker and Swain, 2008), albedo is set to 0.169, water emissivity to 0.95, and for the Penman formula, the aerodynamic resistance is set to 160 divided by wind speed (in meters per second) 

and the bulk stomatal resistance to 150 s/m. Time varying input includes wind-speed, solar radiation, relative humidity, and air 

temperature; all of these data were taken from the Florida Automated Weather Network (FAWN) site in Immokalee, 38.6 km 

north of the TTI domain (fig. 7). The temperature of the surface-water inflows to the model area were estimated due to a lack of 

data. Water temperature can be represented by lagged air temperature (Shoemaker and others 2005), and was estimated as the 

average air temperature for the preceding 24 hours. This average was also used for the temperature of groundwater to surfacewater leakage, even though groundwater temperature should not be as closely related to air temperature as is surface-water 

temperature. As part of the sensitivity analysis, the mass-transfer coefficient and the bulk stomatal resistance values were varied.

Rainfall data collected at the FAWN site at 15-minute intervals were applied to the TTI area averaged over 6-hour intervals. 

Insufficient data exist to define a spatial distribution of rainfall, and values were applied uniformly to the model domain. The 

temperature of incoming rain Train was estimated by the Magnus-Tetens formula (Barenbrug, 1974):

T C rain

o

= −





237 7

17 27

.

. , (18)

where

 = ++

17 27

237 7

.

. ln( ) TC T

RH air

o

air , (19)

where Tair is air temperature in degrees Celsius and RH is relative humidity.

3.1.5–Boundary Conditions

Oceanic water-level boundaries were used for the western and southern boundaries of the model domain. Data for these 

boundaries were obtained from the NOAA tidal stations shown in figure 7 and were obtained through the NOAA Tides and Currents portal, which provides access to its historical oceanographic and meteorological data (National Oceanic and Atmospheric 

Administration, 2008). Tidal water-level data were collected from January 1, 1998, to December 31, 2005. The tidal values were 

varied as part of the sensitivity analysis.

The inland surface-water boundaries are specified as no-flow conditions with the exception of the inflows to Miller, Merritt 

and Faka Union Canals at I–75 along the northern model domain boundary. Because discharge is continuously measured only at 

weir FU1 (fig. 1), the sum of the flows along the northern canal boundaries was initially assumed to be equal to flow at FU1. The 

land surface in the cells between the northern boundaries of the Merritt, Faka Union, and Miller Canals was lowered, and the 

discharge rate at FU1 was distributed across the lowered cells. Adjusting this flow to achieve better simulated flows at FU1 was 

part of the calibration process and was tested as part of the sensitivity analysis discussed later.

3.2–Groundwater Model

The shallow surficial aquifer system in the TTI domain is composed of shelly limestone and shelly, marly sand (Klein, 1980). 

The top of the surficial aquifer system is coincident with land surface, which varies in elevation from zero to about 5 ft above 

NAVD 88 (fig. 4). The bottom of the surficial aquifer system meets a confining layer above the Floridan aquifer system at an 

elevation of about -15 m NAVD 88 (U.S. Army Corps of Engineers, 2004). Therefore, the thickness of the surficial aquifer 

system is about 15 to 20 m onshore, and decreases to about 4 m offshore.

The horizontal configuration of the groundwater model grid matches the surface-water model uniform finite-difference grid 

with 500-m spacing. A daily stress-period length is specified as the time between successive solutions of equation 5. The TTI 

application represents groundwater flow down to the confining layer and leakage to surface-water features (wetlands, canals, 

streams, and the Gulf of Mexico). The bottom of the lowest model layer is a no-flow boundary, as are all of the lateral offshore 

boundaries. The inland groundwater boundaries are defined head boundaries using values spatially interpolated from the stations 

shown in figure 7. Rainfall recharge is accounted for in the surface-water model and can be transferred to groundwater through 

leakage. Evapotranspiration is also computed in the surface-water model and is removed directly from the water table when 

surface water is absent.
16 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 7. Location of Florida Automated Weather Network (FAWN) weather station, National Oceanic and 

Atmospheric Administration (NOAA) tidal stations, stations used to define groundwater boundaries, and stations 

for groundwater comparison in the Ten Thousand Islands domain. Control structure FU1 location is shown for 

reference.

FAWN

weather

station

Study

area

NOAA 8725110

GOLD.W1

GOLD.W2

GOLD151

HEND84

MLRI75 FAKI75

MRTI75

SGT2W2

SGT2W6

SGT4W1

BARRON

EXPLANATION

UTM NORTHING, IN METERS

UTM EASTING, IN METERS

41

75

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

Barron River

Faka Union Canal

FU-1

Miller Canal

Merritt Canal

Prairie Canal

Faka Union Canal

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2855000

2865000

2875000

2885000

2895000

GULF OF MEXICO

0 5 10 MILES

0 5 10 KILOMETERS

NOAA TIDAL STATION

STATION FOR GROUNDWATER BOUNDARY HEAD

GROUNDWATER HEAD COMPARISON STATION

GROUNDWATER BOUNDARY HEAD INTERPOLATED BETWEEN

GROUNDWATER STATIONS

GROUNDWATER BOUNDARY HEAD INTERPOLATED BETWEEN

GROUNDWATER STATIONS AND TIDAL STATION

GROUNDWATER NO-FLOW BOUNDARY

CONTROL STRUCTURE
3–Model Input Development 17

To delineate vertical salinity variation in the aquifer, the groundwater model was divided into three layers of equal thickness. The thickness of each layer varies from 1.3 to 6.7 m, according to the total depth of the aquifer. Horizontal conductivity, 

vertical conductivity, and transmissivity were assigned vertically uniform values, because insufficient information exists to 

define vertical variation in these properties. The specific yield of the upper layer was assigned a value of 0.1 and the confined 

storativity of the lower layers was assigned a value of 0.00001. The cells in the uppermost layer can range between confined and 

unconfined conditions depending on the elevation of the simulated groundwater level relative to the elevation at the top of the 

cell. For unconfined conditions in which the groundwater level is lower than the top of the model cell, the saturated thickness 

of individual cells in the upper layer (head minus cell bottom elevation) is used to calculate transmissivity and storativity; for 

confined conditions in which the groundwater level is higher than or equal to the top of the model cell, the cell thickness is used. 

The lower two layers are always treated as confined, and their transmissivity and storativity values are specified directly in the 

relevant model input files.

Groundwater conductivity and transmissivity values were assigned for the TTI domain based on information compiled from 

Klein (1980) and the U.S. Army Corps of Engineers (2004). The horizontal hydraulic conductivity map by the U.S. Army Corps 

of Engineers (2004) shows a conductivity range within the TTI domain of 0.0035 to 0.132 ft/s (92–3,476 m/d). Previous investigations lack information for the TTI area, but provide transmissivity ranges for nearby areas (Klein, 1980) that can be used to 

estimate TTI values and confirm those used by the U.S, Army Corps of Engineers (2004). The ranges provided in these investigations yield an overall transmissivity range of 8,000 to 134,000 ft2/d (2,438- 40,843 m2/d). This overall range corresponds to a 

horizontal hydraulic conductivity range of about 120 to 2,700 m/d, given the previously mentioned 15- to 20-m thickness range 

for the aquifer in the onshore region. Based on these data, the conductivity within the individual layers of the model was set to a 

constant value of 1,500 m/d, which yields a transmissivity range of about 2,000 to 10,000 m2/d for the individual model layers. 

A map of the resulting aquifer transmissivity values is shown in figure 8. Vertical conductivity in the aquifer was set to 150 m/d, 

reflecting a vertical to horizontal anisotropy ratio of 1:10.

Leakage between the groundwater and surface water is computed in FTLOADDS by the following equation:

q

K

z h h z leak

mean

sw

f

f

= ∆− +

 −









 ∆













 1  1

1

1

 

 , (20)

where 

qleak is leakage per unit area, 

∆z1 is the distance that leakage occurs over, 

hsw is the surface-water equivalent freshwater head, 

h1 is the groundwater top-layer equivalent freshwater head, 

ρ1 is the density of the groundwater in the top layer, and 

	 ρf is freshwater density, and 

Kmean is the mean hydraulic conductivity that leakage occurs over, given by:

K z

z

K

z z

K

mean

tl

tl

tl

= ∆

∆

+

∆ − ∆

1

1

1 , (21)

where ∆ztl is the thickness of a thin surficial layer, Ktl is the hydraulic conductivity of the thin surficial layer, and K1 is the 

hydraulic conductivity of the upper aquifer layer. 

Insufficient information exists to define a spatial distribution for the parameters in equations 20 and 21, and therefore, 

spatially uniform values were assigned for ∆ztl, Ktl, and K1. The formulation in FTLOADDS uses values for conductivity and 

thickness for a leakage layer that is assumed over the top layer. Wetlands typically have a less permeable organic litter layer 

that impedes leakage (Harvey and others, 2000). Due to numerous gaps in this litter layer, a vertical hydraulic conductivity of 

350 m/d and a thickness of 0.5 m were used for the TTI application. The hydraulic conductivity value is higher than what is 

typically used for organic peats, and was borrowed from the TIME application for Everglades National Park (Wang and others, 

2007). The sensitivity of model results to the leakage value was tested and is presented later. 

Groundwater head was assigned an initial elevation of 1.0 m NAVD 88 across the model domain, and represents the 

approximate average head value between the southern part of the model domain (sea level) and the northern model boundary 

(over 2 m NAVD 88). Initial conditions dissipate within the first few months as a result of boundary and leakage flows. Although 

salinity concentration is initially assigned a value of zero across the model domain, the final calibration uses a multiple iterative 

run method (described later) to calculate initial salinities.
18 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 8. Groundwater transmissivity values applied to the Ten Thousand Islands model area. Control structure 

FU1 location is shown for reference.

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

TRANSMISSIVITY, IN METERS SQUARED PER DAY

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

EXPLANATION

1,800 2,850 3,900 4,950 6,000 7,050 8,100 9,150

FU1

GULF OF MEXICO

0 5 10 MILES

Base from U.S. Geological Survey digital data, 0 5 10 KILOMETERS Universal Transverse Mercator Zone 17

CONTROL STRUCTURE
3–Model Input Development 19

3.3–Restoration Changes

To represent the PSRP scenario in the model, the canal topography was modified as shown in figure 5, and barriers representing canal plugs and levees were added. All other parameters, including surface-water inflows, remained unchanged. It was 

assumed that the PSRP-implemented changes do not appreciably affect hydrology external to the model domain.

3.3.1–Canal Plugs

Alternative 3D specifies the placement of 83 canal plugs within the Merritt, Miller, Faka Union, and Prairie Canals (fig. 2). 

As noted earlier, cell topography for the region was altered to represent these canals. The plugs were positioned based on the 

plan and location map included in U.S. Army Corps of Engineers (2004, app. C). In cases where multiple plugs were positioned 

within one model cell or several adjacent model cells, the cell topography was set to land-surface elevation. Where plugs are 

more widely spaced, they are represented by hydraulic barriers (broad-crested weirs). Barriers were assigned elevations equal to 

those of adjacent land cells. 

3.3.2–Pump Stations and Spreader Channels

The proposed restoration incorporates pumping stations on Miller, Merritt, and Faka Union Canals (fig. 2) to provide 

drainage and mitigate adverse backwater effects. Pump flow rates were determined by estimating the maximum output required 

to maintain current headwater levels, even during a flood event with a 100-year return period. Maximum computed pump flow 

rates were 35.4 m3/s for Miller Canal, 74.5 m3/s for Faka Union Canal, and 22.7 m3/s for Merritt Canal; only a fraction of each 

of these flows is required during normal operation. Downstream from these pump stations, spreader channels are proposed with 

the following lengths: 1,372 m (Miller Canal), 2,134 m (Faka Union Canal), and 427 m (Merritt Canal). Pump stations are to 

be operated so that existing headwater levels remain unchanged from current conditions; therefore, existing condition head and 

flow boundaries were used at the top of the model within each of the three major canals to simulate the restoration scenario. The 

flows travel down the remaining canals into spreader channels. These channels are represented by reducing the topography of 

the corresponding model cells. The land-surface elevation of these model cells is set to the area-weighted average of the channel 

bed elevation and the surrounding land surface. The locations for these channels were determined using the plan and location 

map in U.S. Army Corps of Engineers (2004, app. C). 

3.3.3–Levees

Alternative 3D includes the addition of three sets of proposed levees to protect private properties and developments (fig. 2), 

as described in U.S. Army Corps of Engineers (2004). The first set surrounds the Port of the Islands and two adjacent residential areas, and consists of three levees totaling about 6.7 km in length. The second set, referred to as the 6L levees, protects the 

agricultural areas in southern Belle Meade from the expected higher water levels and is about 19.6 km in length. The final set, 

referred to as the Private Lands levees, forms a ring around private properties in northeastern Belle Meade and is about 6.8 km 

in length. The proposed levee heights vary from 1.8 to 2.7 m depending on location, and are typically 4.6 m wide. For the restoration scenario, these levees were represented in the model application by adding cell-face barriers with heights equal to the 

proposed levee heights to impede flow in the same direction.

3.3.4–Road Removal and Culverts

The proposed restoration specifies the removal of about 365 of the 449 km of existing roads within the SGGE area (fig. 2; 

U.S. Army Corps of Engineers, 2004). The goal of this removal is to lower the affected topography to a level that no longer 

impedes surface-water sheetflow through the area. In addition, the construction of nine culverts is proposed to reduce flow 

resistance toward the coast. Some of these culverts would supplement existing bridges and culverts along U.S. 41, and others 

would help distribute water in the upper portions of the SGGE area and under remaining roads in the area. Given the coarseness of the model grid (500 m × 500 m), the effects of smaller scale modifications such as the road removal and addition of 

culverts just described cannot be represented explicitly. Therefore, these changes were not incorporated into the restoration 

simulation. 
20 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

4–Testing and Application of Coupled Model with Heat and Salinity Transport

The TTI application was run in two separate configurations for the 1998–2005 simulation period. The existing conditions 

simulation uses the field-derived data described in sections 3.1 and 3.2, and the restoration simulation uses the PSRP modifications described in section 3.3. Both simulations assumed identical discharge conditions where the canals connect to I–75 at the 

northern domain boundary, as well as identical tidal conditions at the coast. Climatic conditions and the spatially distributed 

parameters are also identical between simulations.

4.1–Existing Conditions Calibration and Results

Iterative least-squares best-fit procedures are used in many parameter-estimation methods (Poeter and others, 2005), especially to calibrate groundwater models. These methods are especially useful for cases involving a large number of parameters 

with a high degree of associated uncertainty. The need for the parameter-estimation method decreases as the amount of field 

information available for a given model area increases. When few parameters are uncertain, a calibration can be based on simple 

adjustment of one parameter at a time. The calibration of the TTI application was greatly simplified by the large amount of field 

data already available, as described earlier. This left relatively few variables with high uncertainty, primarily: (1) the defined 

flows at the I–75 northern canal boundaries, (2) the constituent dispersion coefficient, (3) the evapotranspiration parameters, and 

(4) groundwater conductivity. Although the effects of these variables are interrelated, the first three are surface-water parameters 

that primarily affect the flow, salinity transport, and heat transport computations, respectively. Therefore, inflows were calibrated 

first, followed by the dispersion coefficient and the evapotranspiration parameters, and groundwater hydraulic conductivity.

The extent of salinity intrusion in the shallow aquifer affects inland surface-water salinity because substantial water and 

salt exchange occurs between the aquifer and overlying canals, streams, and wetlands. Because the saltwater interface can take 

decades to reach equilibrium, multiple simulations of the same 8-year period were used to approximate a longer period. The 

8-year simulation was repeated 12 times, with the final simulated aquifer salinity distribution for each simulation used as the 

initial salinity distribution for the following simulation. Salinity is expressed in practical salinity units (PSU), which is the ratio 

of the electrical conductivity to the conductivity of a 1 mg/L salt solution. The resulting salinity distribution after these multiple 

runs was then used as the initial salinity distribution in the existing conditions and PSRP simulations (fig. 9) so that any longterm change in aquifer salinity caused by restoration changes was not included.

4.1.1–Calibration Parameters

The initial run using the northern-boundary surface-water inflows described earlier overestimated the flow at FU1 because 

of additional overland flow, and possibly, groundwater leakage into the canal. Assuming that the frictional resistance of the channel banks is reasonable (as described in section 3.1.2), the error was corrected by altering the inflow at I–75 by the discrepancy 

in flows at FU1. The revised flows at I–75 shown in figure 10 were used in the calibrated simulations. The effects of varying the 

inflows are discussed later as part of the sensitivity analysis.

To obtain an accurate surface-water dispersion coefficient, the most easily measured and sensitive quantity was used for 

matching, namely, the extent of seawater intrusion along the coast. Salinity transport along the Faka Union Canal to the Port 

of the Islands was more sensitive to the dispersion coefficient than other areas along the coast. A dispersion coefficient value 

of 100 m2/s yielded salinity values that nearly coincided with measured values at this location. Reasonable salinities were also 

simulated at other coastal points using this value, and insufficient data were available to justify defining spatial variation in the 

dispersion coefficient.

Parameters for the heat transport external flux terms were developed for the TIME model area (Decker and Swain, 2008) 

and applied to the TTI application. These parameters (and assigned values) include the following: water-surface emissivity 

εw = 0.95; coefficients to compute atmospheric emissivity εa = 0.80 + 0.019 ea

0.5, where ea is atmospheric vapor pressure; in 

equation 8, the mass transfer coefficient C = 0.0020; albedo α = 0.169; and in the Penman formulation equation 9, the stomatal 

resistance term rs = 150. The average simulated temperatures in the TTI area tended to be too low with these parameters. The 

mass-transfer coefficient was considered to be the most viable parameter to adjust, because accepted values tended to be lower 

than those taken from the TIME application, which range from 0.0009 to 0.0018 (Brutsaert, 1982, p. 126). The accepted values 

are established for open-water sites, and the predominance of the canal systems in TTI-area hydrology makes it more likely that 

these may be more suitable than the values from the TIME application, especially because temperatures are mostly compared at 

open-water locations. A mass transfer coefficient value C = 0.0008, slightly lower than the standard range, was found to yield 

mean temperatures relatively close to measured temperatures. The stomatal resistance value and the mass-transfer coefficient 

were both varied in the sensitivity analysis.
4–Testing and Application of Coupled Model with Heat and Salinity Transport 21

Figure 9. Initial groundwater salinity conditions in the Ten Thousand Islands model area. Control structure FU1 

location is shown for reference.

FU1

SALINITY IN PRACTICAL SALINITY UNITS

0 5 10 MILES

0 5 10 KILOMETERS

GULF OF MEXICO

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

EXPLANATION

0 5 10 15 20 25 30

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

CONTROL STRUCTURE
22 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

4.1.2–Existing Conditions Simulation Results

The simulation of existing conditions using the parameters described in section 4.1.1 was evaluated to determine the accuracy 

and reliability of the model application. The flow patterns at structure FU1 were reproduced well, but the model generally underestimated peak flows (fig. 11). The measured flows have a coefficient of variation (standard deviation divided by mean) of 1.15, 

meaning that flow variability was quite high relative to flow magnitude. Table 1 shows the statistics comparing measured flows 

to simulated existing conditions and sensitivity runs at FU1. The difference between results from the existing conditions simulation and the field-measured data is expressed in terms of root-mean-square error (RMSE) and percent explained variance (PEV). 

Figure 10. Comparison of simulated existing conditions flows at the I-75 northern canal boundaries and 

measured flows at structure FU1

0

20

40

60

80

100

1998 1999 2000 2001 2002 2003 2004 2005

YEAR

FLOW, IN CUBIC METERS PER SECOND

EXPLANATION

SIMULATED FLOWS FOR

EXISTING CONDITIONS AT

NORTHERN BOUNDARY

FU1 MEASURED

Figure 11. Comparison of measured and simulated flows at FU1. PSRP is Picayune Strand Restoration Project.

0

10

20

30

40

50

60

70

YEAR

1998 1999 2000 2001 2002 2003 2004 2005

FLOW, IN CUBIC METERS PER SECOND

MEASURED

EXPLANATION

SIMULATED PSRP CONDITIONS

SIMULATED EXISTING CONDITIONS
4–Testing and Application of Coupled Model with Heat and Salinity Transport 23

Table 1. Flow simulation statistics at structure FU1 for existing conditions and sensitivity runs

[Values in bold indicate improvement over existing conditions. RMSE, root mean square error; PEV, percent explained 

variance. Units: m, meter; m3/s, cubic meter per second; %, percent]

Sensitivity run Statistic Flow Stage

Headwater Tailwater

Existing conditions

RMSE 8.39 m3/s 0.16 m 0.10 m

PEV 63.9 % 57.2 % 68.9 %

No northern inflows

RMSE 13.5 m3/s 0.20 m 0.13 m

PEV 23.1 % 31.1 % 52.8 %

Dispersion coefficient plus 50%

RMSE 8.40 m3/s 0.16 m 0.10 m

PEV 63.8 % 57.5 % 68.8 %

Dispersion coefficient minus 50%

RMSE 8.41 m3/s 0.16 m 0.10 m

PEV 63.8 % 57.3 % 69.2 %

Atmospheric mass transfer plus 50%

RMSE 8.40 m3/s 0.16 m 0.10 m

PEV 63.8 % 57.3 % 68.9 %

Atmospheric mass transfer minus 50%

RMSE 8.40 m3/s 0.16 m 0.10 m

PEV 63.8 % 57.2 % 68.9 %

Stomatal resistance plus 50%

RMSE 8.24 m3/s 0.16 m 0.10 m

PEV 64.2 % 54.8 68.3 %

Stomatal resistance minus 50%

RMSE 8.60 m3/s 0.14 m 0.10 m

PEV 63.5 % 62.5 % 69.5 %

Offshore water levels raised 0.15 m

RMSE 8.36 m3/s 0.16 m 0.15 m

PEV 63.9 % 55.2 % 74.6 %

Manning’s n plus 0.05

RMSE 8.47 m3/s 0.16 m 0.12 m

PEV 63.8 % 56.5 % 61.2 %

Manning’s n minus 0.05

RMSE 8.41 m3/s 0.16 m 0.10 m

PEV 63.2 % 59.7 % 71.6 %

No salinity 

RMSE 8.41 m3/s 0.16 m 0.10 m

PEV 63.8 % 58.1 % 68.8 %

No wind

RMSE 8.28 m3/s 0.16 m 0.10 m

PEV 64.5 % 57.2 % 69.0 %

Aquifer conductivity times 10

RMSE 8.43 m3/s 0.15 m 0.10m

PEV 63.9% 59.2% 69.1%

Aquifer conductivity divided by 10

RMSE 8.25 m3/s 0.15 m 0.10m

PEV 65.2% 61.6% 69.4%

No groundwater leakage

RMSE 8.46 m3/s 0.16 m 0.15 m

PEV 59.2% 59.9% 45.1%
24 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

The PEV is a measure of the amount of measured-data variance that is explained by the simulated data and is computed using 

the Nash-Sutcliffe equation:

×

, (22)

where r

2 is the residual variance, and m2 is the measured data variance. For the simulation period, the model reproduced measured 

flows at FU1 with a mean error of -2.75 m3/s and, as shown in table 1, an RMSE of 8.39 m3/s and a PEV of 63.9 percent. The 

high variability of the flow regime represented by the coefficient of variation probably increases the difficulty in reproducing the 

measured data variance with a numerical model.

Continuous water-level data collected at structure FU1 were used for comparison and evaluation of the simulation. Figure 12 

shows measured and simulated existing headwater and tailwater levels at FU1. The largest discrepancy is for headwater stage 

below the weir elevation of 0.20 m NAVD 88. A substantial drawdown in the measured headwater occurs during dry times, and 

it is not represented in the simulation. Leakage through the weir structure would explain this phenomenon, and visual accounts 

have verified the prevalence of leaks around the sides of structure FU1. The statistics comparing measured values to simulated 

existing conditions for headwater and tailwater levels are shown in table 1. The PEV values are similar to the flow statistics at 

FU1, and both RMSE values are on the order of 0.1 m.

Figure 12. Comparison of measured and simulated stage at FU1. PSRP is Picayune Strand Restoration Project.

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1.0

YEAR

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1.0

STAGE, IN METERS NAVD 88 STAGE, IN METERS NAVD 88

YEAR

MEASURED

EXPLANATION

SIMULATED PSRP CONDITIONS

SIMULATED EXISTING CONDITIONS

MEASURED

EXPLANATION

SIMULATED PSRP CONDITIONS

SIMULATED EXISTING CONDITIONS

HEADWATER

TAILWATER
4–Testing and Application of Coupled Model with Heat and Salinity Transport 25

Limited groundwater-level data exist in the simulation area for comparison to simulated values. Because minimal vertical head 

gradient exists in the shallow aquifer, the measured groundwater levels were taken as the average for the three layers computed 

in the model. The groundwater monitoring wells shown in figure 7 were not installed until October 2003, and consequently, the 

coincident period with the simulation is a little over 2 years. 

The simulated groundwater heads do not show as much variability as the measured values (fig. 13), and the inability to represent 

peaks and troughs in groundwater head can result from local phenomena at the wells, such as interaction with surface water or, possibly, aquifer characteristics. Figure 13 shows the computed surface-water stage at the well locations, and the land surface is dry most 

of the time at these sites. However, the greatest discrepancies between measured and computed groundwater head generally occur 

during periods of surface-water inundation (fig. 13). A comparison of measured and computed groundwater heads yield PEV values 

of 43.2 percent for well SGT2W2, 45.1 percent for SGT2W6, and 47.7 percent for SGT4W1. However, if the computed groundwater 

values are replaced by the average of the groundwater and surface-water levels at times of surface-water inundation, the PEV values 

increase to 55.0 percent, 52.4 percent, and 51.2 percent, respectively. The degree to which the water levels in the monitoring wells are 

directly influenced by surface-water inundation is unknown. Another possibility is that the aquifer hydraulic conductivity or the ratio 

of conductivity to storativity (diffusivity) is not correct in the model input. However, conductivity is one of the parameters that was 

allowed to vary in the sensitivity analysis, and varying the hydraulic conductivity did not improve overall simulation results.

Figure 13. Comparison of measured and simulated water levels for selected sites in the Ten 

Thousand Islands area.

LAND SURFACE

0

0.5

1.0

1.5

2.0

2.5

3.0

SGT2W6

2003 2004 2005

YEAR

WATER LEVEL, IN METERS ABOVE NAVD 88

SGT2W2

0

0.5

1.0

1.5

2.0

2.5

3.0

SGT4W1

0

0.5

1.0

1.5

2.0

2.5

3.0 LAND SURFACE

LAND SURFACE

MEASURED

SIMULATED GROUNDWATER

SIMULATED SURFACE WATER

EXPLANATION

MEASURED

SIMULATED GROUNDWATER

SIMULATED SURFACE WATER

EXPLANATION

MEASURED

SIMULATED GROUNDWATER

SIMULATED SURFACE WATER

EXPLANATION
26 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Accurately simulating salinity and temperature in a model is an important step in characterizing the ecological effects associated with hydrologic change. Figure 14 shows the location of salinity and temperature stations used for comparison in the TTI 

domain; data are collected at irregular intervals at these stations. The comparison in figure 15 indicates good agreement between 

simulated and measured existing salinity fluctuations at Port of the Islands, RK451, TTI67, and TTI75, and fair agreement 

at RK457, TTI70, and TTI76. Headwater salinities at FU1 are very low, and BC19 was frequently dry during the simulation 

period, making comparisons difficult. Goodness of fit statistics for salinity data were computed for three of these sites: Port of 

the Islands, RK451, and TT176 (table 2). Port of the Islands salinity matches with an RMSE of 7.56 and a PEV of 71.0 percent. 

The RK451 salinity matches with an RMSE of 3.22 and a PEV of 46.4 percent. Although RK451 salinity has a lower associated 

RMSE than Port of the Islands, this is due primarily to the lower salinity variability at RK451 compared to Port of the Islands; 

the lower PEV indicates that the fit at RK451 is not as good as at Port of the Islands. TTI76 salinity station matches with an 

RMSE of 6.51 and a PEV of 66.9 percent. The goodness of fit for TTI76 salinity is, therefore, apparently between that for Port 

of the Islands and RK451. A number of factors affect the ability of the model to accurately reproduce salinity values at different 

locations, including the spatial and temporal variability in salinity, proximity to boundaries, and accuracy of the flow dynamics.

Figure 14. Surface-water salinity and temperature stations in the Ten Thousand Islands area. 

Picayune Strand

Restoration Project area

0 5 10 MILES

0 5 10 KILOMETERS

Fakahatchee

Strand

TTI70

FU1

TTI75

TTI67

TTI76

RK457 BC19

RK451

EXPLANATION

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

41

75

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

Barron River

Faka

Union Canal

Port of the Islands FU-1

Miller Canal

Merritt Canal

Prairie Canal

Faka Union Canal

UPSTREAM

SALINITY AND TEMPERATURE

MEASUREMENT STATION

CONTROL STRUCTURE

TTI67

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2855000

2865000

2875000

2885000

2895000

GULF OF MEXICO
4–Testing and Application of Coupled Model with Heat and Salinity Transport 27

Figure 15. Comparison of measured and simulated salinity for selected sites in the Ten Thousand Islands area. PSRP is 

Picayune Strand Restoration Project.

PORT OF THE ISLANDS

0

0

0

5

5

5

10

10

10

15

15

15

20

20

20

25

25

25

30

30

30

35

35

35

40

40

40

45

45

45

1998 1999 2000 2001 2002 2003 2004 2005

FU1 UPSTREAM

RK457

SALINITY, IN PRACTICAL SALINITY UNITS

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

RK451

BC19

0

5

10

15

20

25

30

35

40

45

1998 1999 2000 2001 2002 2003 2004 2005

YEAR

0

5

10

15

20

25

30

35

40

45

EXPLANATION

SIMULATED EXISTING CONDITIONS

MEASURED

SIMULATED PSRP CONDITIONS
28 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 15. Comparison of measured and simulated salinity for selected sites in the Ten Thousand Islands area. PSRP is 

Picayune Strand Restoration Project.—Continued

0

5

10

15

20

25

30

35

40

45

0

5

10

15

20

25

30

35

40

45

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

TTI75

TTI70

TTI76

TTI67

0

5

10

15

20

25

30

35

40

45

0

5

10

15

20

25

30

35

40

45

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

YEAR

SALINITY, IN PRACTICAL SALINITY UNITS

EXPLANATION

SIMULATED EXISTING CONDITIONS

MEASURED

SIMULATED PSRP CONDITIONS
4–Testing and Application of Coupled Model with Heat and Salinity Transport 29

Table 2. Transport simulation statistics at field sites for existing conditions and sensitivity runs. 

[Salinity is expressed in practical salinity units. Values in bold indicate improved over existing conditions. PEV, percent explained variance; RMSE, root mean 

square error. Units: °C, degrees Celsius; %, percent; m, meter]

Sensitivity run Statistic

Port of the Islands RK451 TTI76

Salinity Temperature Salinity Temperature Salinity Temperature

Existing conditions

RMSE 7.56 2.11°C 3.22 1.74°C 6.51 1.78°C

PEV 71.0 % 72.7 % 46.4 % 88.1 % 66.9 % 90.4 %

No northern inflows

RMSE 14.1 2.13°C 3.00 1.73°C 8.32 1.78°C

PEV 42.2 % 71.2 % 49.6 % 88.1 % 54.4 % 90.3 %

Dispersion plus 50%

RMSE 9.29 2.12°C 2.95 1.73°C 7.63 1.75°C

PEV 68.7 % 72.3 % 52.7 % 88.2 % 62.0 % 90.3 %

Dispersion minus 50%

RMSE 8.00 2.08°C 3.93 1.74°C 5.22 1.83°C

PEV 70.1 % 73.2 % 27.4 % 87.9 % 67.9 % 90.6 %

Atmospheric mass transfer plus 50%

RMSE 7.57 2.34°C 3.21 1.89°C 6.55 2.13°C

PEV 70.9 % 75.9 % 46.8 % 88.9 % 66.8 % 92.6 %

Atmospheric mass transfer minus 50%

RMSE 7.56 1.94°C 3.22 1.74°C 6.50 1.59°C

PEV 70.9 % 67.2 % 46.4 % 86.5 % 67.2 % 86.4 %

Stomatal resistance plus 50%

RMSE 7.52 1.85°C 3.42 1.63°C 6.15 1.45°C

PEV 70.0 % 69.3 % 44.3 % 88.2 % 66.4 % 89.2 %

Stomatal resistance minus 50%

RMSE 7.75 3.11°C 2.99 2.24°C 7.32 2.65°C

PEV 72.1 % 76.4 % 48.9 % 87.7 % 64.4 % 91.6 %

Offshore water levels raised 0.15 m

RMSE 8.53 2.11°C 2.92 1.72°C 7.76 1.72°C

PEV 70.1 % 72.2 % 54.4 % 88.2 % 63.1 % 90.7 %

Manning’s n plus 0.05

RMSE 7.37 2.10°C 3.67 1.73°C 5.91 1.79°C

PEV 71.8 % 72.6 % 34.9 % 88.2 % 69.0 % 91.0 %

Manning’s n minus 0.05

RMSE 7.82 2.12°C 3.11 1.74°C 6.69 1.79°C

PEV 69.0 % 72.7 % 49.0 % 88.1 % 65.7 % 90.2 %

No salinity 

RMSE — 2.11°C — 1.74°C — 1.79°C

PEV — 72.8 % — 88.1 % — 90.4 %

No wind

RMSE 7.52 2.10°C 3.30 1.73°C 6.46 1.79°C

PEV 71.0 % 72.7 % 43.4 % 88.1 % 66.1 % 90.5 %

Aquifer conductivity times 10

RMSE 7.52 2.11°C 3.31 1.74°C 7.08 1.78°C

PEV 71.6% 72.8% 46.4% 88.1% 59.8% 90.4%

Aquifer conductivity divided by 10

RMSE 7.50 2.11°C 3.46 1.73°C 6.22 1.78°C

PEV 72.0 % 72.8% 41.3% 88.1% 67.6 % 90.4%

No groundwater leakage

RMSE 7.88 2.11°C 4.55 1.74°C 5.70 1.79°C

PEV 66.6% 73.6% 12.4% 88.2% 60.2% 90.5%
30 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

The heat transport capabilities of the model were evaluated by comparing simulated temperatures to those measured in 

the field (fig. 16). Temperature is a more spatially uniform property, and consequently, the associated matches tend to be closer 

than for salinity. TTI76 temperatures have an associated RMSE of 1.78 oC and a PEV of 90.4 percent, and RK451 temperatures 

have an RMSE of 1.74 oC and a PEV of 88.1 percent (table 2). These statistics indicate better agreement than those for Port of 

the Islands temperatures, which have an associated RMSE of 2.11 oC and a PEV of 72.7 percent. This difference may be related 

to the deeper waters at Port of the Islands, where substantial vertical variations in temperature and salinity have been observed. 

Figure 17 shows measured time-series temperature data at the Port of the Islands from two probes at depths of 0.5 and 3.0 m 

below the water surface. The model-computed temperatures for existing conditions are mostly lower than measured values, at 

least for this time period. The deviation between the two measured temperature datasets indicates that, compared to shallow 

water, deeper water tends to remain warmer longer during periods of overall cooling. Vertical temperature stratification in Port 

of the Islands may limit the ability to simulate this feature with a two-dimensional model.

4.2–Sensitivity Analysis

To determine the importance of parameters used in the application and to evaluate the effect of parameter uncertainty, a 

sensitivity analysis was performed on the calibrated existing conditions application. As noted earlier, an iterative parameterestimation technique was not used because only a few parameters were adjusted during model calibration. As with iterative 

parameter estimation techniques, the sensitivity analysis performed here examined the effect of perturbing the model parameters 

by a representative amount from their calibrated values. The parameters chosen for testing include the calibration parameters 

(northern inflows, dispersion coefficient and atmospheric mass-transfer coefficient) and other parameters of importance (stomatal resistance, tidal water levels, frictional coefficient, salinity density effects, wind forcing, leakage coefficient, and aquifer 

hydraulic conductivity).

The comparison of measured and computed flow and water levels at FU1 yielded several lower RMSE or higher PEV

statistics (table 1), which individually indicate a better match to measured data. The most substantial improvements at FU1 

occurred when the aquifer conductivity was reduced by an order of magnitude. The greatest decrease in statistical fit occurred 

when the northern boundary inflows were removed, which is not unexpected for conditions at FU1.

The statistical comparison of measured and computed salinity and temperature at the Port of the Islands, RK451, and 

TTI76 (table 2) shows that none of the parameter perturbations yield an overall improvement at all sites, although some more 

localized phenomena are of interest. For example, the statistics for Port of the Islands improved when aquifer conductivity was 

reduced one order of magnitude, matching the improvement in flow statistics at nearby FU1 under the same conditions (table 1). 

This result indicates that the conductivity assigned to the aquifer system may be too high near the canal system associated with 

FU1 and the Port of the Islands. The leakage coefficient for the canal also affects this interaction; however, and the statistics for 

RK451 (table 2) do not suggest that lowering aquifer conductivity improves results regionally. Increasing the dispersion coefficient by 50 percent improved all statistics at RK451, but degraded all statistics at Port of the Islands (table 2); the same pattern 

occurred when the tidal levels were increased 15 cm. If increasing both dispersion and tidal levels improves results at RK451, 

but not elsewhere, then the connectivity of the area around RK451 (Rookery Bay) is most likely a factor. This conclusion is reinforced by the fact that reducing the Manning’s n frictional resistance improves the salinity statistics at RK451, with no change 

in temperature statistics, and degrades the statistics at the other locations (table 2). The Rookery Bay area was difficult to define 

topographically in the model, because it has numerous small channels and islands that are smaller than individual model cells. 

Localized lowering of the model topography or decreasing of frictional resistance could improve the results in Rookery Bay, but 

the improvements shown in table 2 are minimal. 

Further insight into model parameters can be made from the simulation without leakage statistics for salinity and temperature 

given in table 2. Salinity statistics were not improved, and are in fact substantially worse at Port of the Islands and RK451, but 

the temperature PEV improved slightly at all three sites. Although the RMSE does not reflect this trend, and surface-water temperatures are relatively insensitive to groundwater leakage, using a previous daily average air temperature for groundwater temperature may not be a valid assumption. In fact, field data indicate that groundwater is typically cooler than surface water during 

the summer, with the opposite being true in winter (Harvey and others, 2000). Thus, a more complex approach for simulating 

groundwater heat transport may be justified.

The sensitivity analysis does not indicate any specific improvement in the calibrated model based on changes to the 

parameters examined. Spatially varying changes in parameters could be implemented with more information. The uncertainty 

in the northern boundary inflows must be considered the most important, because it (1) highly affects the response parameters, 

(2) makes the relative contribution of leakage more difficult to define, and (3) was calibrated with model results at the structure 

FU1. Wind forcing appears to have a smaller effect in the TTI application compared to Everglades National Park applications 

(Swain and others, 2004), which is probably due to the greater degree of channelization in the TTI area.
4–Testing and Application of Coupled Model with Heat and Salinity Transport 31

Figure 16. Comparison of measured and simulated water temperature for selected sites in the Ten Thousand Islands 

model domain. PSRP is Picayune Strand Restoration Project.

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

PORT OF THE ISLANDS

RK457

WATER TEMPERATURE, IN DEGREES CELSIUS

FU1 UPSTREAM

10

15

20

25

30

35

1998 1999 2000 2001 2002 2003 2004 2005

RK451

BC19

1998 1999 2000 2001 2002 2003 2004 2005

YEAR

10

15

20

25

30

35

10

15

20

25

30

35

10

15

20

25

30

35

10

15

20

25

30

35

EXPLANATION

SIMULATED EXISTING CONDITIONS

MEASURED

SIMULATED PSRP CONDITIONS
32 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 16. Comparison of measured and simulated water temperature for selected sites in the Ten Thousand Islands area. 

PSRP is Picayune Strand Restoration Project.—Continued

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

TTI75

TTI70

WATER TEMPERATURE, IN DEGREES CELSIUS

TTI67

TTI76

1998 1999 2000 2001 2002 2003 2004 2005

1998 1999 2000 2001 2002 2003 2004 2005

YEAR

10

15

20

25

30

35

10

15

20

25

30

35

10

15

20

25

30

35

10

15

20

25

30

35

EXPLANATION

SIMULATED EXISTING CONDITIONS

MEASURED

SIMULATED PSRP CONDITIONS
4–Testing and Application of Coupled Model with Heat and Salinity Transport 33

4.3–Simulation of Picayune Strand Restoration Project Scenario

The PSRP simulation incorporated parameters identical to those used in the calibrated existing conditions simulation, with 

the exception of changes outlined in section 3.6. The same time period was simulated with the same offshore and inland boundary conditions. This allowed a direct comparison between the existing conditions and PSRP simulations on a time step-by-time 

step basis.

The flow conditions at structure FU1 strongly affect conditions at the Port of the Islands, and the PSRP simulation indicates 

that the change in flow at FU1 resulting from PSRP modifications is minimal (fig. 11). Some reduction of peak FU1 flows in the 

PSRP simulation is indicated, with higher sporadic flows occurring at times when the existing conditions simulation has low 

steady flows. The comparison of FU1 water levels shown in figure 12 indicates that in the PSRP simulation, headwater stage 

drops much more during dry times and the tailwater seems to have greater fluctuations at time scales less than 1 month, compared to the existing conditions simulation. Statistics comparing the simulation of the PSRP scenario with the existing conditions simulation at FU1 indicate only a 2-percent drop in PSRP simulated flow (table 3). This result may indicate that substantial 

water still enters the canal immediately north of FU1, even when the canal network is plugged. However, the root-mean square 

difference (RMSD) of 4.14 m3/s indicates substantial changes in flow timing, consistent with the rerouting of water in the restoration plan. The stage results in table 3 help explain the changes induced by the PSRP. Although average headwater stage at 

FU1 is 8.5 cm lower in the PSRP simulation than the existing conditions simulation, the average tailwater stage is 0.9 cm higher. 

The flow reduction at the structure due to canal plugging reroutes water across the adjacent wetlands, and actually elevates the 

average tailwater stage slightly.

Several PSRP effects on salinity are apparent in figure 15: (1) increased peak salinities at Port of the Islands; (2) marked 

salinity intrusion upstream from FU1; (3) substantial decreases in salinity at RK451, RK457, and TTI75, especially during drier 

periods; (4) minor reductions in salinity at TTI70 and TTI76; and (4) minimal change in salinity at TTI67 and BC19. This spatial 

pattern of salinity change indicates that the redistributed freshwater affects areas west of Faka Union Canal the most, and areas 

east of the canal the least. The reduction of freshwater flow to the Port of the Islands yields periods of increased salinity, but has 

little effect on locations near the coastline (TTI67) and further inland in the wetlands (BC19). Based on the analysis in section 

4.3, these salinity differences can be considered large enough to be an important result of the model.

Figure 17. Comparison of measured and simulated water temperature for the Port of the Islands site.

EXPLANATION

2004 2005

MEASURED, 0.5-m DEPTH BELOW WATER SURFACE

MEASURED, 3.0-m DEPTH BELOW WATER SURFACE

SIMULATED EXISTING CONDITIONS (VERTICALLY AVERAGED)

10

14

18

22

26

30

WATER TEMPERATURE, IN DEGREES CELSIUS

DECEMBER JANUARY FEBRUARY MARCH APRIL
34 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Surface-water temperature differences between existing conditions and the PSRP simulation are minor at the majority 

of sites shown in figure 16. One small but important pattern is evident for Port of the Islands on dates of minimum simulated 

annual temperatures, particularly January 7, 1999, January 9, 2002, and December 21, 2003. The PSRP minimum temperatures 

are not as low as the existing conditions temperatures, and are sometimes warmer by over 1.5 °C. This is important to certain 

marine life, such as manatees, which suffer from cold stress at temperatures below 20 °C. The slightly higher temperatures in 

the simulated PSRP conditions compared to existing conditions for December 21, 2003, are related to the heating of additional 

ponded water in areas that are dry under existing conditions (fig. 18).

The salinity and temperature comparison statistics in table 3 indicate that, in the PSRP simulation, average salinity increases 

at Port of the Islands and decreases at RK451 and TTI76, as would be expected from the PSRP rerouting of water. Because the 

deviation in salinity patterns between the existing conditions and PSRP simulations only occurs part of the time (fig. 15), the 

average deviations listed in table 3 are small. The average Port of the Islands temperature is actually 0.03 oC lower in the PSRP 

simulation, despite the maintenance of higher temperatures during cold periods just described.

4.4–Ecological Comparison: Manatee Locations

The simulation of salinity and heat transport in the existing conditions and PSRP scenarios yields information that can 

be used to examine the effect of hydrology on manatee response, and as input for models of manatee behavior. Aerial surveys 

are used for tracking manatee locations as well as satellite telemetry data from transmitters attached to the manatees (Stith and 

others, 2004). Temperature and salinity outputs from the TTI existing conditions simulation (fig. 19) were compared with the 

manatee tracking data for November 16 and 22, 2002. The affinity of manatees for warmer and fresher water in the Port of the 

Islands area can be inferred from the tracked locations. Although only 6 days apart, the substantially colder water temperatures 

on November 22, 2002, make the warmer inland waters around the Port of the Islands preferable to manatees.

The hydrologic conditions produced by the existing conditions and PSRP simulations can be used in the development of 

the individually based manatee behavior model (Stith and others, 2006). To couple the hydrologic and manatee behavior models, 

temperature and salinity conditions at specified locations of interest (nodes) along the manatee travel corridors (fig. 20) must be 

extracted from the existing conditions and PSRP simulations. The differences between existing and restoration conditions for the 

nodes in figure 20 are summarized statistically in table 4. Average salinity decreases at most locations, except the Faka Union 

Canal Basin and the Collier Seminole Basin. Inland sites tend to be warmer, as with Faka Union River and SR 92 Canal. The 

analysis of the model’s ability to represent changes presented earlier indicates that the majority of the differences in table 4 are 

statistically significant.

Table 3. Simulation statistics at field sites comparing Picayune Strand 

Restoration Project-implemented conditions with existing conditions

[Salinity expressed in practical salinity units. PSRP, Picayune Strand Restoration Project; 

RMSD, root mean square difference; NA, not applicable. Units: m, meter; m3/s, cubic meter 

per second; %, percent]

Measurement Average 

difference RMSD Percent difference

FU1

Flow -0.16 m3/s 4.14 m3/s -2.0%

Headwater stage -0.085 m 0.174 m NA

Tailwater stage 0.009 m 0.061 m NA

Port of the islands

Salinity 1.25 4.52 NA

Temperature -0.03 °C 0.38 °C NA

RK451

Salinity -2.26 2.99 NA

Temperature -0.01 °C 0.06 °C NA

TT176

Salinity -1.98 3.24 NA

Temperature 0.01 °C 0.12 °C NA
4–Testing and Application of Coupled Model with Heat and Salinity Transport 35

Figure 18. Simulated surface-water temperature for December 21, 2003, 

for A, existing conditions, and B, Picayune Strand Restoration Projectimplemented conditions in the Ten Thousand Islands model area.

11 13 15 17 19 21 23 25 27 29 31 33 35

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM NORTHING, IN METERS

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

WATER TEMPERATURE, IN DEGREES CELSIUS

UTM EASTING, IN METERS

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

EXPLANATION

B

A

0 5 10 KILOMETERS

0 5 10 MILES
36 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Figure 19. Manatee locations and measured surface-water salinity and temperature on November 16 and 22, 2002, in the Ten 

Thousand Islands model area.

11 13 15 17 19 21 23 25 27 29 31 33 35

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34

Nov. 16, 2002 Nov. 22, 2002

SALINITY SALINITY

TEMPERATURE TEMPERATURE

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

2895000

2890000

2885000

2880000

2875000

2870000

2865000

2860000

2855000

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

WATER TEMPERATURE, IN DEGREES CELSIUS

SALINITY, IN PRACTICAL SALINITY UNITS

EXPLANATION

MANATEE LOCATION

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000 420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

EXPLANATION Nov. 16, 2002 Nov. 22, 2002

MANATEE LOCATION

0 5 10 KILOMETERS

0 5 10 MILES

0 5 10 KILOMETERS

0 5 10 MILES
4–Testing and Application of Coupled Model with Heat and Salinity Transport 37

Figure 20. Salinity of temperature nodes for manatee travel corridors in the model area. Node number references 

are presented in table 4.

EXPLANATION

UTM EASTING, IN METERS

UTM NORTHING, IN METERS

41

75

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator Zone 17

Barron River

Fakahatchee

Strand

Faka

Union

Canal

FU-1

Miller Canal

Merritt Canal

Prairie Canal

Picayune Strand

Restoration Project area

Faka Union Canal

420000 425000 430000 435000 440000 445000 450000 455000 460000 465000

2855000

2865000

2875000

2885000

2895000

0 5 10 MILES

0 5 10 KILOMETERS

NODE WHERE SALINITY AND TEMPERATURE ARE

DEFINED FOR MANATEE TRAVEL CORRIDOR

CONTROL STRUCTURE

3

GULF OF MEXICO

6

2

12

14

11

18

10

25

19

23

17 15

4

1

3

24

20

22

7

16 21 9

8

5

13
38 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Table 4. Simulation statistics for selected nodes in Manatee travel corridors comparing Picayune 

Strand Restoration Project-implemented conditions with existing conditions.

[Node locations are shown by reference number in figure 20. Salinity is expressed in practical salinity units. PSRP, 

Picayune Strand Restoration Project; RMSD, root mean square difference]

Location of interest (node) Salinity Temperature

(degrees Celsius)

Reference 

number Name Average 

difference RMSD Average difference RMSD

1 Blackwater Bay -3.55 3.24 -0.014 0.110

2 Blackwater River -7.98 7.04 0.015 0.355

3 Buttonwood Bay -2.93 2.94 -0.015 0.124

4 Cape Romano -0.52 0.49 -0.004 0.015

5 Collier Seminole Basin 0.78 1.39 -0.659 4.078

6 East River -2.38 3.49 0.031 0.201

7 Faka Union Bay -1.17 2.39 -0.004 0.130

8 Faka Union Canal Basin 1.25 4.34 -0.030 0.381

9 Faka Union River -0.53 2.90 0.742 1.456

10 Fakahatchee Bay -1.84 2.41 0.010 0.104

11 Fakahatchee Pass -0.39 0.53 -0.002 0.021

12 Fakahatchee River -2.71 3.55 0.033 0.243

13 Ferguson Bay -1.59 1.94 0.006 0.077

14 Ferguson River -1.73 2.10 0.007 0.090

15 Goodland Bay -3.22 2.68 -0.016 0.084

16 Little Wood River -2.52 3.52 -0.007 0.208

17 Marco Island Basin -0.37 0.38 -0.009 0.038

18 No Name River -1.76 2.17 0.007 0.089

19 Palm Bay -4.47 3.99 -0.007 0.135

20 Pumpkin Bay -1.78 2.38 -0.009 0.113

21 Pumpkin River -1.91 2.55 -0.010 0.127

22 Santina Bay -1.35 2.16 -0.005 0.103

23 SR 92 Canal -0.34 1.68 1.602 1.579

24 Whitney River -1.87 2.41 -0.009 0.119

25 Wood River -0.21 3.45 -0.012 0.260
4–Testing and Application of Coupled Model with Heat and Salinity Transport 39

4.5–Model Limitations and Capability to Predict Restoration Changes

The TTI model has several limitation which must be considered when applying the model and utilizing the results. The 

boundary conditions are specified time series of water-level, flow, or salinity values, so the model cannot represent any new 

scenario that would cause these parameters to differ from user-input values. The tidal levels and salinity are the boundary conditions that can be considered most realistically unaffected by systematic changes, as the open ocean provides a very large source 

and sink. However, at locations where the tidal boundary is close to the shoreline, such as the northwestern corner of the model 

domain (fig. 1), changes in near-shore salinity due to modified coastal outflows, for example, may not be properly represented 

due to boundary dominance. The inland surface-water flow boundaries represent water delivered across I–75 (fig. 1), and could 

differ with changes to the downstream system. Different water-level gradients should require different boundary flow rates. 

However, if a new simulation assumes a similar amount of water available north of I–75, and it is not diverted elsewhere, using 

the same flow-rate time series could be an acceptable approximation. The groundwater head boundaries on the northern and 

eastern sides of the model also should be affected by different surface-water conditions and could require modifications in a new 

scenario, perhaps using model-generated surface-water levels as a guide.

The representation of the canal system in the TTI model with model cells is an approximation with limitations that must be 

considered. As the model cells are much larger than the canals, values of evapotranspiration are computed for a larger inundated 

area than the actual canal. Rainfall is also applied for the whole cell area, although this can be seen as the rain on either side of 

the canal being put directly in the canal. The land-surface elevation of the model cells are higher than the canal bottoms, so the 

cells can become dry at water levels that would still allow actual canal flow. This does not occur in the simulations described in 

this report, but the possibility exists at very low water levels.

The limitations in model parameter accuracy and output uncertainties is discussed thoroughly in sections 4.1 and 4.2. The 

ability of the model to reproduce measured conditions is sufficient to develop insight into hydrologic processes in the TTI area, 

but the model cannot be expected to produce parameter values more exact than the error bounds indicate. One of the primary 

limitations noted earlier is the assumed temperature of the surface-water inflows and the groundwater. Particularly, the effect of 

groundwater leakage on surface-water temperature is not well represented, although this seems to be a minimal effect.

Using the current model to examine how restoration efforts affect hydrology in the TTI area requires simulating the PSRP 

scenario and comparing the results to the calibrated simulation of existing conditions. Often the ability of a model to predict 

differences between simulations is equated to its ability to match field data (Nandakumar and Mein, 1997), but a comparison of 

two model simulations from the same calibration does not yield the same errors as a comparison of a model simulation with collected data. In fact, the model errors tend to be similar in different simulations, so errors cancel out when taking the difference 

between the calibrated existing conditions and the PSRP scenario.

The difference between the TTI model existing conditions simulation and the PSRP simulation can be visualized as the difference in the mean of two model-generated datasets. A statistical Z-test can be used to determine if the difference in the mean 

of two datasets is significant. The uncertainty in each dataset is characterized by a distribution about the mean of the dataset. The 

RMSE of the existing condition simulation fit to measured data expresses the standard deviation of the uncertainty distribution. 

The Z-statistic is calculated as follows (Ott, 1993, p. 220):

μ

, (23)

where μdiff is the difference in means of the two datasets, σ is the standard deviation (the RMSE), and n is the number of datapoints. Using a Z-statistic of 1.28 for a 90-percent confidence limit (Ott, 1993, p. A–3), a salinity RMSE of σ = 7.56 at the Port 

of the Islands from table 2 and n = 11,680 for the 8-year simulation using 6-hour output yields a significant salinity μdiff value 

of 0.09. This value is the minimum significant difference in average salinity, over the entire simulation period, that the model 

can predict between scenarios with 90-percent confidence. Application of equation 23 to the temperature prediction at Port of 

the Islands using the RMSE of σ = 2.11 oC from table 2 yields a significant temperature μdiff value of 0.025 oC. This method 

indicates that the model can predict differences between scenarios that are only Z n / = 1.2 percent of the error in matching 

field data.

There are some assumptions associated with using the Z-test to evaluate the ability of the model to predict differences 

between scenarios. First, the restoration simulation is assumed to have the same RMSE for predicted quantities as the existing 

conditions simulation. If the differences between the existing conditions and restoration simulations are small, this is typically a 

valid assumption. Second, the Z-test statistic is based on a normal distribution. The distribution of salinity error for Port of the 

Islands data does not pass a Kolmogorov-Smirnov test for normality at a 90-percent confidence level. For a sample size greater 

than 30 to 40, however, the Z-test is considered representative, even for a non-normal distribution.
40 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

5–Summary 

A coupled surface-water/groundwater simulation was created for the TTI area using the FTLOADDS model code. This application includes the simulated transport of salinity and heat in the surface-water system and salinity in the groundwater system, making it 

useful for examining the ecological effects of hydrologic change. Recent advancements in the heat-transport algorithm allow a choice 

of techniques to represent latent heat and evapotranspiration. The representation of flow at hydraulic barriers has been simplified and 

was incorporated into the solution matrix, because weir flow into the Port of the Islands area is an important feature of the model area.

The simulation input dataset utilized a number of different sources. Data from NOAA and LABINS were combined to 

create the grid of offshore bathymetry and inland topography at 500 × 500-m spacing. The canal system was represented by 

computing a composite frictional resistance term for canals and altering the grid topography so that the model cells have the 

proper storage volume. The flow coefficients for surface-water structures were determined from field data, and atmospheric variables and parameters for heat transport, rainfall, wind friction, and evapotranspiration were obtained or derived from the FAWN 

network and results of the TIME model application of FTLOADDS. Boundary tidal levels were obtained from the NOAA Tides 

and Currents Portal, and inland flows were initially equated to measured flows at weir structure FU1 and altered for calibration.

The simulation was executed for the 1998–2005 period, and parameters having the greatest uncertainty were varied to calibrate the model. The northern boundary surface-water inflows were modified to minimize flow errors at the FU1 site on the Faka 

Union Canal. The dispersion coefficient for salt and heat transport in surface water was set to a value of 100 m2/s to obtain mean 

simulated surface-water salinity values that matched measured values, especially at the Port of the Islands. The mass-transfer 

coefficient for the computation of sensible and latent heat was reduced from a value of 0.0020, used for Everglades National 

Park wetlands, to a value of 0.0008 in order to fit measured surface-water temperatures.

The simulation of existing conditions was evaluated with statistics comparing measured and simulated values of surfacewater flow and stage at structure FU1, salinity and temperature at various surface-water sites, and groundwater levels. Because 

flow at structure FU1 is highly variable, the PEV of 63.9 percent is considered to be a reasonable reproduction of measured flows 

at FU1. Sensitivity analyses were performed by varying the important model input values and examining the output statistics. 

Several model behavior characteristics are evident based on sensitivity analysis results. The fact that increasing tidal levels 

and dispersivity only improved results at Rookery Bay indicates that the connectivity of Rookery Bay to the offshore area is 

probably underrepresented. Neglecting groundwater leakage yielded slightly improved surface-water temperatures statistics, 

which indicates that the average air temperature for the previous day may not be a reasonable estimate of groundwater temperature. This estimated temperature was also used for the northern boundary surface-water inflows, which may be a better assumption than for groundwater temperature, although this has not been determined. The simulation of heat transport in groundwater 

may improve the representation of temperature, and the latest version of SEAWAT, which has this capability, could be incorporated into FTLOADDS. Because no parameter modification was shown to improve overall model performance from the 

calibrated simulation, the sensitivity analyses did not necessitate a modification of the model calibration.

The proposed PSRP changes were implemented in the model and simulation results were compared to those for existing 

conditions. Although flow, water level, salinity, and temperature did not change greatly, substantial trends are evident. Coastal 

wetland areas, especially west of Faka Union Canal, remained inundated longer into the dry season in the PSRP simulation 

compared to the existing conditions simulation. Coastal embayments remained less saline at the beginning of the dry season, 

especially in the western areas of the model domain around Rookery Bay. The Port of the Islands area experienced higher peak 

salinities, and salinity was substantially greater upstream from the FU1 structure. Temperatures did not change greatly, although 

the Port of the Islands temperatures remained a few degrees warmer than existing conditions during periods when annual minimums occurred. Apparently, these higher temperatures occur because the PSRP simulation retains a greater amount of warm, 

standing surface water that enters the Port of the Islands area during these times.

Field telemetry data for manatee locations at specific times were compared with model results for surface-water salinity and 

temperature. The manatee preference for the warmer waters of the Port of the Islands area was apparent during colder periods. The 

presence of fresher water north of FU1 coincides with the movement of the manatees toward FU1 from an offshore area where they 

feed. Statistics were generated at locations of importance in an individually based manatee behavior model that requires temperature 

and salinity values as input to determine manatee migration patterns. The model results indicate that the PSRP modifications reduced 

salinity in most coastal locations, increased salinity in the blocked canal system, and increased water temperature slightly in inland 

areas. These patterns are exemplified, respectively, by (1) a decrease in average salinity of 7.98 in Blackwater River, (2) an increase 

in average salinity of 1.25 in the Faka Union Canal Basin; and (3) an increase in average temperature of 0.74 °C in Faka Union River.

The application of the FTLOADDS code to the TTI area has produced a tool that can be used to evaluate a number of 

factors that affect flows, hydroperiods, salinities, and temperatures in the TTI coastal area. The changes proposed for the PSRP 

are one example of how the application can be used to evaluate the effects of changes to the hydrologic system. The linkage of 

the TTI application with the manatee telemetry data and the individually based manatee behavior model input shows how the 

TTI application can be used to study the ecological effects caused by hydrologic change. The application could also be used to 

generate information for other species and interests that are affected by salinity, temperature, water depth, or flow.
6–References Cited 41

6–References Cited

Barenbrug, A.W.T., 1974. Psychrometry and Psychrometric Charts: 3rd Edition, Cape Town, 484 S.A.: Cape and Transvaal 

Printers Ltd.

Brutsaert, W., 1982, Evaporation into the atmosphere: The Netherlands, Kluwer Academic Publishers.

Decker, Jeremy, and Swain, Eric, 2008. Hydrologic Modeling of South Florida Environmental Parameters and Application to 

Ecology, Salinity, and Heat Transport: 2008 Greater Everglades Ecosystem Restoration Conference, Naples, Florida, July 28–

August 1, 2008, p. 81.

Donders, T.H., Wagner, Friederike, and Visscher, Henk, 2005, Quantification strategies for human-induced and natural hydrological 

changes in wetland vegetation, southern Florida, USA: Quaternary Research, v. 64, p. 333–342

ESRI, 2008, ArcMap: Mapping and Editing Made Easy: Accessed on Nov. 25, 2008, at http://www.esri.com/news/arcnews/

summer99articles/ai8special/ai8_arcmap.html

Florida Department of Environmental Protection, 2008, LABINS: Land Boundary Information System: Accessed on Nov. 21, 2008, 

at http://data.labins.org/2003/

French, R.H., 1985, Open Channel Hydraulics: New York, McGraw-Hill.

Guo, Weixing, and Langevin, C.D., 2002, User’s guide to SEAWAT: A computer program for simulation of three-dimensional variable-density ground-water flow: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A7, 77 p.

Harvey, J.W., Jackson, J.M., Mooney, R.H., and Choi, J., 2000, Interaction between ground water and surface water in Taylor 

Slough and vicinity, Everglades National Park, South Florida: Study methods and appendixes: U.S. Geological Survey Open-File 

Report 00–483, 67 p.

Klein, Howard, 1980, Water-resources investigations, Collier County, Florida: U.S. Geological Survey Open-File Report 

80–1207, 29 p.

Langevin, C.D., Thorne, D.T., Dausman, A.M., Sukop, M.C., and Guo, Weixing, 2008, SEAWAT Version 4: A computer 

program for simulation of multi-species solute and heat transport: U.S. Geological Survey Techniques and Methods book 6, 

chap. A22, 39 p.

Langtimm, Cathy, Reid, James, Slone, Dan, Stith, Brad, Swain, Eric, Doyle, Terry, Snow, Skip, and Butler, Susan, 2006, Effects 

of hydrological restoration on manatees: A research program to integrate data, models and long-term monitoring across the 

Ten Thousand Islands and Everglades, in 2006 Greater Everglades Ecosystem Restoration Conference proceedings, Lake 

Buena Vista, Florida, June 5–9, 2006, p. 126, available online at http://sofia.usgs.gov/geer/2006/index.html

Lee, J.K., and Carter, Virginia, 1999, Field measurement of flow resistance in the Florida Everglades, in Proceedings of the 

South Florida Restoration Science Forum: May 17–19, 1999, Boca Raton, Florida: U.S. Geological Survey Open File Report 

99–181.

Leendertse, J.J., 1987, Aspects of SIMSYS2D, a system for two-dimensional flow computation: Santa Monica, Calif., Rand 

Corporation Report R–3572–USGS, 80 p.

Nandakumar, Nandakumar, and Mein, Russell G., 1997, Uncertainty in rainfall—Runoff model simulations and the implications 

for predicting the hydrologic effects of land-use change: Journal of Hydrology, v. 194, p. 211–232.

National Oceanic and Atmospheric Administration, 2008, Center for Operational Oceanographic Products and Services: 

Accessed on Dec. 3, 2008, at http://tidesandcurrents.noaa.gov/

Neelz, Sylvian, and Pender, Gareth, 2007, Sub-grid scale parameterization of 2D hydrodynamic models of inundation in the 

urban area: Acta Geophysica, v. 55, no. 1, p. 65–72.

Otero, J.M., 1994, Computation of flow through water control structures: West Palm Beach, South Florida Water Management 

District Technical Publication WRE–328, 88 p.

Ott, R.L., 1993, An Introduction to Statistical Methods and Data Analysis: Belmont, CA, Duxbury Press.
42 Development, Testing, and Application of a Coupled Hydrodynamic Surface-Water/Groundwater Model (FTLOADDS) 

Poeter, E.P., Hill, M.C., and Banta, E.R., 2005, UCODE_2005 and six other computer codes for universal sensitivity analysis, 

calibration, and uncertainty evaluation: U.S. Geological Survey Techniques and Methods, book 6, chap. A11.

Schaffranek, R.W., 2004, Simulation of surface-water integrated flow and transport in two dimensions: SWIFT2D user’s 

manual: U.S. Geological Survey Techniques and Methods, book 6, chap. B–1, 115 p.

Schaffranek, R.W., Ruhl, H.A., and Hansler, M.E., 1999, An overview of the Southern Inland and Coastal Systems project of the 

U.S. Geological Survey South Florida Ecosystem Program, in Proceedings of the 3rd International Symposium on Ecohydraulics, Salt Lake City, Utah, July 13–16, 1999.

Shoemaker, W.B., Sumner, D.M., and Castillo, Adrian, 2005, Estimating changes in heat energy stored within a column of 

wetland surface water and factors controlling their importance in the surface energy budget, Water Resources Research, v. 41, 

p. W 10411.

South Florida Water Management District, 2005, Documentation of the South Florida Water Management Model Version 5.5, 

November 2005: Accessed on July 20, 2008, at https://my.sfwmd.gov/portal/page?_pageid=1314,2556275,1314_2554821:1

314_2554861&_dad=portal&_schema=PORTAL

Stith, Brad, Reid, Jim, Butler, Susan, Doyle, Terry, and Langtimm, Cathy, 2004, Predicting the effects of hydrologic restoration 

on manatees along the southwest coast of Florida: U.S. Geological Survey Fact Sheet 2004–3137, 4 p.

Stith, Brad, Swain, Eric, Reid, Jim, and Langtimm, Catherine, 2006, Linking a manatee individual-based model to assess restoration effects in the Everglades and Ten Thousand Islands, in 2006 Greater Everglades Ecosystem Restoration Conference 

proceedings, Lake Buena Vista, Florida, June 5–9, 2006, p. 127, available online at http://sofia.usgs.gov/geer/2006/index.html

Stith, Brad, Swain, Eric, Langtimm, Catherine, and Reid, Jim, 2007, Integrating water temperature into hydrologic models to 

assess restoration effects in southwest Florida: Manatees as mobile temperature probes, in Second National Conference on 

Ecosystem Restoration, April 23–27, 2007, Kansas City Missouri, p. 334.

Swain, E.D., 2005, A model for simulation of surface-water integrated flow and transport in two dimensions: User’s guide for 

application to coastal wetlands: U.S. Geological Survey Open-File Report 2005–1033, 88 p.

Swain, E.D. and Decker, J. D., 2007, Developing a Heat-Transport Formulation for a Two-Dimensional Hydrodynamic Model of 

Coastal Wetlands: Second National Conference on Ecosystem Restoration, April 23–27, 2007, Kansas City Missouri p. 340.

Swain, E.D., and Stith, B.M., 2006, Numerical modeling of heat and salinity transport for West Indian Manatee habitats in 

Southwest Florida, in 2006 Greater Everglades Ecosystem Restoration Conference proceedings, Lake Buena Vista, Florida, 

June 5–9, 2006, p. 224: Available online at http://sofia.usgs.gov/geer/2006/index.html

Swain, E.D., Wolfert, M.A., Bales, J.D., and Goodwin, C.R., 2004, Two-dimensional hydrodynamic simulation of surface-water 

flow and transport to Florida Bay through the Southern Inland and Coastal Systems (SICS): U.S. Geological Survey WaterResources Investigations Report 03–4287, 56 p., plus 6 plates.

Swayze, L.J., and McPherson, B.F., 1977, The effect of the Faka Union Canal system on water levels in the Fakahatchee Strand, 

Collier County, Florida: U.S. Geological Survey Water-Resources Investigations Report 77–61, 19 p.

Thorne, D., Langevin, C.D., and Sukop, M.C., 2006, Addition of simultaneous heat and solute transport and variable fluid 

viscosity to SEAWAT: Computer and Geosciences v. 32, 1758–1768.

U.S. Army Corps of Engineers, 2004, Comprehensive Everglades Restoration Plan, Picayune Strand Restoration, Final 

Integrated Project Implementation Report and Environmental Impact Statement: Accessed on May 21, 2008, at 

http://www.evergladesplan.org/pm/projects/docs_30_sgge_pir_final.aspx

Wang, J.D., Swain, E.D., Wolfert, M.A., Langevin, C.D., James, D.E., and Telis, P.A., 2007, Applications of Flow and Transport 

in a Linked Overland/Aquifer Density Dependent System (FTLOADDS) to simulate flow, salinity, and surface-water stage in 

the southern Everglades, Florida: U.S. Geological Survey Scientific Investigations Report 2007–5010.