A Water Balance Model to Estimate Flow Through the Old and Middle River Corridor

We applied a water balance model to predict tidally averaged (subtidal) flows through the Old River and Middle River corridor in the Sacramento–San Joaquin Delta. We reviewed the dynamics that govern subtidal flows and water levels and adopted a simplified representation. In this water balance approach, we estimated ungaged flows as linear functions of known (or specified) flows. We assumed that subtidal storage within the control volume varies because of fortnightly variation in subtidal water level, Delta inflow, and barometric pressure. The water balance model effectively predicts subtidal flows and approaches the accuracy of a 1–D Delta hydrodynamic model. We explore the potential to improve the approach by representing more complex dynamics and identify possible future improvements.


INTRODUCTION
Tidally averaged (hereafter referred to as subtidal) flow through the Old River and Middle River (OMR) corridor is an important metric for describing hydrodynamics in the interior Sacramento-San Joaquin Delta (Delta). As a result of south Delta water diversions, net flow through the corridor is typically in a landward (southerly) direction, except during times of high San Joaquin River inflow to the Delta. This so-called "reverse flow" affects Delta transport patterns and water residence times and thus has implications for water quality and ecology in the region (Glibert et al. 2014). Movement of water from north to south generally improves water quality in the OMR corridor by pulling high-quality water from the Sacramento River into the interior Delta. However, during periods of low net Delta outflow, this flow pattern tends to pull saline water from the western Delta into the interior. Salvage of the federally threatened Delta Smelt (Hypomesus transpacificus) in export facilities has been correlated with reverse OMR flows (Grimaldo et al. 2009). As a result, restrictions have been imposed on OMR flows as part of the U.S. Fish

Term Biological Opinion Reasonable and Prudent
Alternative (USFWS 2008) to limit the potential for smelt entrainment. Similarly, OMR flow restrictions are incorporated in the National Marine Fishery Service Biological Opinion to limit exposure of outmigrating winter-and spring-run Chinook Salmon (Oncorhynchus tshawytscha) and juvenile steelhead (Oncorhynchus mykiss) to the southern export facilities (NMFS 2009).
Because of the aforementioned restrictions, water managers need fast and relatively simple methods to accurately estimate OMR flows for short-and long-term operations and facility planning. A water balance model is an efficient and conceptually clear approach to meet these needs. Hutton (2008) developed a water balance model to estimate OMR flows, and provided a comparison with previously available statistical models (Snow 1986, unreferenced, see "Notes";Ruhl et al. 2006, unreferenced, see "Notes"). The Hutton (2008) model is coded in the Central Valley reservoir operations model (CalSim II) for long-term planning studies, and was adopted by the U.S. Bureau of Reclamation and California Department of Water Resources (CDWR) for operations planning and regulatory compliance on a demonstration basis (USBR 2014).
Water balance models are available in the CDWR DAYFLOW program to describe a variety of subtidal flows in the Delta (CDWR 1986). Notably, a water balance method estimates Delta outflow to comply with flow requirements that the State Water Resources Control Board (SWRCB 1999) imposes. This approach, referred to as the Net Delta Outflow Index, neglects changes in subtidal storage in the Delta by assuming that inflows and outflows balance on a daily basis. Oltmann (1998) compared this index with net flows estimated by flow monitoring and found it to be accurate at moderate to high flows but less accurate at low flows. Potential sources of error in the water balance method were cited as effects of the spring-neap cycle, variability in barometric pressure, and uncertainty in net channel depletions, herein referred to as Delta NCD (Oltmann 1998). Although not directly mentioned by Oltmann (1998), measurement error is also inherent in determining statistically significant small net flows in the presence of much larger estuarine tidal flows .
Simple water balance methods that assume a balance between inflows and outflows may be improved if changes in subtidal storage are considered. In addition to the spring-neap cycling and variation in barometric pressure Oltmann mentioned (which is also reported by Walters and Gartner [1985]), the magnitude of river inflows, regional and local winds, and hydraulic structure operations can also influence subtidal water levels. Because of the complex channel connectivity of the Delta, and the depth-and flowdependent effects of bottom friction, these subtidal water level forcing factors can interact to affect subtidal flows in non-linear ways.
The spring-neap cycle of subtidal water level, also referred to as spring-neap "filling and draining" (e.g., Stacey et al. 2010), has been widely observed in estuaries (LeBlond 1979). This variability is associated with compound tides, and occurs at frequencies related to those of astronomical tidal constituents (Parker 2007). For example, variation at the frequency of the principle lunar tide (M 2 ) minus the frequency of the principle solar tide (S 2 ), referred to as the compound tide constituent MS, is related to variations in tidal range over the spring-neap cycle and associated changes in subtidal friction (Buschman et al. 2009). This constituent has the same frequency as the astronomical constituent MSf (corresponding to a period of 14.77 days), but is created by shallow-water hydrodynamic effects, not astronomical forcing. Since these hydrodynamic effects are generated by bottom friction, they depend on river flow, tidal amplitude, and the non-linear interactions that develop between the two (Buschman et al. 2009;Godin 1999). A simple description of the spring-neap cycle of subtidal water levels is that the higher flow velocities during spring tides result in increased friction; an increased subtidal water level slope is, therefore, required to transport river water seaward (Buschman et al. 2009).
A related process that results in subtidal water level variability is Stokes drift, a landward flow which occurs in progressive wave systems from the temporal correlation between tidal currents and water depth (Stacey et al. 2010). This process and the associated compensation flow are substantial in the northern San Francisco Estuary (Stacey et al. 2010 tides results in an increased landward Stokes drift balanced by a seaward Stokes drift compensation flow which, like river flow, is driven by a subtidal water level slope (Jay and Flinchem 1997). Like river flow, Stokes drift compensation flow is also impeded by increased friction that results in spring-neap variability. Furthermore, Stokes drift and Stokes drift compensation flow are not always in balance. Sassi and Hoitink (2013) found substantial landward water flux in periods of peak Stokes drift that corresponded to spring tides. When river flow is present, its magnitude similarly influences water levels, subtidal friction, and the generation of compound tides. The overall subtidal water level variation, therefore, cannot be represented solely by a harmonic analysis (Jay and Flinchem 1997). Hydrodynamic models and complex analytical models (e.g., Buschman et al. 2009) can estimate these effects to accurately predict water levels. However, Buschman et al. (2009) reported that the pragmatic approach Godin (1999) proposed, whereby subtidal water level is estimated as a linear function of tidal range and net river flow, was also able to accurately reproduce observations. Not only do river flows influence the generation of compound tides, they also directly influence water level even in the absence of tides. Subtidal water level is further influenced by barometric pressure, local and coastal wind, and operations within the Delta. In South San Francisco Bay, Walters (1982) found that subtidal water level variations were generated by non-local coastal forcing, primarily related to barometric pressure, and that local wind contributed only a small amount of setup. A similar study by Walters and Gartner (1985) found similar forcings for subtidal water levels in San Pablo and Suisun bays. Operations that may influence water levels in the interior Delta include temporary barrier installation, diversions for the Central Valley Project (CVP), State Water Project (SWP), Contra Costa Water District (CCWD), and Delta NCD. Delta NCD is particularly uncertain (Siegfried et al. 2014) and may constitute a substantial portion of net flows during low inflow conditions. The water balance model for OMR flow requires estimation of subtidal flow division at channel junctions. Observed subtidal flow division depends on local water surface slopes, channel geometry and friction, and tidal amplitude (Buschman et al. 2010).
One reason for tidal variation is that Stokes drift and Stokes drift compensation flow, both of which vary with tidal amplitude, can be distributed unevenly in branching channels (Sassi et al. 2012). A portion of the water volume transported landward by Stokes drift in one channel may flow into an adjacent channel at a junction and return by a different pathway as Stokes drift compensation flow.
Observed flow divisions at a junction can change dramatically from temporary barrier installation. A barrier is typically installed at the head of Old River (HOR) in the fall and spring and is intended to benefit migrating San Joaquin River Chinook Salmon. When the HOR barrier is not in place, the net downstream flow at the Old River-San Joaquin River junction predominantly travels down Old River at low San Joaquin flows and is split approximately evenly at higher San Joaquin flows. With the HOR barrier in place, flow into Old River is restricted, and about 80% of the flow continues in the San Joaquin River. Temporary agricultural barriers are typically placed at three locations (on Old River, Middle River, and Grant Line Canal) during the summer months in order to raise water levels and keep local agricultural intakes underwater. These structures restrict flow, but allow some water over and through them, altering local water surface slopes and affecting flow splits.
In the remainder of this paper, we propose simple water balance models to estimate net flows in the Old and Middle rivers. To develop these models, we define a control volume in which the only unknowns are the flow divisions at two junctions, variation in storage, and OMR flow. We perform a water balance for this control volume in which the independent variables are south Delta diversions, barrier installation status, and San Joaquin River flow. The work we present in this paper builds on that of Hutton (2008) by adding the effects of subtidal storage in the water balance, and analyzes the model performance in detail. We examine the assumptions associated with the proposed approach and suggest possible future improvements. 4 VOLUME 14, ISSUE 2, ARTICLE 2 Q div = south Delta diversions: the sum of Delta exports and NCD in the control volume ΔV = change in water volume over time Δt All terms in Equation 1 have units m 3 s -1 . Flow is considered positive in the seaward direction; this is north for most channels and west for Indian Slough. South Delta diversions are considered positive when water is removed from the control volume.
Since long-term flow records are not available for Indian Slough and the San Joaquin River at Lathrop, we estimated them by linear regression with Q omr (Indian Slough), and Q vns and Q div (San Joaquin River at Lathrop), using the results of simulations performed with the Delta Simulation Model-2 (DSM2). The linear regression equations can be expressed as (2) (3) where a, aʹ, and bʹ are dimensionless fitting parameters, and c and cʹ are fitting parameters with units m 3 s -1 .
where A wb and B wb are dimensionless parameters, and C wb has units m 3 s -1 .
We performed two long-term DSM2 simulations, and used the results to estimate the fitting parameters in Equations 2 and 3. The first DSM2 simulation was historical, using observed values for boundary inflows and major diversions. We chose the timeperiod 1990 through 2012, the longest period with CDWR-verified boundary flow records and input files available at the time of this work. From this record, two time-periods were excluded: Jan-Feb 1997, because of flooding conditions on the San Joaquin River around Vernalis, and Jun-Dec 2004, because of the Jones Tract levee failure and subsequent pumpout. Both time-periods include anomalous flows into and out of the control volume that are not accounted for in Equation 1. The second DSM2 simulation was similar to the historical case, but did not include SWP and CVP diversions. Our intention for including this simulation data was to encompass a broader range of operational conditions in the regressions, so we could

Control Volume Approach to Estimating OMR Flow
We calculated OMR flow as the residual flow in a control volume centered on the south Delta ( Figure 1). Flow may enter or exit the control volume through river channels at the San Joaquin River at Vernalis, the San Joaquin River downstream of the HOR split (near Lathrop), Indian Slough, and Old River and Middle River at Bacon Island. Our motivation for defining the control volume in this way was to make use of the long-term flow records measured by the U.S. Geological Survey ( (CDWR 1995) drove the magnitude of these sources and sinks. Additionally, aggregate estimates of Delta NCD are made in real-time for compliance with net Delta outflow standards; these estimates are based on the difference between assumed Delta-wide gross channel depletions and real-time estimates of Delta precipitation, and are archived in CDWR's DAYFLOW program (CDWR 1986).
Conservation of fluid volume within the control volume dictates that, at a given time-step, inflows must be offset by outflows and changes in storage.
( evaluate non-historical operational regimes without relying on extrapolation in the water balance model. The specific version of DSM2 used in this study was v8.1.2, which underwent a full recalibration effort in 2009 (CDWR 2013). Model flow data were output at Indian Slough at Old River (DSM2 channel node 236); Old and Middle River at Bacon Island (channel nodes 106, 144, and 145); and the San Joaquin River at Lathrop (channel node 8). We tidally filtered raw 15-minute output data using a Godin filter to obtain net flows, and then daily averaged them to create a manageable number of data points for the two full 23-year time-series. We regressed Indian Slough flow against OMR flow because of their proximity and similarity in hydraulic behavior. We statistically related San Joaquin River flow at Lathrop to San Joaquin River flow at Vernalis by linear regression.
At high flows, a portion of the San Joaquin River flow upstream of the Old River junction spills over an overflow weir that connects the San Joaquin River to Paradise Cut. Because of the presence of this weir, we developed relationships of San Joaquin River flow at Lathrop to San Joaquin River flow at Vernalis for multiple ranges of San Joaquin River inflow. We accounted for the effect of barrier operation on the San Joaquin River-Old River flow split by 6 VOLUME 14, ISSUE 2, ARTICLE 2 obtaining different best-fit regression equations for filtered time-series when barriers were in place or absent. We considered specific cases with all barriers out, with the Grant Line Canal (GLC) barrier in and the HOR barrier out, and with the HOR barrier in. Because of different prevailing hydraulic conditions and construction designs, we treated the fall HOR barrier installation separately from the spring HOR barrier. We also included south Delta diversions in the regression because of their effect on local water surface slopes in all cases except the highest San Joaquin flows and when the spring HOR barrier is installed.
CDWR estimated Delta island diversion and return flows, and provided them as DSM2 boundary conditions (CDWR 1995). NCD in the south Delta control volume consistently averaged around 20% of the total Delta NCD ( Figure 2).

Subtidal Water Level Analysis
The final term in Equation 1 takes into account changes in subtidal storage in the control volume. Our approach followed Godin (1999), who estimated a linear effect of river flow on subtidal water level in addition to a periodic spring-neap influence. A linear effect of barometric pressure, acknowledged by Godin (1999) to have a significant effect on subtidal water level, was also included, resulting in: where η 0 is the subtidal water level in m, N i is the number of compound tidal constituents fit, A i is the amplitude of the ith tidal constituent, s i is the frequency of the ith constituent in radians day -1 , φ i is the phase of the ith constituent in radians, Q inflow is Delta inflow in m 3 s -1 , P is barometric pressure in millibars, a 0 is a fitting parameter in m -2 s, b 0 is a fitting parameter in m millibars -1 , and c 0 is a fitting parameter in m. We determined the empirical coefficients by fitting Equation 5 to water level records predicted by the DSM2 simulations at Old River at Bacon Island using a nonlinear least squares optimization approach (Levenberg 1944). While b 0 is known in theory and may be specified a priori by assuming hydrostatic pressure, Walters (1982) reported that the effect of barometric pressure observed on south San Francisco Bay subtidal water level is greater than expected. Therefore, this parameter was determined by fitting. We obtained pressure data from NOAA Station 9414290, located on the south side of the Golden Gate inlet. We supplemented this data series with measurements from the San Francisco International Airport before 1996, and where data gaps in the NOAA record existed. We developed the Delta inflow record by summing daily river inflows provided as DSM2 boundary conditions. Using a wind record from NOAA Station 9415144 we found wind effects to be negligible based on a correlative analysis of wind speed and stress components to the DSM2modeled water level. We note here that wind effects in DSM2 are accounted for implicitly through the use of observed stage data at Martinez for the water surface boundary condition, but no explicit wind surface stress is applied in DSM2. Therefore, we performed an additional correlative analysis using wind velocity components and USGS-observed subtidal water level at the Old River at Bacon Island station. A small correlation was found with east-west wind, indicating slightly increased water level with increased westerly wind. The effect was not large enough to justify the inclusion of wind effects in Equation 5.
Before fitting, a power spectrum analysis of water levels indicated three distinct amplitude peaks at periods greater than 25 hours, corresponding to the shallow water interactions of the K 1 and O 1 tides (constituent KO, period 328 hours), the M 2 and S 2 tides (constituent MS, period 354 hours), and the M 2 and N 2 tides (constituent MN, period 661 hours). Therefore, three amplitudes (A i ) and three phases (φ i ) were determined in Equation 5. For convenience the time origin (time = 0 days) of the estimated phases was taken to be Jan 1, 1900 at 00:00 in Pacific Standard Time (PST).
Equation 5 does not consider nodal factors to account for variations in tidal amplitude during the 18.61-year lunar node cycle. These node factors are important for the primary astronomical tidal constituents but are more ambiguous for compound tides. For simplicity, we neglect them here. To examine this assumption, we performed a harmonic analysis using the Vtide tidal harmonic analysis and prediction package (Foreman et al. 2009). We ran the Vtide package in analysis mode to calculate amplitudes for the spring-neap constituents. We then isolated these (all other tidal constituent amplitudes were set equal to zero) and the Vtide package was run in prediction mode to construct a water surface elevation time-series reflecting only the spring-neap tidal cycle. The improvement in fit to subtidal water level using Vtide was negligible, so we retained Equation 5 for conceptual simplicity and to allow simultaneous fitting of water level as a function of compound tides, Delta inflow, and barometric pressure.
The specific water surface elevation time-series used for the harmonic analysis was a 23-year (1990 through 2012) DSM2-predicted stage record at the Old River at Bacon Island station. We tested other locations throughout the control volume, and found results at this station similar to results at other stations located downstream of the temporary agricultural barriers. We analyzed the DSM2 stage time-series instead of the observed USGS stage at that location because of its long-term record without the complications of missing data.
We converted the subtidal water level predicted using Equation 5 to water volume using a relationship derived from a hypsographic curve of the southern Delta control volume: V = 14.916 × 10 6 * η 0 + 28.845 × 10 6 where V is the water volume in m 3 and η 0 is water surface elevation in m, NAVD88.
This approach implicitly assumes that subtidal water level is constant through the control volume.
The bathymetry data we used to derive this relationship were aggregated by CDWR from multiple bathymetric surveys (Wang and Ateljevich 2012). We calculated change in storage using centered differences.

A Direct Fit Approach to the Water Balance Model
The water balance method as presented is conceptually clear and founded on physical principles, including discrete configurations of the physical system (e.g., barrier installation or channel connectivity) and subtidal storage in the control volume. We performed separate statistical regressions for Equations 2, 3, and 5, which we then substituted into Equation 1 to obtain Q omr predictions for different flow and barrier installation cases. We did not fit the subtidal storage parameters to match Q omr directly, but to match a subtidal water level record, which we then converted to storage volume using Equation 6. Eighteen parameters are present before subtidal storage is considered, and we introduce an additional nine parameters to account for subtidal storage.
An alternative approach is to fit parameters to directly optimize the fit to OMR flow instead of developing regressions at individual junctions. In this direct fitting approach, nonlinear optimization of a single equation is used to estimate all relevant parameters. To derive this equation, we substituted Equations 2, 3, 5, and 6 into Equation 1, and combined parameters. We incorporated the effect of Paradise Cut by inclusion of a threshold flow at Vernalis, above which the slope of the Q omr dependence on Q vns changes. This guarantees a continuous relationship between Q omr and Q vns . We represented barrier effects as stepwise changes in the slope of the relationship between Q vns and Q omr . The resulting equation is: where A, Aʹ, B, A S , A f , A g , are dimensionless fitting parameters, A i ʹ and φ i ʹ are the unknown amplitudes in m 3 s -1 and phases in radians of the compound tide constituents, a 1 is a fitting parameter in days, b 1 is a fitting parameter in m 3 s -1 millibars -1 day, and C and D are fitting parameters in m 3 s -1 .
I s , I f , and I g are indicator functions for different barrier operations. I s takes a value of 1 during periods of spring HOR barrier installation and zero otherwise. I f is the analogous indicator function for fall HOR barrier installation, and I g is the indicator function for periods when only the GLC barrier is installed. OMR flows related to subtidal changes in control volume storage now depend directly on changes in Delta inflow and atmospheric pressure.
We made several simplifications in the derivation of Equation 7. We neglected terms, including barrier indicator function effects on the slope of the relationship between Q omr and Q div , and considered only one breakpoint in the slope of the Q vns and Q omr relationship. These simplifications reduced the number of fitting parameters from 27 to 16 by removing terms implicitly included in Equations 2 and 3, which are expected to have small effects on the OMR flow prediction. We fit the parameters of Equation 7 using the differential evolution optimization approach (Storn and Price 1997).

Model Performance Metrics
We compared predictions from the water balance approach and the DSM2 hydrodynamic model to USGS observed OMR flow. We obtained 15-minute discharge data directly from the USGS in April 2015 for the Old River at Bacon Island (USGS station number 11313405) and Middle River at Middle River (11312676) stations. We tidally filtered these data using a Godin filter and then daily-averaged them. Both records include periods during which a sensor was malfunctioning and no data were recorded. We generated a more complete record by developing a piecewise linear regression between the two stations and using it to fill in missing data. We also excluded from the regression analysis periods excluded from the DSM2 analysis for the same reasons. To account for different prevailing hydraulic conditions during strongly negative flows, we used a piecewise linear relationship (Figure 3). We determined the slopes and location of the breakpoint using the differential evolution non-linear optimization method (Storn and Price 1997).
We compared predicted and observed data on both a 5-day and 14-day running-average basis. For the empirical models, we used period averages of San Joaquin River flow at Vernalis and south Delta diversions to compute average OMR flow. For DSM2, we averaged predicted OMR flow over the period. The purpose of this averaging was to examine the accuracy of each method without the additional scatter caused by large day-to-day flow variations. Additionally, both averaging periods are important in regulatory contexts, including the Reasonable and Prudent Actions under the USFWS and NMFS longterm biological opinions.

Statistical Relationships for Ungaged Control Volume Flows
The relationship between Indian Slough and OMR flow is shown in Figure 4A. Subtidal Indian Slough flow averages about 6% of OMR flow. An implication of this is that a slightly less than 1:1 relationship exists between negative OMR flows and the magnitude of south Delta diversions. An improved regression could be developed by accounting for the effect of local NCD in Indian Slough (Hutton 2008). However, the local diversion term is omitted for simplicity and because its contribution to the accuracy of the Indian Slough flow estimate little affects the estimate of OMR flow.
The relationship between San Joaquin River flow at Vernalis and Paradise Cut is shown in Figure 4B.
There is a change point in the relationship at 467 m 3 s -1 (16,500 ft 3 s -1 ) when flow begins to spill over the overflow weir that connects the San Joaquin River to Paradise Cut. At 818 m 3 s -1 (28,900 ft 3 s -1 ), expansions in the San Joaquin River flow area and the geometry of the weir lead to smaller increases in Paradise Cut flow with increases in San Joaquin flow. We determined the exact location of these change points by fitting a piecewise linear function to the data using a non-linear least squares optimization method (Levenberg 1944). The resulting relationship forms the basis for different Q vns -Q lrp regressions based on San Joaquin River flow. During high flow conditions when Paradise Cut is active, the temporary barriers are typically not installed because of concerns of localized flooding. Figure 5 shows the effect of different barrier configurations and diversion levels on the San Joaquin-Old River flow split. Figure 5A shows the effect of the head of Old River (HOR) barrier. For San Joaquin flows at Vernalis below 468 m 3 s -1 , an approximately even flow split occurs when the HOR barrier is not installed. When the barrier is installed, flow into the Old River is restricted and San Joaquin River flow past Stockton is higher for a given San Joaquin River flow at Vernalis. The magnitude of this effect differs between the spring and fall barrier installations; a "full" barrier implementation is typically installed in the spring, and a "partial" barrier is typically installed in the fall.
During periods when the HOR barrier is not installed, the GLC barrier affects the San Joaquin-Old River flow split ( Figure 5B). When the GLC barrier is installed, it raises water levels in the south Delta, which influence the water surface slope near the junction and cause more water to flow down the San Joaquin River. All of the temporary agricultural barriers are usually installed and removed within a month of one another. Limited data during periods when the GLC barrier was installed and the Old and Middle river barriers were not installed suggests that the GLC barrier has a much larger effect on the split than either of the other two. The GLC barrier is also closer to the junction (~14 river km) than the Old River (~29 river km) or Middle River (~26 river km) barriers. For these reasons, the installation of the GLC barrier is treated as a different case in the Q vns -Q lrp regressions while the remaining agricultural barriers are not.
When neither the HOR barrier nor the Grant Line Canal barrier is installed, south Delta diversions noticeably influence the San Joaquin River-Old River flow split ( Figure 5C). As diversions increase, more flow is pulled into the Old River channel from the San Joaquin River. At very low Vernalis flows, diversions may even cause reverse flows in the San Joaquin River downstream of HOR. The influence of diversions on the flow split is also important for low San Joaquin River flow conditions during which the fall HOR barrier or the GLC barrier are installed (dependence not shown in Figure 5), and we considered this influence for the regression analysis. At higher San Joaquin River flow conditions, the influence of south Delta diversions on the flow split is less important.  Table 2.
Error metrics for DSM2 and the water balance model in predicting observed OMR flow are presented in Figure 7 and Table 3. DSM2 shows the highest accuracy, with 71% of 5-day-average predictions falling within ±15 m 3 s -1 (530 ft 3 s -1 ) of observed. The water balance approach without subtidal flow has 65% of predictions falling within ±15 m 3 s -1 . Both methods are off by greater than 35 m 3 s -1 (1,200 ft 3 s -1 ) only a small percent (1% to 4%) of the time. DSM2 predictions are generally more negative than observed, with 61% of model predictions having a negative residual. The water balance approach is not as biased, with 46% of predictions having a negative residual. When comparing to observed data on a 14-day-average basis, the short-term variations in subtidal storage are averaged out, and the water balance model without subtidal storage approaches DSM2 accuracy.
12 VOLUME 14, ISSUE 2, ARTICLE 2   Table 1 are compared against observed flow measured at Lathrop and Garwood Bridge

Subtidal Water Level Analysis
We calculated the amplitudes for the MS, MN, and KO spring-neap tidal constituents at the Old River at Bacon Island as 0.043 m, 0.026 m, and 0.047 m, respectively. Their phase lags, relative to Jan 1, 1900 at 00:00 PST, were 0.314 rad, 1.531 rad, and 2.880 rad. Estimated values for a 0 , b 0 , c 0 and were 0.0001 m -2 s, -0.0123 m millibar -1 and 13.67 m, respectively.
We examined the suitability of Equation 5 as a model equation. To confirm substantial effects of river flow and atmospheric pressure, we fit each parameter individually. After we fit subtidal harmonic variability (the summation term in Equation 5), we compared the residual subtidal water level to Delta inflow ( Figure 8A); a strong correlation (R 2 = 0.556) was found. After fitting both harmonic and flow variability, we then compared the residual subtidal water level to barometric pressure ( Figure 8B), and found a similarly strong correlation (R 2 = 0.388), indicating that both Delta inflow and pressure were suitable in Equation 5. Figure 8C shows the DSM2predicted and Equation 5-predicted subtidal water level for a representative year of the 23-year period. Each term of Equation 5 is added incrementally and has a significant effect on the estimated subtidal water level. The standard error of the estimated water level is 0.132 m when we considered only tidal harmonic variability, 0.088 m when we considered tidal harmonic variability and the Delta inflow effect, and 0.069 m for the complete Equation 5 including the barometric pressure effect. Several alternative parameters were considered, including natural logarithm and power law expressions for flow effects, regional and local wind, and south Delta diversions. None of these produced a substantial improvement in subtidal water level fit.
We used Equation 6 to convert the subtidal water levels predicted by Equation 5 to control volume storages and differenced them to calculate the final flow term in Equation 1. Table 3 and Figure 7 show results from the water balance model with the inclusion of the subtidal storage term. Fiveday-average model accuracy within ±15 m 3 s -1 is improved 3% to 68%, and standard error is reduced in all cases except for the highest San Joaquin flows. Fourteen-day-average model accuracy is very close to the water balance method without the inclusion of subtidal storage, since the change in storage term approaches zero as longer averaging periods are considered. Table 4 shows the parameters estimated by the direct fit optimization approach. Fitted slopes for the Q omr dependence on Q vns and Q div are similar to those derived for the incremental (multiple step) fit water balance model described previously. We also found a similar flow cutoff (parameter D in Equation 7) and change in Q omr dependence on Q vns at high flows.  Figure 7 suggest that the direct fitting approach, despite having a reduced set of parameters (16 instead of 27), has similar accuracy to the incremental fitting approach. The least accurate predictions were for periods with GLC barrier installation. In contrast to the incremental fit water balance, the direct fit, similar to DSM2, showed a tendency to underpredict OMR flows.

DISCUSSION
Quantitatively, the most important improvement we present over previous statistical models of OMR flow is the development of distinct flow division ratings for conditions with and without barrier operations. We found that ratings for the San Joaquin River-Old River junction varied with San Joaquin River flow and south Delta diversions, and could be represented well by piecewise linear functions. We also found that linear and continuous piecewise linear fits  by Gaspar and Ponte (1997) who found substantial correlations between wind-driven sea level variation and barometric pressure. Walters (1982) also reported a stronger than expected barometric pressure effect in south San Francisco Bay, consistent with our fitting results.
Significant error (a standard error of 0.069 m) remains in estimating the subtidal water level predicted by DSM2. Part of this error is likely from the simple relationships used to represent complex interactions between river flow and tidal flow. In addition, the Delta is a highly modified environment, and anthropomorphic effects are significant. The detailed timing of operations at Clifton Court Forebay and Jones Pumping Plant is not currently accounted for in the water balance approach, which uses dailyaverage boundary conditions. Furthermore, even if the subtidal water level in Old River at Bacon Island were predicted perfectly, other sources of error would remain in predicting subtidal storage. One is the assumption that the water level in Old River at Bacon Island represents water level in the south Delta. This appears to be a good first approximation, but some landward regions of the control volume are more fluvially influenced than Old River at Bacon Island, and areas near the southern export locations may experience more water level drawdown. A more accurate relationship between subtidal water level and volume derived, for example, from regressing control volume storage obtained from DSM2 to subtidal water level, could improve the storage estimate. The main advantage of using the DSM2 provided a good approximation of the DSM2predicted flows. There is evidence of a spring-neap signal in the flow residuals (the difference between DSM2 predicted flows and flows predicted by flow division regressions) at Indian Slough in particular. This is consistent with the finding of Sassi and Hoitink (2013) that Stokes drift and the Stokes compensation flow can be distributed unevenly in individual channels; one channel can feed water volume from Stokes drift into Stokes compensation flow in an adjacent channel. This spring-neap cycle can be substantial in some tidal rivers, but the flow residual at Indian Slough is typically less than 5 m 3 s -1 (Figure 4), suggesting that spring-neap effects in the relationship of OMR flow to flow at this location are weak.
Accounting for changes in subtidal storage in the south Delta control volume improves the prediction of subtidal OMR flow (Table 3). We found tidal harmonic variability in subtidal water level to depend primarily on three compound tide constituents, similar to studies in other estuaries (Godin 1999). Significant additional variability in water level is contributed by Delta inflow and barometric pressure. We represented each of these effects with a linear relationship, following Godin (1999). The estimated coefficient of proportionality for water level variability with barometric pressure, b 0 , was 0.0126 m millibar -1 , similar to the 0.01 m millibar -1 expected from the "inverted barometer" effect (Gaspar and Ponte 1997). Large variation from the expected value of 0.01 m millibar -1 with latitude was reported hydrodynamic model instead of the proposed water balance approach is improved prediction of storage in the south Delta.
The water balance approach of applying known flows and estimating unknown flows into the control volume using Equations 2 and 3 is conceptually clear. But noting that the individual regressions are then substituted into Equation 4, and that Equations 5 and 6 can also be substituted into Equation 4, a single equation can be derived and directly optimized to fit observed or predicted OMR flow. This approach allows all parameters to be fit in a single optimization step instead of through a series of linear regressions.
Equation 7 has conceptual advantages compared to the incremental (multiple-step) fitting approach.
Notably, this single equation shows the effect of all relevant parameters. The barrier effects can immediately be seen to be represented by a change in slope in the relationship of Q omr to Q vns . In addition, the relationship of Q omr to Q vns is continuous in Equation 7, while the incremental fitting approach has discontinuous relationship at two values of Q vns . Table 3 indicates that despite several simplifications introduced in Equation 7, which decrease the total number of parameters from 27 to 16, the overall standard error of the OMR flow predictions differs little from the incremental fitting approach. These simplifications include a single change in slope of Q omr to Q vns (as opposed to the two changes used in the incremental fitting) and barrier effects which were assumed to alter the Q omr to Q vns relationship but not Q omr to Q div . The performance of the simplified approach supports the assumption that the primary effect of varied barrier and flow conditions is change in the slope of the relationship of Q omr to Q vns . However, because the least accurate predictions of Equation 7 are for conditions with the GLC barrier, we conclude that the GLC barrier has some influence on the slope of the relationship between Q div and Q omr . A sensitivity test indicates that adding one additional parameter to represent this slope change can indeed slightly improve overall accuracy. However, since the difference was not large, we retained Equation 7 for conceptual simplicity.
Though the water balance approach can still be improved, the methods we present here are accurate. The water balance approach is a marked improvement over previous empirical approaches (evaluated in Hutton 2008), and we recommend its adoption in place of those currently in use. Its accuracy in predicting 5-day-average OMR flows approaches that of DSM2 and does not require a full hydrodynamic simulation of the Delta. Practically, the regressions in Table 2 provide a straightforward approach for managers to estimate OMR flows. The incorporation of the subtidal flow term requires more information, including forecasts of Delta inflow and barometric pressure, but these forecasts are typically available, and the subtidal flow can be readily estimated using Equations 5 and 6.
Because the errors shown in Figure 7 are approximately Gaussian, confidence intervals can be estimated using the standard errors shown in Table 3. The water balance approach is less accurate in predicting OMR flows during high San Joaquin inflow conditions. This is largely a result of error in predicting the Old River-San Joaquin River flow split. At more typical San Joaquin River inflows the 5-dayaverage standard error in the water balance estimate, without the subtidal flow term, is approximately 16.5 m 3 s -1 , indicating that OMR flows can be predicted with 95% confidence to within ± 33 m 3 s -1 . The inclusion of the subtidal flow term narrows the confidence interval to approximately ± 31 m 3 s -1 .
The proposed approach of analyzing flow divisions and accounting for subtidal storage has broad applicability to the Delta. We can readily envision several applications. One is to improve Delta outflow estimates by accounting for subtidal storage in the Delta. Another is checking the accuracy of estimated subtidal flow at monitoring stations by forming control volumes and accounting for subtidal storage within these volumes. This procedure could identify flow stations that require improved calibration, and quantify the uncertainty of observations. This would be particularly useful for calibrating hydrodynamic models. A demanding application that would only be possible with a highly accurate flow observation network is estimation of south Delta NCD using a water balance approach that incorporates observed flows and estimated storage.