Three-dimensional Modeling of Tidal Hydrodynamics in the San Francisco Estuary

Simulations of circulation in the San Francisco Estuary were performed with the three-dimensional TRIM3D hydrodynamic model using a generic length scale turbulence closure. The model was calibrated to reproduce observed tidal elevations, tidal currents, and salinity observations in the San Francisco Estuary using data collected during 1996-1998, a period of high and variable freshwater flow. It was then validated for 1994-1995, with emphasis on spring of 1994, a period of intensive data collection in the northern estuary. The model predicts tidal elevations and tidal currents accurately, and realistically predicts salinity at both the seasonal and tidal time scales. The model represents salt intrusion into the estuary accurately, and therefore accurately represents the salt balance. The model’s accuracy is adequate for its intended purposes of predicting salinity, analyzing gravitational circulation, and driving a particle-tracking model. Two applications were used to demonstrate the utility of the model. We estimated the components of the longitudinal salt flux and examined their dependence on flow conditions, and compared predicted salt intrusion with estimates from two empirical models.


INTRODUCTION
Abundance or survival of several estuarine biological populations in the San Francisco Estuary is positively related to freshwater flow (Jassby and others 1995).Freshwater flow into the estuary in spring is regulated to control the position of 2 psu salinity at the bed, or X2 (Jassby and others 1995).This regulation is based on the observed relationships of abundance to flow, although some of these relationships have changed (Sommer and others 2007).The high cost of the water (Kimmerer 2002) has stimulated interest in investigating the mechanisms underlying the "fish-X2" relationships.
Of the many proposed mechanisms for the observed fish-X2 relationships (Kimmerer 2002), two are particularly appropriate to explore with a hydrodynamic model.One hypothesis proposes that transport processes may partially explain the observed fish-X2 relationships.In particular, for species that recruit from the ocean, increased residual circulation, including gravitational circulation, with seaward X2 may result in more rapid landward transport and retention of organisms in the low-salinity zone, possibly resulting in higher survival and subsequent abundance.Another hypothesis proposes that fish-X2 relationships can be explained by the covariability of habitat area with X2.The habitat area for some pelagic species can be described by a combination of commonly-measured variables such as depth, salinity, turbidity, and temperature (Feyrer and others 2007).
A hydrodynamic model used to test hypotheses regarding estuarine circulation must realistically represent estuarine circulation.This requires a threedimensional model that accurately predicts the horizontal and vertical distributions of salinity, particularly during periods of high freshwater flow when gravitational circulation is strong (Monismith and others 2002).A hydrodynamic model used to explore variability of habitat area should represent bathymetry and horizontal variability of salinity adequately.
The hydrodynamic model TRIM3D was applied to conduct the above investigations.This paper presents the calibration and validation of TRIM3D for the San Francisco Estuary as far landward as the western Sacramento-San Joaquin Delta.The model was calibrated to data from 1996-1998, and validated with data from 1994-1995.Because salinity is central to these fish-X2 hypotheses, we emphasize comparisons of predicted salinity to salinity observations.The numerical modeling study focuses on San Pablo Bay, Carquinez Strait, and Suisun Bay.The interest in this region is motivated by the presence of the 2 psu salinity contour in this region during typical flow conditions, the location of important habitat in this region for several fish species (Kimmerer 2004), and because transport processes in this region may be particularly important for the retention of larval fish in the San Francisco Estuary.
Two analyses were conducted to demonstrate the utility of the model.In the first, we used the model to determine salt flux under conditions of steady flow.Longitudinal dispersion coefficients were determined and compared with results from previous studies, with particular emphasis on the contribution of gravitational circulation to salt flux.We also compared the relationship of freshwater flow to X2 from steady state runs of the model with two alternative statistical models based on historical data.

METHODS
Several numerical hydrodynamic models have been applied to the San Francisco Estuary as part of research and environmental impact studies.One model that has been widely applied to the San Francisco Estuary studies is the Tidal Residual Intertidal Mudflat (TRIM) model (Casulli 1990;Cheng and others 1993;Casulli and Cattani 1994).Cheng and others (1993) performed depth-averaged simulations of hydrodynamics and achieved excellent calibration of tidal predictions against an extensive set of harmonic constants derived from observations of tidal elevation and tidal currents.Three-dimensional simulations of salinity in South San Francisco Bay were performed with the TRIM3D model (Gross and others 1999b) which achieved good comparison to harmonic constants and salinity data; it also showed substantial influence of baroclinic pressure gradients on salt transport even during conditions of minimal stratification.

Model Formulation
To derive the governing equations of the TRIM3D model, the Navier-Stokes equations are Reynoldsaveraged to remove fluctuations at the turbulent time scale (e.g., Hirsh 1988) and a standard approximation to simplify the influence of density differences, known as the Boussinesq approximation, is applied: (2) (3) sediment-water interface the bottom friction is specified by (7) where τ x and τ y are the bottom stress components in the x and y direction, respectively.
A quadratic stress formula is applied at the bottom boundary as follows: (8) where u b and v b are the horizontal velocity components in the lowest layer in the water column and C d is the coefficient of drag, computed as: (9) where z b is the height of the center of the bottom cell above the bed, z o is the roughness coefficient, and κ is the von Kármán constant, a dimensionless parameter typically given as 0.40 (e.g., Sturm 2001).This is the standard formulation in three-dimensional models for flows that are "hydraulically rough," meaning that the roughness elements are larger than the thickness of the viscous sublayer (e.g., Sturm 2001).At the free surface, the same quadratic stress formula is applied, but the coefficient of drag is specified as a function of wind speed using the formulation of Large and Pond (1981).Free slip conditions were applied at lateral boundaries.
where u, v, and w are the velocity components in the x, y, and z directions, respectively; t is time; ζ is the height above a reference elevation; f is the Coriolis parameter; g is gravitational acceleration; r o is the reference density and r' is the variation from reference density; and v h and v v are the coefficients of horizontal and vertical eddy viscosity, respectively (Casulli and Cattani 1994).Conservation of mass is equal to conservation of volume for incompressible fluids: (4) The free surface equation is obtained by integrating the continuity equation (Equation 4) over depth and using a kinematic condition (the condition that the vertical velocity at the free surface equals the time derivative of the free surface elevation) at the free surface, which yields (Casulli and Cattani 1994) (5) where H 0 is the depth measured downward relative to the reference elevation.Thus, the total water depth is given by H 0 + ζ.
The boundary conditions at the free surface are specified as a function of wind stresses as (Casulli and Cattani 1994) (6) where w x τ and w y τ are the wind stress components in the x and y direction, respectively.Similarly, at the The governing equation for salt transport is: (10) stratification is present, a mild restriction on the time step applies (Gross and others 1998).Advection of momentum is performed by an unconditionally stable Eulerian-Lagrangian Method (ELM) (Casulli and Cattani 1994).The numerical method conserves water volume, allows wetting and drying of tidal flats and marshes, and reduces to the two-dimensional TRIM2D method (Casulli 1990) if one vertical layer is specified.
The numerical method used for scalar transport includes a widely used flux-limiting method (Van Leer 1974) for horizontal advection and an unconditionally stable semi-implicit method for vertical advection (Gross and others 1998).The excellent local and global mass conservation properties of the method are documented by Gross and others (2002).
The numerical method used to solve the equations of the turbulence closure is a semi-implicit method that results in tridiagonal positive-definite matrices in each water column and ensures that the turbulent variables remain positive (Deleersnijder and others 1997).

MODEL INPUT AND BOUNDARY CONDITIONS
The model is driven by water surface elevation at the seaward boundary and freshwater inflows at several landward boundaries.Wind forcing, baroclinic forcing, precipitation, and evaporation are also represented.

Bathymetric Grid
The model domain includes the western Sacramento-San Joaquin River Delta, the remainder of the estuary seaward of the delta, and a portion of the coastal ocean extending to approximately 22 km west of the Golden Gate (Figure 1).The primary data source for the bathymetric grid of South San Francisco Bay, Central San Francisco Bay, and San Pablo Bay was National Oceanic and Atmospheric Administration (NOAA) Digital Elevation Model (DEM) data.The where s is the scalar concentration; ε h is the horizontal eddy diffusivity coefficient; and ε v is the vertical eddy diffusivity coefficient.A linear equation of state is used to relate salinity to density, allowing a substantial reduction in computational time relative to the use of a nonlinear relationship.Salinity was simulated but temperature was neglected because density gradients in most of the San Francisco Estuary depend predominantly on salinity gradients (Cloern and Nichols 1985).
The turbulence closure model used in this study includes a turbulent kinetic energy equation and a generic length-scale equation.The parameters of the generic length-scale (GLS) equation are chosen to yield the "gen" closure proposed by Umlauf and Burchard (2003) and the Kantha and Clayson (1994) quasi-equilibrium stability functions are used.This closure has several advantages over the commonly used Mellor-Yamada level 2.5 closure (Umlauf and Burchard 2003), and performs similarly to the GLS versions of the k-ε and k-ω closures (Warner and others 2005).

Numerical Method
TRIM3D is a semi-implicit finite-difference model which solves the three-dimensional Navier-Stokes equations on a staggered grid.More specifically, the staggered grid used is an Arakawa C grid (Arakawa and Lamb 1977), meaning that the free surface elevations and scalar concentrations are defined at cell centers and normal velocities are defined at cell sides.The central feature of the TRIM3D method (Casulli 1990) is a highly efficient solution method that is stable even at a large computational time step (Casulli and Cattani 1994).A preconditioned conjugate gradient method is used to solve the matrix of equations at each time step.For uniform density simulations this iterative matrix solver is guaranteed to converge without any restriction on the computational time step (Casulli and Cattani 1994).When The hydrodynamics of San Francisco Estuary are simulated on a 200-m horizontal grid with two major approximations (Figure 1).First, a portion of the coastal ocean bordering San Francisco was excluded to save computational time.Second, because representation of hydrodynamics and transport processes in the Sacramento-San Joaquin Delta would require roughly 25 meter resolution on a Cartesian grid (Monsen 2001), and representation of these processes in the delta was not of interest for the current purposes of the model, the majority of the delta is represented by a "false delta" consisting of two rectangles.These were sized to allow approximately correct tidal flow past Rio Vista and Jersey Point compared with ultrasonic velocity meter (UVM) flow data (Oltmann 1998) collected by the USGS.Such a false delta approach has been applied in several previous modeling efforts (e.g., Ganju and Schoellhamer 2007).The shape and bathymetry of the rectangles is somewhat arbitrary.The rectangle width and orientation was chosen to allow flow across the entire cross-section of the resolved portion of the Sacramento River and San Joaquin River while the length of each rectangle was adjusted to achieve a realistic tidal prism for the unresolved portions of the delta.The bathymetry of the rectangles slopes up linearly from the depths of the resolved channels to a maximum elevation of -1.2 meters NGVD.
The vertical grid spacing was 1 meter and the model bathymetric grid consists of 42,424 active water columns and 569,602 active grid cells.The time step chosen was 120 seconds, based on a Courant-Friedrichs-Lewy (CFL) condition on internal wave propagation (Gross and others 1998).

Boundary and Initial Conditions
The western boundary of the model domain is a long open boundary in the coastal ocean, located approximately 22 km seaward of the Golden Gate.Water level observations were available for the model calibration period at three stations near the model boundary: NOAA station 9415020 at Point Reyes, 49 km northwest of the Golden Gate, station 9414290 at Fort Point in San Francisco, and station 9413450 at Monterey, 143 km south of the Golden Gate.The observed tidal range (mean high water -mean low water) is 1.19 m at Point Reyes, 1.25 m at Fort Point, and 1.08 m at Monterey.The form number, the ratio of the sum of the amplitudes of diurnal tidal constituents to the sum of the semidiurnal tidal constituents, is typically approximated as (O 1 + K 1 )/ (M 2 + S 2 ) (Cheng and Gartner 1984).The observed form number is 0.83 at Point Reyes, 0.90 at Fort Point, and 0.96 at Monterey, based on harmonic constants reported by NOAA (2003).Simulations were performed to evaluate the suitability of observations from each station for forcing the model.The observations from each station were used to drive a simulation and the predicted form number was calculated at Fort Point and other locations for each simulation.Forcing the model with the Fort Point observations and Point Reyes observations yielded realistic predictions of form number inside the San Francisco Estuary while forcing the model with Monterey tides did not.Therefore either Fort Point or Point Reyes water level observations could be applied as a boundary condition with a simple amplification around mean sea level and phase lag applied for tuning, as opposed to forcing with harmonic constituents.Forcing with water level observations retains some non-tidal boundary forcing including effects of barometric pressure and coastal winds.
Tidal observations from both Point Reyes and Fort Point are available for the 1996-1998 calibration period.During the high delta outflow periods in the winter of 1997-1998, the observed tidally-averaged water level was significantly higher at Fort Point than at Point Reyes, presumably due to some combination of delta outflow and coastal winds.Therefore the Point Reyes observations were used to specify the boundary condition for the 1996-1998 simulation.No amplification or phase lags were applied in this simulation.The water level along the open boundary was assumed to be uniform.Observations from Fort Point were used for the 1994-1995 validation period, when observations were not available from Point Reyes.To account for the increase in average tidal range from the Point Reyes station to the Fort Point station, the observed water levels at Fort Point were multiplied by a factor of 0.9544.A phase lead of 30 minutes was applied to account for the phase difference between Fort Point and Point Reyes.
The salinity on the ocean boundary was specified as 33.5 psu.This boundary condition is approximately correct for most conditions (Dever and Lentz 1994) but salinity can be depressed substantially by winter storms (Wilkerson and others 2002).
Freshwater inflow from several rivers, creeks, and water pollution control plants (WPCPs) was represented in the simulations (Figure 1).Net outflows from the Sacramento River and San Joaquin River, computed by the "DAYFLOW" program (CDWR 1986) were specified daily.These delta outflows were inserted into the model domain at the landward (eastern) edges of the Sacramento River and San Joaquin River false delta rectangles (Figure 1).The salinity of all tributary inflows was assumed to be 0 psu.Napa River, Petaluma River, Alameda Creek, Guadalupe River, and Coyote Creek were represented explicitly in the simulations.Flow was estimated based on USGS and water district data and, in most cases, the gaged flows were adjusted to account for the ungaged flows that enter the streams downstream of the gaging stations.Daily-averaged flows from the San Jose/Santa Clara, Sunnyvale, and Palo Alto WPCPs were included in the simulations.
Wind forcing was applied at the water surface as a wind stress.Hourly wind speed and direction were obtained from the Bay Area Air Quality Control District from San Carlos, Point San Pablo, and Pittsburg (Figure 1).The San Carlos wind observations were used for South San Francisco Bay and a portion of Central San Francisco Bay.The Point San Pablo wind observations were used for the coastal ocean, Central San Francisco Bay, and San Pablo Bay.The Pittsburg wind observations were used for Carquinez Strait, Suisun Bay, and the western delta.Wind was assumed to be uniform within each of these regions.

Daily precipitation measured at the San Pablo
Reservoir in Berkeley was used to specify spatially uniform precipitation.Monthly evaporation data from Newark were used to specify spatially uniform evaporation.Newark was chosen because evaporation has a significant effect on salinity in South San Francisco Bay (Denton and Hunt 1986) but a smaller effect in other portions of the San Francisco Estuary.Pan evaporation was multiplied by a factor of 0.695 to convert pan evaporation to evaporation from a water body (Linsley and others 1982).
The initial water surface elevation was set at a constant and uniform 0 m NGVD.The initial velocities were specified as zero.The initial salinity field was specified based on the USGS transects.For the calibration simulation, salinity data were collected on December 17, 1996 (Baylosis and others 1998) from 6:22 am to 5:19 pm.These quasi-synoptic salinity data were interpolated linearly along the axis of the estuary and the initial salinity field was assumed to be laterally uniform.The initial salinity field was specified at 10:15 am on December 17, 1996, corresponding to the collection of a salinity profile in Central San Francisco Bay (USGS Station 18).
For the validation simulation, salinity data were set based on synoptic salinity data collected on March 16, 1994 using the same approach as for calibration.The observations were made from 6:33 am to 6:05 pm from South San Francisco Bay to the western delta.The initial salinity field was specified on March 16, 1994 at 11:20 am, corresponding to the collection of a salinity profile in Central San Francisco Bay.
Following Cheng and others (1993), the bottom roughness parameter z 0 was specified according to the bed elevation of the water column (Table 1).The z 0 values used in all simulations are realistic in shoals and intertidal regions and unrealistically small in channel regions.As noted by Gross and others (1999a), the low z 0 values in the channels are probably required to compensate for loss of momentum from the channel to the shoal resulting from numerical diffusion inherent in the Eulerian-Lagrangian Method (ELM) used for advection of momentum.

Model Parameters
Horizontal and vertical eddy viscosity and eddy diffusivity are often used as calibration parameters in estuary models.Horizontal eddy viscosity and diffusivity were set to zero in the simulations because realistic coefficients are the same order of magnitude as numerical diffusion associated with the TRIM3D model.The model results were insensitive to the eddy diffusivity parameter within the range from 0 to 10 m 2 s -1 .
Baseline vertical eddy viscosity and diffusivity coefficients are used in turbulence closure models so that stratification does not lead to unrealistically low vertical mixing.Some vertical mixing is always present due to breaking internal waves, shear instabilities, and other processes not accounted for by typical turbulence closures.A baseline vertical eddy viscosity and eddy diffusivity value of 10 -4 m 2 s -1 was applied, which is on the high end of typically applied baseline eddy viscosity and eddy diffusivity (Li and others 2005).A sensitivity test using baseline vertical turbulent eddy viscosity and diffusivity of 10 -6 m 2 s -1 suggested that predicted salinity had minimal sensitivity to this parameter within the range of 10 -6 m 2 s -1 to 10 -4 m 2 s -1 , except at extremely high delta outflows.

Observations Used in Calibration and Validation
The calibration simulation period of December 16, 1996 to April 1, 1998 was chosen because of the large range of freshwater flow during this period (Figure 2).An independent dataset for markedly different delta outflow was used for model validation.The validation period of April 1994 through March 1995 was chosen because of the availability of extensive field data from the Entrapment Zone Study in spring of 1994 in Carquinez Strait and Suisun Bay (Burau and others 1998;Kimmerer and others 1998).This data set includes high-frequency velocity data from acoustic doppler current profilers (ADCPs), salinity from continuous monitoring stations, and salinity transects (Figure 3).Eight salinity transects were collected during spring tides on April 27 and 28, and nine salinity transects were collected during neap tides on May 17 and 18. Water year 1994 was classified as a "critical" year, the lowest flow classification, on both the Sacramento River and San Joaquin River.
Salinity predictions during the validation period were also compared with additional USGS continuous monitoring data and synoptic salinity transect data (Burau and others 1998).
For both calibration and validation, the quality of fit between predicted and observed data was assessed by a cross-correlation procedure (RMA 2005).The mean values of predicted and observed continuous monitoring data over the calibration period were compared, and the cross-correlation analysis provided computed phase lag between predicted values and observed values.A linear regression of observed values and lagged predicted values was performed to calculate the amplitude ratio (the slope of the best linear fit), the offset (intercept) of the fit, and the coefficient of determination which represents the degree of scatter in the regression.Salinity transect data were compared with model output obtained from the grid cell containing each station for the time of the individual observation.
ADCP data were compared with model predictions using two approaches.First the predicted and observed depth-averaged current speeds were compared using the cross-correlation procedure described above.Second, velocity profiles during the peak of each flood and each ebb tide were averaged to calculate the average maximum flood and ebb tidal velocity profiles from observations and model predictions.
The USGS synoptic salinity observations were used to estimate X2 to compare with predicted X2 during the calibration and validation periods.Because X2 is the location where daily-averaged bottom salinity is 2 psu, and the USGS synoptic salinity observations do not directly indicate daily-averaged salinity, a method to estimate X2 values from the USGS observations was developed.First, the distance between the predicted 2 psu bottom salinity at the time of the USGS observations and the predicted daily-averaged 2 psu bottom salinity was calculated.Then this computed distance "shift" was applied to the synoptic observations to estimate observed X2.The comparison of observed and predicted X2 is particularly relevant to the intended applications of the model because the variability of physical processes and properties with X2 will be explored using the model.

Salt Flux Analysis
The "salt balance" equation (Fischer and others 1979) is a simplified but useful description of salt transport: (11 where Q is the tidally-averaged flow, S is tidally and cross-sectionally averaged salinity, K is the longitudinal dispersion coefficient, A is the cross-sectional area, and x is the longitudinal position.The salt balance equation applies to the longitudinal salinity distribution under tidally averaged steady state conditions.If these conditions are met, Equation 11 can be used to estimate longitudinal dispersion coefficients. Estimating the portion of the total dispersion coefficient associated with gravitational circulation and other individual processes requires detailed analysis of simulation results.The salt flux associated with individual physical processes can be estimated at any cross-section by an analysis method outlined in Fischer and others (1979).The longitudinal velocity where x is the longitudinal position of a crosssection, y and z are the lateral and vertical distances within a cross-section, and t is time.The velocity components are the cross-sectional and tidallyaveraged velocity (U a = Q/A), the deviation of the cross-sectional average from the cross-sectional and tidally-averaged velocity (U c ), the deviation of the tidally-averaged velocity from the cross-sectional and tidally-averaged velocity (u s ) and the remaining variability (u').The capital letters refer to depth-averaged quantities.The last two terms of Equation 12 are further decomposed into lateral and vertical variability ( 13) The same decomposition approach is followed for salinity.The cross-sectional area is decomposed into a tidal cycle average and variation from this average, !
The salt flux through a cross-section at any time is where g is gravity, r is density, x is longitudinal position, H is water column depth and u* is friction velocity.
To achieve steady state conditions necessary for an unambiguous result, six steady delta outflow scenarios were simulated with TRIM3D.All scenarios used an idealized tide based on the M2 and K1 harmonic where A is the cross-sectional area and the square brackets represent a cross-sectional average.The average salt flux during a tidal cycle is determined by averaging over the tidal cycle: (17) where the angle brackets represent a tidal cycle average.This notation follows Fischer and others (1979) closely except that square brackets are used instead of an overbar to represent cross-sectional averages.The decomposed velocity, salinity and area are substituted into Equation 17.Many product terms are zero or negligible (Dyer 1973).Keeping all terms that are expected to be significant in any part of the San Francisco Estuary results in the equation (18) constituents at the Point Reyes station (NOAA station 9415020).The M2 period was modified to 12.0 hours so that exactly two M2 cycles occur per K1 cycle (1 day).All other harmonic constituents and nontidal forcing at the seaward boundary were neglected, and oceanic salinity was held constant.The only difference among the six scenarios was the delta outflow which ranged from 55 m 3 s -1 to 2,810 m 3 s -1 .This range of flow extends from summer flows during a critically dry year to high delta outflow associated with peak winter flows.
Each simulation was run for several months until tidally-averaged salinity reached steady state.Q, S and A values for Equation 11 were calculated for each cross-section and each scenario analyzed.The salinity gradient (dS/dx) was estimated along the axis or main channel of the estuary.The local longitudinal salinity gradient was determined by a linear fit to the depth-averaged and tidally-averaged salinity over 12 km, which roughly corresponds to one tidal excursion.Dispersion coefficients were not calculated in regions with tidally-averaged salinity gradients of less than 0.05 psu km -1 or salinity less than 0.4 psu during any portion of the tidal cycle.
The fluxes at each cross-section were examined to evaluate the variability of the gravitational circulation dispersion component with Ri x .The root mean square (RMS) depth-averaged velocity along the axis of the estuary and an assumed coefficient of drag of 0.0025 were used in the estimate of friction veloc-ity.The tidally-averaged and depth-averaged salinity gradient and the channel depth along the centerline transect were also used to calculate Ri x .

Flow Dependency of X2
An X2 value was calculated for each steady delta outflow and these data were used to revisit regressions of X2 with delta outflow based on historical data (Jassby and others 1995; Monismith and others 2002).X2 was initially calculated from continuous salinity observations at several locations (Jassby and others 1995).That relationship was a time-series regression of X2 on the log of outflow, developed empirically from daily data and currently used to estimated X2 for regulatory and analytical purposes.This relationship is (20) where X2(t) is X2 at time t in days, Q is net delta outflow, and A, B, and C are coefficients (A = 10.2,B = 0.945, and C = 2.30 with Q in m 3 s -1 ).Monismith and others (2002) developed an alternative time-series model: (21) where α = 0.919, β = 13.57, and γ = -0.14 are the reported best-fit coefficients.Although the former model fit the data slightly better than the latter, only the power law function has a theoretical basis (Monismith and others 2002).

RESULTS
Results of the calibration and validation are summarized here, with example graphs indicating the com-parisons between model output and data.Complete sets of figures are presented in Appendix A for calibration and in Appendix B for validation.

Water Surface Elevation
Predicted water levels were compared to observations at five stations (Figure 3).Predicted mean water levels were very close to observed mean values, except at Antioch (Table 2).The amplitude ratio (Table 2) was close to 1 at Fort Point, Alameda, and Richmond, but substantial amplitude errors were present at Port Chicago and Antioch.The phase lags between predicted and observed water levels indicated accurate phase propagation in Central Bay and San Pablo Bay but substantial phase lags occurred in Suisun Bay and the western delta (Table 2).
The observed and predicted tidally-averaged water levels at Port Chicago were similar but the tidal amplitude was significantly under-predicted and some phase lag is noticeable (Figure 4).Similar figures for the other stations show excellent comparison of observed and predicted water levels in Central Bay and San Pablo Bay but limited correspondence at Antioch (Figures A4-A8, Appendix A) where water level is affected by reflection from the false delta regions.

Tidal Flows in the Western Delta
The geometry of the false delta regions was adjusted to tune the tidal prism to roughly match flow observations.A 29-day period with continuous UVM observations was used for calibration (Figure 5).

Synoptic Salinity Calibration
The predicted salinity along the axis of the San Francisco Estuary was compared with USGS synoptic sampling observations (Baylosis and others 1998) during all San Francisco Bay cruises between January 1, 1997 and April 1, 1998, except cruises that were limited to South San Francisco Bay.The average errors and standard errors are small relative to the large range of salinity conditions that occurred during the calibration (Table 4, Figures A12-A25, Appendix A).For example, on January 28, 1997, when the observed salinity field was strongly depressed as the result of the January storm, the predicted and observed salinity fields contained low salinity in Carquinez Strait and strong stratification though San Pablo Bay and Central Bay (Figure 6).During the following summer and fall, no substantial flow events occurred and both observed and predicted salinity increased slowly.
Predicted and observed salinity on November 6, 1997 (Figure 7) indicate substantial salt intrusion into the western delta and limited stratification.

Continuous Salinity Calibration
Predicted salinity was compared with salinity observed at continuous monitoring stations (Table 5, Figures A27-A40, Appendix A) and generally predicted observed trends.For example, the seasonal trends in salinity were predicted well at Martinez with salinity underestimated by an average of 1.4 psu at the bottom sensor and 0.5 psu at the top sensor (Figure 8 and Figure 9).The cross-correlation statistics (Table 5) suggest that salinity variability at the tidal time scale was predicted accurately at both sensors and that stratification is underestimated by the model on average at this location.The most significant errors in tidally-averaged salinity were underestimation of salinity at both sensors during spring and summer of 1997 (Figure 8 and Figure 9).

Tidal Velocity Validation
Depth-averaged current speeds were reasonably well predicted by the model at most stations (Table 6, Figures B1-B10, Appendix B).For example, the mean predicted speed at the Martinez ADCP was 0.65 m s -1 and the mean observed speed was 0.57 m s -1 , though the highest peak currents were typically underestimated by the model (Figure 10).The average peak ebb and flood velocity profiles at the Martinez ADCP (Figure 11) indicate that both the observed and predicted velocity profile during ebb was more sheared than the velocity profile on flood tides.

Synoptic Salinity Validation
Salinity was predicted accurately by the model on most dates (Table 7, Figures B11-B21, Appendix B).For example, the observed and predicted salinity profiles along the axis of the estuary on April 19, 1994 (Figure 12), slightly more than a month after initial conditions were applied, agreed well.Both observed and predicted salinity showed much greater penetration up the estuary during the validation period than in the calibration period.Note that small deviations between actual locations sampled by the USGS and the reported station locations caused the reported depth at each station to vary among profiles.
No substantial flow events occurred over the next several months, and both observed and predicted salinity profiles showed a gradual trend of increasing salt intrusion (Appendix B).On October 26, 1994, both observed and predicted salinity values were high and X2 was in the delta (Figure 13).The observed salinity for this cruise was underestimated by an average of 0.4 psu and the standard error was 0.87 psu (Table 7).
A flow event with peak outflow of over 7,000 m 3 s -1 occurred during January 1995, causing nearly complete flushing of salt from the delta and Suisun Bay by the time of the USGS cruise on January 18, 1995 (Figure 14).Both observed and predicted salinity were below 0.5 psu in the western delta and Suisun Bay and sharp horizontal salinity gradients had formed in Carquinez Strait and San Pablo Bay.All portions of the estuary with substantial salinity were strongly stratified in both observations and predictions.The salinity for this cruise was overestimated by an average of 0.45 psu and the standard error was 2.74 psu.Much of the error in salinity predictions occurred in South San Francisco Bay, most likely because inflows from several tributaries were not included in the model's boundary conditions.

Continuous Salinity Validation
During the validation period, salinity observations were available at 5 of the 8 stations used for the model calibration (Figure 3).The predicted salinities are typically similar to the observed salinity (Table 8, Figures B23-B30, Appendix B).The observed salinity at the bottom sensor at Point San Pablo was overestimated by an average of 0.8 psu (Figure 15) while the salinity at the top sensor was underestimated by an average of 1.0 psu (Figure 16).Therefore the stratification at this station was generally overestimated during this period.The amplification ratio of 0.72 for the bottom sensor indicates that the range of salinity was underestimated by the model, while the amplification ratio of 0.92 at the top sensor (Table 8) suggests that the salinity range at that location was predicted accurately.

Entrapment Zone Synoptic Salinity
The eight salinity transects collected on April 27 and April 28, 1994, show the evolution of the salinity field in Suisun Bay during 2 tidal cycles of a spring tide.Approximately two hours were required to complete sampling of each salinity transect.In the following discussion, the time associated with each transect is the time at which the first salinity profile of the transect was collected.
Observed and predicted salinity transects are similar, although the model consistently underestimated salinity during these spring tide conditions (Table 9 Figures B31-B47, Appendix B).Because of the strong longitudinal salinity gradients in Suisun Bay, these salinity errors corresponded to relatively small errors in isohaline location.The salinity transect collected on April 27 at 9:44 am (Figure 17) was collected early in a flood tide with vertically well-mixed conditions which were predicted accurately by the model though the predicted salinity was underestimated by an average of 1.48 psu.The predicted isohalines were typically located 3 to 6 km seaward of the observed isohalines for this transect.The fourth transect, collected during an ebb tide, showed increased longitudinal salinity gradients at the seaward end of Suisun Bay and some stratification in Suisun Bay (Figure 18).The average underestimate of salinity was 0.64 psu and the standard error was 0.63 psu.
Salinity transects collected during neap tides on May 17 and May 18 indicate somewhat greater stratification in both observed and predicted profiles than was present during the spring tide period.The first transect during this period (Figure 19), was collected early during an ebb tide at 6:46 am, with salinity ranging from below 0.5 psu to over 18 psu and stratification ranging from 2 psu to 8 psu in the seaward portion of Suisun Bay.The predicted salinity shows a similar distribution, though salinity is underestimated by an average of 0.7 psu and stratification is also under-predicted.In addition, the predicted 2 psu isohaline is approximately 5 km too far seaward.The third transect (Figure 20), was collected near low water in fairly well-mixed conditions with salinity ranging longitudinally from below 0.5 psu to over 8 psu.The predicted stratification was also minimal and salinity was underestimated by an average of 0.76 psu.
The average error for each transect during the validation period was negative, indicating underestimation of salinity (Table 9).The mean error in predicted salinity of all the spring tide transects was -1.5 psu and the mean error of the neap tide transects was -1.0 psu.The salt flux analysis results in one value of calculated K gc and Ri x at each cross-section and flow scenario which can be used to develop a relationship between K gc and Ri x .The fitting equation K gc /u*H = a Ri x 2 used by Monismith and others (2002) was generalized to the form ( 22) and the coefficients a and b were determined by a best fit approach at each cross-section.The relationship between K gc and Ri x varied substantially among cross-sections (Table 10).An example of the variability of K gc with Ri x (Figure 25) is provided for a cross-section located 1 km west of the entrance of Carquinez Strait.The curve fit equation describes the variability in dispersion coefficients accurately (r 2 = 0.998) and the exponent (1.9) calculated at this location was similar to the exponent found by Monismith and others (2002).

Salt Flux Analysis
Dispersion coefficients were calculated at 28 different cross-sections in the San Francisco Estuary (Figure 21) for 6 different flow scenarios.The salinity field along the axis of the estuary is shown for 4 of these scenarios in Figure 24.The portion of dispersion coefficients associated with gravitational circulation (K gc ) drops substantially in Suisun Bay (Figure 23 and Figure 24) relative to coefficients in Central Bay and San Pablo Bay.These results are consistent with field observations which suggest that the reduced depth at and east of the Benicia shoal (km 55) reduced the strength of gravitational circulation (Burau and others 1998;Schoellhamer 2001).Note that K gc increases with flow in Central San Francisco Bay and San Pablo Bay but decreases with increased flow in the western delta.This occurs because the longitudinal salinity gradients increase with flow in Central San Francisco Bay and San Pablo Bay but decrease with increased flow in the western delta and a portion of Suisun Bay as the increased flow flushes salt from those regions (Figure 22).

Flow Dependency of X2
Observed X2, calculated from the USGS synoptic salinity observations during the calibration period, ranged from 40 km to 92 km and X2 was predicted accurately by TRIM3D through this large range of observed X2 (Figure 26).The average error in X2 predictions was -1.5 km and the standard error was 3.3 km.Observed X2 during the validation period ranged from 53 km to 96 km (Figure 27).The average error of X2 predictions was -0.6 km and the standard error was 3.3 km.A linear regression of the predicted X2 and observed X2 shows limited scatter and virtually no bias in the TRIM3D predictions (Figure 28).In contrast, the X2 values predicted by Equation 21show more scatter around the best-fit line and a bias toward overestimating X2 at low observed X2 values and underestimating X2 at high X2 values.
The TRIM3D predictions were used to develop a bestfit power law function of the form shown in Equation 21 (Figure 29).A best-fit power law function was also developed from interpolated X2 values calculated using continuous salinity observations.The interpolated X2 dataset was used previously by Monismith and others (2002).In our analysis, the interpolated X2 values were used for the simulation periods only (Figure 29), so that the regression based on interpolated X2 values could be directly compared to the regression based on TRIM3D results.Interpolated X2 values less than 55 km were not included because these were calculated by extrapolation beyond the westernmost continuous monitoring station at Martinez, and are therefore less reliable than the interpolated values.The coefficients of the regression of interpolated X2 on flow are somewhat different from the coefficients reported by Monismith and others (2002), primarily because their observation record was 1967-1990.Lastly, a steady state power function was developed using the results of the steady flow TRIM3D simulations (Figure 29).The coefficients (Table 11) are different for the different X2 datasets.Some differences between the steady and unsteady simulations are likely to arise because unsteadiness in X2 cannot be represented accurately by Equation 21, which has a single time constant of variability, 1/(1-α) (Monismith and others 2002).Table 11 Results of analysis of X2 with flow.The α, β and γ values are the curve fit parameters for the equation .For steady state conditions, the equation reduces to .The " Monismith and others (2002)" row shows parameters reported by Monismith and others (2002).The "Predicted X2" row shows regression results using X2 predicted by TRIM3D during the two simulation periods.The "Interpolated X2 (> 55 km)" row shows regression results for X2 values calculated by interpolation of continuous salinity observations using only the observations collected during the TRIM3D simulation periods.The "Predicted Steady X2" row shows regression results from steady flow scenarios.The first hypothesized mechanism for the observed fish-X2 relationships to be explored using this model involves the covariability of habitat area with X2 (Kimmerer and others 2009).The analysis presented here indicates that water levels are predicted accurately through the tidal cycle, the spring-neap cycle, and due to episodic events such as large storms (Figure 4).Accuracy decreases in and near the delta (Table 2) due to artifacts that arise from not resolving most of the delta geometry.

Regression
Salinity is generally predicted well for a large range of delta outflows, in particular the extreme flows of January 1997 and winter 1997-1998 (Figure 2).Uncertainties in the habitat use by most species are likely to be much larger than uncertainties in modelbased predictions of salinity or water level.Therefore, to the extent that habitat area is defined by a combination of salinity and water depth, the present model is appropriate for analysis of variability of habitat area with X2.However, it should also be noted that the salinity calibration and validation data were collected predominantly in channels.The accuracy of salinity predictions in shoal and intertidal regions is less well established than the accuracy of salinity predictions in the main channels of the estuary.
A time constant for variability of X2 with flow can be determined for each steady delta outflow simulation by fitting the convergence to a steady state X2 with the equation ( 23) where X2 steady is the steady state X2, and A and B are fitting parameters.The time constant for convergence is then B -1 .This time constant decreased as the specified steady flow increased among scenarios with a time constant of 18 days for delta outflow of 260 m 3 s -1 and a time constant of roughly 1.4 days for delta outflow of 2810 m 3 s -1 .Reliable time constants for the scenarios with flow less than 260 m 3 s -1 could not be estimated, due to influence of the false delta regions, but should presumably be longer than 18 days.

DISCUSSION
The San Francisco Estuary TRIM3D model was developed to predict the spatial distribution of salinity, analyze gravitational circulation, and drive a particle-tracking model.The model calibration and validation indicates that the model predicts hydrodynamics and salinity fairly accurately over a large range of freshwater flow, suggesting that transport processes are represented adequately.The deviations between the predictions and observations could have The second hypothesized mechanism for the observed fish-X2 relationships being explored with the TRIM3D model involves the effect of transport processes on recruitment and retention of some organisms into the estuary.The velocities, water levels, salinity, and turbulence properties output from the model are currently used to drive a particle tracking model to investigate the recruitment of species from the ocean into the estuary for different outflow conditions and vertical migration behavior.
The results presented indicate that the model represents transport processes adequately and, in particular, accurately predicts the variability of X2 with delta outflow (Figure 26 and Figure 27).Gravitational circulation is an important transport process in much of the estuary, and the particular form and slope of the X2 relationship to outflow are tied theoretically to the strength of gravitational circulation and its response to compression of the salinity field and depth of the channels (Monismith and others 2002).Gravitational circulation is of particular interest for the fish-X2 studies because it is hypothesized that recruitment of bay shrimp or starry flounder may increase with increasing gravitational circulation.The strength of gravitational circulation varies with outflow by orders of magnitude (Monismith and others 2002), so the fairly accurate prediction of salinity over a large range of outflow is compelling evidence that gravitational circulation is represented adequately in the model.The accuracy of predictions of spring-neap, seasonal, and spatial variability in stratification and tidally asymmetrical flow also suggest that the model is representing gravitational circulation accurately, although stratification is often under-predicted in Suisun Bay.Again, the uncertainties associated with ecology, such as fish behavior, are likely to be much larger than uncertainties related to representation of transport processes.
Though the present model performs well and is an appropriate tool for the intended applications, substantial uncertainty is associated with the model predictions.Much of this uncertainty is associated with boundary conditions.The specified flows are particularly uncertain during peak flow periods when reported outflow is outside the range for which rating curves of flow measurement have been developed and tested.The flows are also quite uncertain during summer and fall when consumptive use within the delta, estimated on a monthly time step for use in the DAYFLOW program, can be of the same order of magnitude as delta tributary flows.The exclusion of many small tributaries in South Bay is a likely source of significant error in salinity predictions during peak flow periods.The assumed salinity on the model boundary (33.5 psu) may be a substantial overestimate during and following large outflow events when observed salinities in the coastal ocean can be depressed substantially (Wilkerson and others 2002).The current representation of bottom friction is also approximate and does not explicitly represent the effect of wind waves on bottom friction which is known to be substantial in shallow regions of the San Francisco Estuary (Bricker and others 2004).
The uncertainty associated with resolution of the geometry and bathymetry of the estuary also has a substantial effect on hydrodynamic and salinity predictions.The most obvious approximation is the representation of the delta with the false delta geometries.This approximation results in errors in water levels and tidal flows and is likely to affect the timing with which delta outflows arrive in the western delta.The use of the false delta geometries calls into question the validity of predicted hydrodynamics and transport within one tidal excursion, roughly 5-10 km, of these regions.Recent unstructured grid modeling of the San Francisco Estuary on a grid extending from the coastal Pacific Ocean through the entire Estuary, including all of the delta, predicted tidal elevations accurately through Suisun Bay and in the delta (MacWilliams and others 2007; MacWilliams and others 2008).In addition to the obvious geometric approximations at the eastern boundary of the model domain, predictions are also affected, particularly locally, by limited resolution of small-scale bathymetric features, some of which are located near monitoring stations.Given the substantial uncertainty associated with boundary conditions and limited model resolution, the relatively small errors in the salinity predictions are encouraging.

Transport Analyses
Previous modeling and analytical efforts assumed uniform bathymetry in both the longitudinal and lateral directions to derive the relationship K gc /u*H = a Ri x 2 (Monismith and others 2002).In the San Francisco Estuary, depth, velocity and salinity vary both in the lateral and longitudinal directions.For this reason, it is not clear if the thalweg depth is always appropriate to use in calculating Ri x .An additional complication in this analysis is longitudinal variability of bathymetry.Some distance is required to allow development of strong gravitational circulation, so sills along the axis of the estuary can limit the development of gravitational circulation.Because longitudinal and lateral bathymetric variability are substantial in the San Francisco Estuary, the relationship used by Monismith and others (2002) has been generalized Equation 22 with all variables defined at the location of the channel thalweg.This generalized equation is somewhat arbitrary; however, the high r 2 values at most sections in 10 suggest that this equation generally captures much of the variability of K gc with Ri x .The exponent b is close to 2.0 at many cross-sections, but both the exponent and proportionality constant vary substantially, presumably as a result of the complex bathymetry of the San Francisco Estuary.
The relationship between flow and predicted X2 is significantly different from the Monismith and others (2002) relationship determined from empirical data, (Table 11).The reported X2 values were determined by interpolation from a limited set of continuous monitoring stations.The most seaward continuous monitoring station used in the analysis by Monismith and others (2002) was Martinez (km 55), so values west of that station were undetermined and had to be extrapolated from flow.This undoubtedly introduced error in the Monismith and others (2002) regression since the water depth changes abruptly at that point and stratification is much more common in and west of Carquinez Strait than to the east.The use of a false delta limits the extent of valid X2 prediction to roughly one tidal excursion seaward of the false delta, which may a substantial source of error in the regressions based on predicted X2.
The steady-flow TRIM3D results yielded a different exponent to the relationship (Table 11) than the unsteady flow scenarios, suggesting that unsteadiness is not well represented by Equation 21.The adjustment time for salinity in the San Francisco Estuary should increase with decreased flow (MacCready 1999).An analysis of the response of predicted X2 to changes in flow shows large variability in the time constant.In contrast, a single time constant is used in Equation 21.An additional issue suggested by the steady flow X2 predictions is that the exponent in Equation 21 could change as X2 moves seaward, becoming smaller in magnitude (less negative) as gravitational circulation increases in importance (Figure 24), increasing resistance to seaward movement of the salinity field with further increases in flow.These and other possibilities could be investigated through additional modeling with a wider range of flow scenarios.

CONCLUSIONS
The results presented demonstrate that the TRIM3D model for the San Francisco Estuary is an appropriate tool for the intended applications in the study of hydrodynamic influences on habitat and movements of organisms.The model has been applied to estimate habitat area based on water depth and salinity for a range of freshwater flows (Kimmerer and others 2009).In addition, the velocities, water levels, salinity, and turbulence properties output from the model are currently used to drive a particle tracking model to investigate the recruitment of species from the ocean into the estuary for different outflow conditions and vertical migration behavior.
An analysis of salt flux indicates large variability in the magnitude of transport mechanisms both spatially and with delta outflow.Seaward of Martinez, gravitational circulation is the dominant transport mechanism, while several transport mechanisms are important landward of Martinez.
An analysis of X2 variability with flow suggests that improvements could be made to previously developed regressions by accounting for variability in the response time of X2 to changes in flow, allowing response time to decrease as flow increases.In addition, further analysis could allow representation of the spatial variability in the response of X2 to changes in steady flow, presumably with a weaker response in the deeper seaward portion of the estuary where gravitational circulation is a dominant transport mechanism.

Figure 1
Figure 1 Model domain, bathymetry, and locations of freshwater input with distance from the Golden Gate along the northern axis of the estuary (yellow line) Both 1997 and1998 were classified as wet years by the Interagency Ecological Program (IEP); the largest daily delta outflow on record occurred on January 3, 1997, and water year 1998 had the second highest annual mean flow of any water year since 1955.Two weeks of model spin-up were allowed before model predictions were compared to observations, starting on January 1, 1997.The calibration dataset included water level observations collected by NOAA, UVM flow observations collected by the USGS, salinity data from continuous monitoring sites operated by the USGS, United States Bureau of Reclamation (USBR), and the California Department of Water Resources (CDWR), and synoptic salinity observations by the USGS, consisting of vertical profiles of salinity at 1 meter vertical resolution at 38 sampling locations along the axis of the San Francisco Estuary (Baylosis and others 1998; Figure3).Near-surface and near-bottom salinity observations were available from each continuous monitoring station except Fort Point and Martinez, for which only near-surface salinity data were available.Substantial data gaps are present in the continuous observations due to instrument failures and the removal of suspect data by the agencies.

Figure 2 Figure 3
Figure 2 Net delta outflow during the calibration and validation simulation periods (16)   Many terms in this equation are associated with one or more physical processes.Some particularly important terms for the analysis of transport in the San Francisco Estuary are:• -advective salt flux (QS in the salt balance equation)•-gravitational circulation • -the primary term associated with unsteady vertical shear•-the primary term associated with tidal dispersion Previous studies have indicated that the strength of gravitational circulation increases with the horizontal Richardson number, Ri x(Monismith and others 1996,  Stacey and others 2001):   (19)

Figure 5
Figure 5 Observed and predicted tidal flows in the San Joaquin River at Jersey Point and Dutch Slough.Top panel, tidal variability of flows during two spring-neap cycles; lower left panel, tidally-averaged variability of flows; lower right panel, predicted and observed flows during simulation period with cross-correlation statistics.

Figure 7 Figure 6
Figure 7Observed and predicted salinity profiles at synoptic sampling stations, interpolated along the axis of the San Francisco Estuary on November 6, 1997

Figure 9
Figure 9 Observed and predicted salinity at the Martinez top sensor during the calibration period.See caption for Figure 8.

Figure 11
Figure 11 Observed and predicted average peak flood and average peak ebb velocity profiles at the Martinez ADCP.

Figure 12 Figure 13
Figure 12 Observed and predicted salinity profiles at synoptic sampling stations, interpolated along the axis of the San Francisco Estuary on April 19, 1994

Figure
Figure Observed and predicted salinity profiles at synoptic sampling stations, interpolated along the axis of the San Francisco Estuary on January 18, 1995

Figure 16
Figure 16 Observed and predicted salinity at the Point San Pablo top sensor during the validation period.See caption for Figure 15.

Figure 17 Figure 18 Figure 19
Figure 17 Observed and predicted salinity profiles at synoptic sampling stations for the Entrapment Zone Study, interpolated through Suisun Bay and the western Delta on April 27, 1994 at 9:44 am, early in a flood tide

Figure 20
Figure 20 Observed and predicted salinity profiles at synoptic sampling stations for the Entrapment Zone Study, interpolated through Suisun Bay and the western delta on May 17, 1994 at 12:30 pm, near low water

Figure 21 Figure 23
Figure 21 Cross-sections used in salt flux analysis

Figure 24
Figure 24Estimated dispersion coefficients associated with gravitational circulation (K gc ) for all flow scenarios.The horizontal scale is distance along the estuary from the Golden Gate.

Figure 26 Figure 28 Figure 29
Figure 26 Observed X2, predicted X2, X2 estimated from the Monismith and others (2002) regression relationship and net delta outflow during the calibration period

Table 2
Observed and predicted water levels and cross-correlation statistics for water level monitoring stations

Table 4
Average error and standard error for each synoptic sampling cruise covering the axis of the San Francisco Estuary

Table 3
Predicted and observed mean flow and cross-correlation statistics for USGS flow monitoring stations

Table 5
Predicted and observed mean salinity and cross-correlation statistics for continuous salinity monitoring stations

Table 7
Average error and standard error for each synoptic sampling cruise covering the axis of the San Francisco Estuary

Table 6
Predicted and observed mean depth-averaged speed and cross-correlation statistics for depth-averaged speed at ADCP stations

Table 8
Predicted and observed mean salinity and cross-correlation statistics for continuous salinity monitoring stations

Table 9
Average error and standard error for each entrapment zone study synoptic salinity sampling cruise

Table 10
Results of analysis of K gc with Ri x .The a and b values are the curve fit parameters for the equation Figure 25 Estimated dimensionless dispersion coefficients for gravitational circulation at section 9 as a function of horizontal Richardson number with corresponding best fit equation and r 2 value.