**Towards the Development of the National Ocean Service San Francisco Bay Operational Forecast System**

#### **Machuan Peng, Richard A. Schmalz Jr., Aijun Zhang and Frank Aikman III**

**Abstract:** The National Ocean Service (NOS), Center for Operational Products and Services installed a Physical Oceanographic Real Time System (PORTS) in San Francisco Bay during 1998 to provide water surface elevation, currents at PORTS prediction depth as well as near-surface temperature and salinity. To complement the PORTS, a new nowcast/forecast system (consistent with NOS procedures) has been constructed. This new nowcast/forecast system is based on the Finite Volume Coastal Ocean Model (FVCOM) using a computational domain, which extends from Rio Vista on the Sacramento River and Antioch on the San Joaquin River through Suisun and San Pablo Bays and Upper and Lower San Francisco Bay out onto the continental shelf. This paper presents the FVCOM setup, testing, and validation for tidal and hindcast scenarios. In addition, the San Francisco Bay Operational Forecast System (SFBOFS) setup within the NOS Coastal Ocean Model Framework (COMF) is discussed. The SFBOFS performance during a semi-operational nowcast/forecast test period is presented and the production webpage is also briefly introduced. FVCOM, the core of SFBOFS, has been found to run robustly during the test period. Amplitudes and epochs of the M2 S2, N2, K2, K1, O1, P1, and Q1 constituents from the model tide-only simulation scenario are very close to the observed values at all stations. NOS skill assessment and RMS errors of all variables indicate that most statistical parameters pass the assessment criteria, and the model predictions are in agreement with measurements for both hindcast and semi-operational nowcast/forecast scenarios.

Reprinted from *J. Mar. Sci. Eng.* Cite as: Peng, M.; Schmalz Jr., R.A.; Zhang, A.; Aikman III, F. Towards the Development of the National Ocean Service San Francisco Bay Operational Forecast System. *J. Mar. Sci. Eng.* **2014**, *2*, 247-286.

#### **1. Introduction**

In 1998, the National Oceanic and Atmospheric Administration (NOAA)/National Ocean Service (NOS) installed a Physical Oceanographic Real Time System (PORTS) in San Francisco Bay to provide water surface elevation, current, near-surface and near-bottom temperature and salinity, and meteorological information to promote safe and efficient navigation in this area [1]. As PORTS only supplies measured data at selected stations, NOS' National Operational Coastal Modeling Program (NOCMP) is developing an Operational Forecast System (OFS) to complement the service.

Many researchers have employed different versions of TRIM (Tidal, Residual, Intertidal, Mudflat Model) [2–4] for their purposes in the bay. For example, Cheng and Smith [5] employed the TRIM2D (two-dimensional TRIM) for the San Francisco Bay Marine Nowcast. The TRIM3D was recently applied by Gross *et al.* [6] to the entire San Francisco Bay. The UnTRIM, an unstructured version of TRIM3D, has also been applied to San Francisco Bay by MacWilliams and Cheng [7].

Two- and three-dimensional models have been extensively applied to investigate the hydrodynamic and morphologic processes in San Francisco Bay. Barnard *et al.* [8–10] report the existence of sand waves with heights in the order of 2 meters at the entrance of the Bay and consider coastal process evolution and the numerical prediction of severe storms on the coastline initially using the two-dimensional vertically integrated mode of the Delft3D-FLOW model [11]. Uslu *et al.* [12] developed a very high resolution two-dimensional vertically integrated model for tsunami forecasts in this region.

Fringer *et al.* [13] developed the non-hydrostatic option SUNTANS (Stanford Unstructured Nonhydrostatic Terrain-following Adaptive Navier-Stokes) model, which has also been applied in San Francisco Bay by Chua and Fringer [14].

Some of these structured and unstructured models [15,16] have been employed to understand the role of stratification and baroclinic circulation on salt intrusion in the northeastern part of Figure 1. They focus on the dynamic interactions between fresh water from the Sacramento and San Joaquin Rivers and salt water from the open ocean. Their studies, as will be seen later in this paper, have great value in evaluating the advantage and disadvantage of using a flow or stage river boundary condition for these two rivers.

**Figure 1.** The bathymetry of San Francisco Bay Operational Forecast System and the major gauge stations. The definition of "delta" in this paper is the region with complicated water channels to the east of Antioch and Rio Vista. The locations of the three major bays are indicated. Note: the domain in this figure is a little larger than the model grid domain as shown in Figure 2.

**Figure 2.** The SFBOFS grid structure and the open boundaries. Measured river flow data from USGS are used as river forcings for the five small rivers in blue. Measured river flow or river stage data at Rio Vista and Antioch may be used as river forcings, respectively, for Sacramento River and San Joaquin River. Technically, if river stage data are employed, the grid points across these two major rivers are treated the same way as those on the Pacific Open Boundary.

The primary objective of the NOCMP, however, is to develop and operate a national network of OFSs to support NOAA's mission goals and priorities. This ongoing San Francisco Bay Operational Forecast System (SFBOFS) will become a new member of the existing OFS family. Up to now, NOCMP has successfully developed CBOFS (Chesapeake Bay), DBOFS (Delaware Bay), TBOFS (Tampa Bay), NGOFS (Northern Gulf of Mexico), CREOFS (Columbia River Estuary) and other OFSs which along the Atlantic, Gulf of Mexico, Great Lake and Pacific Coasts.

With the use of the COMF (Coastal Ocean Model Framework) on NOAA's High Performance Computer (HPC), each OFS automatically integrates NOAA's observing system's data streams and the forecast output from meteorological and basin scale ocean models to generate necessary model input forcings, and then perform hydrodynamic model predictions with such forcings. Also with COMF, these OFSs perform nowcast and short-term forecast predictions (48 hours in most case) of pertinent parameters which include water levels, currents, salinity, and temperature and disseminate them to users. A state-of-the-art numerical hydrodynamic model driven by real-time data and meteorological, oceanographic, and river flow (or stage) forecasts forms the core of the end-to-end system. For detailed information on the COMF refer to Zhang *et al.* [17]. NOS CO-OPS is evolving to support two hydrodynamic models: ROMS for structured grid applications and FVCOM for unstructured grid applications.

As San Francisco Bay (Figure 1) has complex topography and shallow water features (the average water depth is less than 5 meters in the Bay), the well tested unstructured Finite Volume Coastal Ocean Model (FVCOM) [18–20] is employed as the core of the SFBOFS. Another reason to employ FVCOM lies in the fact that this model has already been ingested into the COMF-HPC as one of the major core models.

This paper shows the major steps in how SFBOFS has been developed, assessed and put into quasi-operational status. First, FVCOM is briefly reviewed in Section 2, followed by an overview of the model's setup in Section 3 for the tide and hindcast cases. Section 4 presents the model's astronomical tide-only scenario simulation evaluation, while the hindcast skill assessments are described in Section 5. The COMF setup and assessment of the quasi-operational nowcast/forecast test are discussed in Section 6. Conclusions and discussion are given in Section 7.

#### **2. The Model Overview, its Grid, and Subsequent Revisions**

The physics of the FVCOM model and many aspects of the computational scheme are equivalent to the widely used Princeton Ocean Model. The FVCOM model solves the three-dimensional, vertically hydrostatic, free surface, turbulent averaged equations of motions for a variable density fluid. The model uses a triangular unstructured horizontal grid with a generalized sigma vertical coordinate. Dynamically coupled transport equations for turbulent kinetic energy, turbulent length scale, salinity and temperature are also solved. The two turbulence parameter transport equations implement the Mellor-Yamada level 2.5 turbulence closure scheme. The numerical scheme employed in FVCOM to solve the equations of motion is summarized in [18–20].

The FVCOM application to San Francisco Bay uses external forcing of water level, ocean density, wind, sea level atmospheric pressure, air temperature, relative humidity, downward long wave radiation, short wave radiation, and fresh water discharges entering the model domain. The model calculates water levels, three dimensional velocity, salinity, and temperature.

The horizontal grid structure and open boundaries are shown in Figure 2. This grid was developed using the Surface Water Modeling System (SMS) Version 10.1 as described by Brigham Young University Surface Modeling Laboratory and was based on the VDATUM grid developed for the coastal waters of North/Central California, Oregon and Western Washington [21]. The open boundary of the San Francisco Bay grid was developed from this grid in the near shelf region external to the Bay. It was necessary to modify the VDATUM grid such that the outer boundary of the San Francisco Bay grid follows an approximate circular arc with one of the element sides near orthogonal to the boundary arc. The grid contains 102264 elements and 54120 nodes with a minimum depth of 0.2 m and maximum depth of 106.8 m [22]. A uniform 20-layer sigma level vertical discretization was considered.

The following element quality checks were used: (1) minimum and maximum interior angles of 10 and 130 degrees, respectively, (2) maximum slope of 0.1, (3) maximum adjacent element area change ratio of 0.5, and (4) maximum number of elements connected to a node of 8. Note the slope corresponds to the maximum allowed gradient of the edge length inside the domain. The slope determines how fast the mesh size will increase toward the middle of the region. A small slope order of 0.1 means small meshes. The paving method was used, which uses an advancing front technique to fill the polygon with elements. Based on the vertex distribution on the boundaries, equilateral triangles were created on the interior to define a smaller interior polygon. Overlapping regions were removed and the process is repeated until the region is filled. Interior nodal locations are relaxed to create better quality elements. Several triangles were adjusted such that the minimum interior angle was at least 30 degrees to improve FVCOM stability. In addition, along the open boundaries, the element topology was adjusted such that each boundary element contained only one boundary side. The triangle lengths are sufficiently small, that a reasonable M2 wavelength to grid size is obtained as shown in Figure 3, where element lengths decrease from 400 to 1700 m along the open ocean boundary to a near uniform resolution of order 150 m throughout the interior bays and into the lower delta.

Several modifications were made in the development of SFBOFS to Version FVCOM 3.1.6. It should be noted that if the HEATING\_CALCULATED\_ON options is selected then the AIR\_PRESSURE\_ON option must be selected. While the sea level atmospheric pressure field is needed for the heating calculations, its gradient does not need to be applied in the momentum equations. In fact, for tidal simulations this is not correct. For tidal simulations with the heat flux calculations selected, it is necessary to provide a constant sea level atmospheric pressure field (1013 mb). Also, if one selects AIRPRESSURE\_ON = F in namelist, the flag FLAG\_28 = íDAIR\_PRESSURE in file make.inc should be noted.

The bottom roughness fix reported by Warner [23] for wetting/drying was added in file brough.F. In model testing, with the min\_depth as 0.05 m, the model ran successfully and works for the wetting/drying case in San Francisco Bay. A Newtonian damping sponge layer was implemented by Lettmann [24], which provides a more robust implementation of the clamped water level open ocean boundary condition.

In the shallow mud flat regions of the Bay, there was also an issue with overheating. As a result, subroutine vdif\_ts.F was modified to limit the short wave radiation and total heat flux as a function of depth. For depths less than 10 m, the fluxes were set to zero. In this manner, the heat transfer is due to only advection and diffusion. There, the zeta1\_eff and zeta2\_eff parameters which control the attenuation of the short wave radiation are set never to be less than 30% of the water depth and therefore always allow attenuation. In total, the following routines are involved in the above modifications:


The interaction between the hydrodynamic and the sediment-water interface, particularly in the shallow water mudflat areas, which occupy some 16% of the Bay surface area, is an area where further research is needed. Fang and Stefan [25] considered the dynamics of heat exchange between the sediment and the bottom boundary layer for several hypothetical lakes. They found that the direction of the heat transfer reverses frequently on daily timescales as well as following an overall seasonal cycle based on weather conditions at Minneapolis-St. Paul, MN. Smith [26] performed a series of heat budget studies in Indian River Lagoon, FL, to estimate the water-sediment heat exchanges using assumed values for conductivity and density. The study sought to characterize subseasonal heat fluxes and temperature changes in the sediment and overlying estuarine waters.

**Figure 3.** The SFBOFS grid resolution structures are depicted. The element length sizes range from order 400–1700 m along the ocean boundary to order 200 m at the Bay Entrance as shown in the top panel. Within the Bay region the element lengths are reasonably uniform of order 150 m as shown in the lower panel with finer resolution around a few small islands of order 50 m.

The bottom stress formulation in shallow water for wetting and drying has received continuing interest. Research by Xue and Due [27], Uchiyama [28], Oey [29,30], and Oey *et al.* [31] has indicated that the bottom drag coefficient must be adjusted if the water depth approaches the bottom roughness height. How to perform this adjustment is an area for further consideration. In the present version of FVCOM, the effective water depth used in the bottom friction formulation is limited to 3 m; e.g., when the actual water depth is less than 3 m, the depth used in the bottom friction formulation is set to 3 m.

#### **3. Model Setup**

Basically, the model needs reasonable specifications of the following four items to obtain skillful predictability. They are (1) River boundary forcing conditions, (2) open ocean boundary conditions, (3) initial conditions, and (4) surface forcings. Each of these model elements is discussed below.

### *3.1. River Boundary Forcing Condition Specification*

There are seven rivers considered in the model. Traditional river discharge condition is used for the five small rivers that are not in the delta (the rivers with names in blue in Figures 2 and 3) area. These five rivers are the Petaluma River (2 m3 /s), Alameda Creek (3 m3 /s), Napa River (100 m3 /s), Coyote Creek (2 m3 /s) and Guadalupe River (3 m3 /s) with approximate mean annual flows in parentheses. The two rivers in the delta are Sacramento River at Rio Vista (1000 m3 /s) and San Joaquin River at Antioch (100 m3 /s) with approximate mean annual flows in parentheses. Two different upstream boundary condition types were considered for the Sacramento and San Joaquin Rivers forcing specification.

In type one, the average daily flow of Sacramento River was used to specify the flow at Rio Vista (RIO), while the San Joaquin River flow was estimated as the total delta outflow (OUT) minus the Rio Vista flow (RIO). The measured data are from the California Department of Natural Resources' DAYFLOW [32]. The average daily flows (note: negative flow indicates flow into the Delta from the Bay) are used only for the hindcast scenario testing. Minimum inflow and zero salinity were set up for the two rivers in low flow period when DAYFLOW's estimates may be suspect as noted by Oltmann [33]. In type two, the water level surface elevations (stage) were specified at Rio Vista and Antioch for Sacramento River and San Joaquin River, respectively, similar to the previous work in this region by MacWilliams *et al.* [16].

Both flow and stage river boundary conditions were used in the hindcast for comparison purposes. Please see Schmalz [22] for details. However, in the nowcast/forecast system we only use a river stage forcing for these two rivers. This is because NOCMP's top priority is to support PORTS for navigation safety, and water level prediction is paramount. Previous studies and personal communication with Michael MacWilliams [34] have found that flow boundary condition may be more suitable if the focus is on salinity prediction. However, our major concern is surface water level as in PORTS, and the stage boundary condition is more suitable.

#### *3.2. Open Ocean Boundary Condition Specification*

The open ocean boundary of the grid (see Figures 2 and 3) is forced with a superposition of the subtidal water levels and predicted tides. The harmonic constants of M2, S2, N2, K2, K1, O1, P1, and Q1 that are used to predict tide are derived from the Oregon State University Tidal Inversion Software (OTIS) for the West Coast (WC2010 1/30°) [35].

In the hindcast scenario, the subtidal water level signal at Point Reyes (see Figure 1 for its location) is used to prescribe the subtidal water level along the outer boundary. A revised sponge layer treatment at the open ocean boundary was considered. The salinity and temperature at the open boundaries were determined, with nudging, from NOAA's World Ocean Atlas 2001 [36]. No velocities are prescribed along the open ocean boundary.

In the nowcast and forecast scenarios, subtidal water level open boundary conditions are generated from the NCEP's (National Centers for Environmental Prediction) G-RTOFS (Global Real-Time Ocean Forecast System) gridded operational products. The temperature, salinity and baroclinic current open boundary conditions are also generated from G-RTOFS. The most recently available products for the given time period are searched and used for adjustments of the open boundary conditions using the COMF-HPC. Several horizontal interpolation methods are implemented, and a linear method is used for vertical interpolation from G-RTOFS vertical coordinates to model vertical coordinates. Measured real time sea surface elevation data at Point Reyes are used for the subtidal water level adjustment along the open boundary. Similarly, measured temperature data at San Francisco (see Figure 1 for its location) are used for boundary temperature adjustment. The adjustment is the difference between the observation and the G-RTOFS prediction at the start of the nowcast/forecast cycle and is discussed in greater detail in Section 6.1.
