Apparent Seasonal Bias in Delta Outflow Estimates as Revealed in the Historical Salinity Record of the San Francisco Estuary: Implications for Delta Net Channel Depletion Estimates

Accurate estimates of freshwater flow to the San Francisco Estuary are important in successfully regulating this water body, in protecting its beneficial uses, and in accurately modeling its hydrodynamic and water-quality transport regime. For regulatory purposes, freshwater flow to the estuary is not directly measured; rather, it is estimated from a daily balance of upstream Delta inflows, exports, and in-Delta water use termed the net Delta outflow index (NDOI). Field research in the 1960s indicated that NDOI estimates are biased low in summer–fall and biased high in winter–spring as a result of conflating Delta island evapotranspiration estimates with the sum of ungauged hydrologic interactions between channels and islands referred to as net channel depletions. In this work, we employed a 50-year observed salinity record along with gauged tidal flows and an ensemble of five empirical flow-salinity (X2) models to test whether a seasonal bias in Delta outflow estimates could be inferred. We accomplished this objective by conducting statistical analyses and evaluating whether model skill could be improved through seasonal NDOI flow adjustments. Assuming that model residuals are associated with channel depletion uncertainty, our findings corroborate the 1960s research and suggest that channel depletions are biased low in winter months (i.e., NDOI is biased high) and biased high in late summer and early fall months (i.e., NDOI is biased low). The magnitude of seasonal bias, which can reach 1,000 cfs, is a small percentage of typical winter outflow but represents a significant percentage of typical summer outflow. Our findings were derived from five independently developed models, and are consistent with the physical understanding of water exchanges on the islands. This work provides motivation for improved characterization of these exchanges to improve Delta outflow estimates, particularly during drought periods when water supplies are scarce and must be carefully managed.


INTRODUCTION
The Mediterranean climate of California is relatively dry and contributes to periodic shortages among the urban, agricultural, and environmental water-using sectors of the state. As a result, water managers widely acknowledge that scarcity of this valuable commodity calls for its careful management through the use of best available data and decision-support tools. The most recent multi-year drought (2012 through 2016) brought attention to limitations in data availability for decision-making and general understanding of the state's water availability, water use, and water rights. In 2016, the Public Policy Institute of California published a report (Escriva-Bou et al. 2016) that identified several gaps in the state's water information system, recommended that the state adopt an overarching goal of modernizing its water accounting, and outlined priority actions to accomplish this goal. In the same year, legislation was enacted (Assembly Bill 1755), which requires the California Department of Water Resources (CDWR) to develop protocols for water data sharing, documentation, quality control, public access, and promotion of open-source platforms and decision-support tools. Additional water datarelated legislation was enacted (Senate Bill 88) that requires reporting water diversions of 10 acrefeet (af) or more. Recent work, led by the state's Office of the Delta Watermaster, has focused on characterizing consumptive use of water in the Sacramento-San Joaquin Delta (Delta) (Medellín-Azuara et al. 2018).
Freshwater flow to the San Francisco Estuary (estuary), commonly referred to as Delta outflow, is a data metric of great significance to the state (see Figure 1). Control of the estuary's freshwater flow is a key aspect of Delta water management and affects many regions of California that depend on imported sources of water. Delta salinity control, implemented in part through management of freshwater flows to the estuary, has been a focus of attention for nearly a century (CDPW 1930(CDPW , 1931. Direct measurement of freshwater flow to the estuary is confounded by relatively large tidal flows and wide channels near the confluence of the Sacramento and San Joaquin rivers where flows exit the Delta and exchange with Suisun Bay. While advancements in direct outflow measurements will certainly continue, the current regulatory framework relies on a water balance metric known as the "Net Delta Outflow Index" (NDOI). NDOI is calculated as the sum of daily river inflows along the Delta's periphery minus water exports and in-Delta channel depletions (CDWR 2016). Channel depletions, which represent the sum of diversions from Delta channels onto adjacent islands minus in-Delta return flows and precipitation, is an ungauged quantity and is the focus of this work. The lack of channel depletion measurements is crudely resolved in the water balance methodology by assuming equivalency between crop evapotranspiration (or consumptive use) and channel depletions. This assumption, however, ignores the role of root zone soil moisture and other hydrologic processes in the overall Delta island water balance; i.e., crop consumptive use is not withdrawn instantaneously from surrounding channels. Determining the corresponding uncertainty in Delta outflow estimates is elusive (Monismith 2016). Because of the intimate relationship between outflow estimates and operational and regulatory actions, this uncertainty has statewide implications for the careful management of its water resources.
Following the work of Owen and Nance (1962), we hypothesize that the NDOI-based estimate of freshwater flow to the estuary may be seasonally biased because of its inherent conflation of in-Delta consumptive use and channel depletions. Figure 2 shows conceptually how a systematic offset between these flow terms would translate into a seasonal bias in the Delta water balance. In this work, we examine the occurrence of seasonal outflow bias through direct comparison with measured Delta outflow, and through indirect comparison with the estuary's historical salinity record. This evaluation was conducted with standard statistical measures of bias and utilized a 5-decade-long record of salinity and outflow data in conjunction with a suite of published flowsalinity models that were re-calibrated with a consistent data set (Rath et al. 2021, this volume).

BACKGROUND Physical Setting
The geographic focus of this paper is the central and northern reaches of the estuary and the Delta (Figure 1). Compared to other deltaic regions of the world, the Delta is unusual in having formed significantly inland from its ocean outlet. It is the entry point of over 90% of the freshwater inflow to San Francisco Bay (Cheng et al. 1993) and is drained by a vast 194,000 km 2 watershed, including the high-elevation Sierra Nevada mountain range and the Central Valley (US EPA 2012). The Delta is a major source of freshwater for the state of California, and water is withdrawn from several intakes to supply agricultural and urban areas throughout the state. San Francisco Estuary, the largest on the Pacific Coast of the western hemisphere, is home to one of the region's most important ecosystems (with 750 species of animals and plants) and serves as a major stopping point for migratory birds on the Pacific Flyway (Luoma et al. 2015).
For millennia, river flows and tidal action have deposited sediments in the Delta, the low point of the Central Valley. As these sediment deposits covered tules and other plants, thick organic peat soils were formed. Beginning in the mid-to late-19th century, most of the region (3,000 km 2 ) was drained, diked, and developed for agriculture. After these reclamation activities, wind erosion and oxidation of organic peat soils have caused  VOLUME 19,ISSUE 4,ARTICLE 4 a steady loss of surface elevation. The elevations of most interior Delta islands are now below the level of the surrounding waters, promoting seepage from adjacent channels and posing a major flood risk during levee failures events.
Under normal conditions, the low land elevations permit irrigation diversions through siphons with minimal energy use. To drain water from the crop root zone, these low-elevation islands use pumps to return water to adjacent channels. Delta islands are also hydrologically influenced by crop evapotranspiration, seepage through levees, and interactions between local groundwater and the root zone. Historically, few of these hydrologic exchanges have been directly measured. However, based on recent state legislation that required measurement of all diversions greater than 10 af, recent efforts have been undertaken to quantify the diversions to selected islands (2019 draft report by Heidel et al., unreferenced, see "Notes"; 2019 draft presentation by MBK Engineers prepared for the Metropolitan Water District of Southern California, unreferenced, see "Notes"). Also, the Delta Watermaster is leading an effort to estimate Delta crop evapotranspiration through remotesensing techniques (Medellín-Azuara et al. 2018). Based on annual totals, channel depletions (the difference between island diversions and returns) appear to be a relatively small part of the average Delta water balance, accounting for approximately 6% of the total output over water years (WYs) 1998-2009(Ariyama et al. 2019. The California water year begins on October 1 of the preceding calendar year. However, the relative Figure 2 We hypothesize that the NDOI-based estimate of Delta outflow may be seasonally biased because of its inherent conflation of in-Delta consumptive use and channel depletions. This conceptual chart shows how a systematic offset between these flow terms would translate into a seasonal bias in the Delta water balance. https://doi.org/10.15447/sfews.2021v19iss4art4 influence of channel depletions on the Delta water balance is magnified during drought years, representing 17% of Delta water use in WY 2015 ( Figure 3). Even in more ordinary years, channel depletion estimates for summer and fall months are significant in magnitude (typically 30% to 80%) compared to Delta outflow.
Salinity intrusion is a natural phenomenon in estuaries and depends on freshwater flow, tidal variations, regional wind patterns, and the system's shape and bathymetry. The seasonal pattern of early 20th-century Delta salinity intrusion was indirectly documented by a local industrial facility through its records of barge travel up the Sacramento and San Joaquin rivers to collect freshwater for its operations . Freshwater flow is the most important factor in repelling salinity intrusion and drives the salinity gradient across the estuary. Salinity ranges from near-ocean concentrations at the mouth of San Francisco Bay at Golden Gate to near-freshwater concentrations in the Delta channels, with an intermediate zone of salinity that moves landward and seaward seasonally as a function of freshwater outflow. The estuary and upstream watersheds are highly developed and, as a result, freshwater outflow is heavily influenced by withdrawals and upstream water management and use. The influence of secondary factors on Delta salinity intrusion, including tidal variations and regional wind patterns, has been examined using threedimensional (3-D) hydrodynamic models (e.g., Martyr-Koller et al. 2017;Gross et al. 2009).

Regulatory Setting & Definition of Net Delta Outflow Index (NDOI)
The first Delta flow regulation, proposed in 1931 (CDPW 1931), was crafted in response to water quality degradation (i.e., salinity intrusion) that had been observed in the early years of the 20th century. This water quality degradation was associated with upstream agricultural development and other regional hydrologic alterations that had taken place and was exacerbated by an extended period of drought. At that time, it was assumed that a chloride concentration greater than one part per thousand (ppt) was unsuitable for irrigation, leading state engineers to recommend that a minimum Delta outflow of 3,300 cubic feet per second (cfs) be maintained year-round 1 . The proposal was incorporated into the California State Water Plan but was not incorporated into the later operations plan of the federal Central Valley Project (Jackson and Patterson 1977).
After the original 1931 proposal, flow regulations, with the goal of protecting multiple beneficial uses of water in the Delta, have continued to evolve. Regulatory activity in the region is led by the California State Water Resources Control Board (CSWRCB), an agency concerned with water quality and water rights adjudication in California. In August 1978, the CSWRCB adopted its Delta Plan and Decision 1485, which set objectives for Delta outflow as well as flow and salinity in interior Delta channels 1. Placing this recommendation in a modern context, a one ppt chloride (equivalent to a specific conductance of 3.6 mS/cm) isohaline indicates a relatively high salinity, and its longitudinal position would be downstream of the congruent X2 isohaline, which is defined as 2 ppt bottom salinity or 2.64 mS cm -1 surface salinity. A year-round Delta outflow objective of 3,330 cfs would hold a steady X2 position of approximately 95 km, as estimated by the empirical relationship reported in Hutton et al. (2015). Such an objective would allow salinity intrusion much greater than current regulations, which hold X2 between 75 and 81 km, with approximate outflows between 7,100 cfs and 11,400 cfs in spring months.

Figure 3
The relative influence of channel depletions on the Delta water balance is magnified during drought years, representing 17% of Delta water use in WY 2015 VOLUME 19, ISSUE 4, ARTICLE 4 (CSWRCB 1978). The CSWRCB continued to pursue additional regulations for protecting the ecosystem, and through a cooperative effort with the US Environmental Protection Agency, nongovernment organizations, and other federal and state agencies, developed an agreement in 1994 termed the Bay-Delta Accord. Consistent with this accord, the CSWRCB updated its Delta Plan in 1995 and adopted Decision 1641 in 2000 (CSWRCB 2000), which still applies today. According to this regulation, the low-salinity zone in the estuary, specifically the location of the 2 ppt bottomsalinity isohaline or X2, is used as the basis of outflow management in the Delta (Jassby et al. 1995;CSWRCB 2000) during several months of the year. It is worth noting that this regulation may be met through maintaining either salinity or outflow constraints, with the latter constraint based on an assumed relationship with Delta outflow and therefore net channel depletions. X2 position is measured as the distance in kilometers upstream from the Golden Gate. Additional X2 restrictions were imposed through a Biological Opinion rendered by the US Fish and Wildlife Service (USFWS 2008) to protect endangered species. The regulatory environment continues to evolve through updated Biological Opinions and anticipated updates to CSWRCB standards (Fleenor et al. 2016;CSWRCB 2018;USFWS 2019).
The current regulatory definition of Delta outflow relies on the NDOI water balance metric, which is calculated as the sum of daily river inflows along the Delta's periphery minus water exports and in-Delta channel depletions (CSWRCB 1978(CSWRCB , 2000:

Net Delta Outflow Index (NDOI) = Delta Inflows -Delta Exports -Net Delta Channel Depletions
(1) The first two terms in the water balance equation above (i.e., inflows and exports) are directly measured, albeit with some imperfections and omissions. The third term (net Delta channel depletions or NCD) is an ungauged quantity and represents the difference between gross channel depletions and precipitation:

Net Delta Channel Depletions = ∑ (Channel Diversions -Channel Returns) -Delta Precipitation
(2) The magnitude of NCD is positive (by definition) when water is removed from channels onto islands in aggregate. NCD can be a significant component of the NDOI calculation, particularly during periods when Delta inflows are relatively low. Like NDOI, challenges associated with direct measurement have rendered NCD an estimated quantity. But unlike NDOI, the water balance terms used to estimate NCD are virtually ungauged. Channel diversions and returns (i.e., drainage) occur in a highly distributed manner, including diversions from Delta channels through approximately 1,800 siphons (CDWR 2016), drainage flows that return to Delta channels from approximately 200 pump stations (CDWR 1995), and precipitation falling onto approximately 3,000 km 2 calculated from measurements at a single location. Additional ungauged hydrologic processes contribute to channel depletions, including seepage from channels through levees, upward flow from groundwater into the root zone, and evaporation from Delta channels.
The lack of directly measured NCD water balance terms is crudely resolved in the NDOI-based methodology by assuming that a single, annually repeating pattern of crop evapotranspiration (or consumptive use) represents gross channel depletions (i.e., channel diversions minus channel returns). We use the terms evapotranspiration and consumptive use interchangeably in this work because crop water needs are generally satisfied in the Delta. The assumed consumptive use pattern is based on a data set developed in the 1960s and thus reflects land-use and wateruse conditions of that era. Land use, particularly in the Delta uplands (i.e., periphery areas within the legal Delta), has changed dramatically since the 1960s (2021 email correspondence between PHH and G. Gartrell, unreferenced, see "Notes"), underscoring the need to reassess Delta-wide consumptive use patterns. Assuming equivalency between consumptive use and channel depletions ignores the role of root zone soil moisture and other hydrologic processes in the overall Delta island water balance; i.e., crop consumptive use needs are not withdrawn instantaneously from surrounding channels.
Historical Efforts to Understand the Inter-Connected Role of Net Delta Channel Depletions, Outflow, and Salinity Intrusion Since the 1920s, if not earlier, California water scientists and engineers have understood the fundamental relationships between net Delta channel depletions, outflow, and salinity intrusion. The California Department of Public Works (CDPW 1931), a predecessor to the CDWR, began collecting data and conducting studies to estimate flows required to repel salinity intrusion into the Delta. Bulletin 27, a synopsis of this work, was devoted to compiling the best available understanding of how salinity intrusion was influenced by Delta outflow, which is influenced by in-Delta water use. As reported in Bulletin 27, detailed data analyses attempted to correlate flows and salinity at discrete locations in the Delta. Several relationships were evaluated, including maximum and minimum annual salinity versus annual unimpaired runoff, date of salinity advance versus annual runoff, maximum annual salinity versus summer Delta inflow volume, and maximum annual salinity versus Delta inflow rate. None of the relationships explicitly recognized the role of antecedent flow conditions, contributing to their modest explanatory power.
Bulletin 27 also reported work undertaken to quantify consumptive use in the Delta, arriving at an annual Delta-wide volume estimate of 1.1 million af, and seasonal flow estimates varying between 400 cfs in the winter and 3,700 cfs in the summer. This early work, undertaken in cooperation with the US. Department of Agriculture, included land-use surveys and evapotranspiration studies using tanks and field studies (CDPW 1930). The Bulletin recommended a flow of 3,000 cfs (3,300 cfs with a safety factor) past Antioch to control salinity at a maximum of 1 ppt chloride. Thus, accounting for Deltawide consumptive use, the Bulletin concluded that required Delta inflows would vary between 3,700 cfs in winter and 7,000 cfs in summer.
This conclusion implicitly assumes equivalence between consumptive use and net channel depletions. Owen and Nance (1962) acknowledged that the best available science of their era was unable to discern consistent relationships between salinity intrusion and computed Delta outflow. While they had limited understanding of how salinity intrusion responded to a complex time history of outflow, they correctly intuited that the lack of correspondence between salinity intrusion and computed outflow resulted (at least in part) because consumptive use does not fully represent the way water is depleted from Delta channels.
To test their hypothesis, the authors conducted a field study on Twitchell Island (see Figure 1) and showed that soil moisture must be accounted for to correctly estimate channel depletions. They monitored surface inflow, drain discharge, precipitation, changes in soil moisture content, weather data, and cropping patterns and methods as part of the field study. They identified siphon inflow, seepage, and precipitation sources of water for the island water balance; crop evapotranspiration and pump return flows were identified as sinks of water for the island water balance; soil moisture and ground water storage were identified as additional components of the island water balance. They concluded that (1) monthly consumptive use is not synonymous with monthly channel depletions in the Delta, and (2) the effect of change in soil moisture on computed Delta outflow is to increase the computed value of outflow in low flow months and to decrease the computed value of outflow in high flow months.

Critique of Net Delta Outflow Index (NDOI)
Although NDOI remains the prevailing outflow estimate for regulatory purposes in the Delta, a variety of field and modeling studies have been performed since the mid-1990s to examine its effectiveness over different time-periods and flow conditions. The US Geological Survey (USGS) has installed a flow monitoring network in the western Delta that allows direct measurement of Delta outflow and NDOI to be compared. Four gauging stations-installed at Rio Vista, Threemile Slough, San Joaquin at Jersey Point, VOLUME 19,ISSUE 4,ARTICLE 4 and Dutch Slough ( Figure 1)-report tidal flow at 15-minute intervals using a suite of hydroacoustic tools (including acoustic Doppler current profilers and ultrasonic velocity meters) (Burau et al. 2016). This monitoring network provides an indirect measurement of the net Delta outflow just upstream of the confluence of the Sacramento and San Joaquin rivers. Monismith (2016) compared NDOI and the direct USGS flow measurements for WYs 2008 through 2014. Over the full range of flows, he found the two flow estimates to be in good agreement. He found a linear relationship between the two flows with an R 2 of 0.94 and a root mean square difference on the order of 5,000 cfs. The author considered the difference-of similar magnitude or larger than reported low-flow NDOI values-to be related to the error in the NDOI net channel depletion estimates. Importantly, however, when the fitting was limited to flows less than 10,000 cfs, there was virtually no correlation between NDOI and direct measurement (R 2 = 0.05). Monismith (2016) noted that because the actual Delta outflow is unknown, it was difficult to say which approach (direct measurement versus water balance) was more accurate. One contribution of the present work is to show that the residuals between NDOI and USGS flow measurements exhibit considerable coherence and seasonality when organized as time-series; these observations are not apparent in the scatter view. CDWR (2016) opined that, despite improvements in the science of estuarine flow measurements, a water balance approach (e.g., NDOI) remains the best available science and most appropriate for estimating Delta outflow within a regulatory framework. This opinion was supported by acknowledging practical challenges associated with interpreting temporally resolved outflow measurements in a management setting that requires multi-day travel times between upstream reservoirs and the river confluence; these challenges are confounded by wide channels, large tidal oscillations, and tidallyinduced filling and draining events. We note that the USGS flow data in the western Delta are actually derived from velocity and water-level measurements using an involved calibration process that is affected by tides and wind-strength differences. This process requires only a modest bias between flood and ebb to produce errors on the same order as the mean measurement. Furthermore, CDWR (2016) concluded through an analysis of data that actual outflow may be higher during the irrigation season and lower during winter months. A peer review of CDWR (2016) acknowledged the complexity of Delta outflow measurement and recommended a broad collaborative monitoring and analysis framework to advance the science (Fleenor et al. 2016). We note that the NCD estimate used in the NDOI water balance calculation is not widely used in other Delta water-management activities; most contemporary hydrodynamic modeling done in support of permitting and planning in the Delta uses updated methodologies (as described in the following section), which incorporate changing land use and atmospheric forcing.

Estimation of Delta Island Hydrologic Exchanges
The water balance terms for a typical Delta island are shown conceptually in Figure 4. Although the availability of field water-balance measurements is limited, methodologies have been developed to refine the NCD estimate used in the NDOI water balance calculation. CDWR (1995) developed the Delta Island Consumptive Use (DICU) model to estimate Delta agricultural diversion and return volumes on a fine spatial level compatible with available hydrodynamic and water-quality transport models. Discretizing the Delta into 142 sub-areas, consumptive use and diversion/return volumes are calculated by region using monthly time-steps that account for crop evapotranspiration, crop type and acreage, precipitation, seepage, applied irrigation water, soil moisture storage, applied leach water, and surface runoff. The DICU model assigns channel diversion and return volumes to nodes associated with the CDWR's Delta Simulation Model (DSM2), a process-based, 1-D hydrodynamic and waterquality transport model (CDWR 2020 This model builds upon the DICU framework and adopts the daily calculation of crop evapotranspiration as incorporated in CDWR's Delta Evapotranspiration of Applied Water (DETAW) model (Kadir 2006). Development of the DCD model was motivated by the finding that DSM2 model simulations, when informed by the DICU model, do not adequately capture salinity intrusion patterns during certain periods. Based on the insight that NCD is not synonymous with consumptive use (both in timing and magnitude), the DCD model considers other flow terms in the calculation of net channel depletions, notably seepage through levees and interaction with groundwater. The role of groundwater in contributing to NCD has been recognized previously, but not separately quantified. As part of the CDWR's ongoing model evaluation, the DCD and DSM2 models are being integrated to better represent salinity at key locations across the Delta (Liang and Suits 2018b). Other islandscale modeling efforts have also been proposed (Siegfried et al. 2014; 2019 draft report by Heidel et al., unreferenced, see "Notes"), although these are limited by the availability of measured data for calibration.

Other NDOI Accounting Issues
Critiques of the NDOI water balance calculation focus primarily on channel depletions as a source of systematic error and seasonal bias. However, other elements of the water balance have not been fully modernized and may contribute to NDOI error and bias. An example is Yolo Bypass inflow to the Delta, which is computed as the sum of perimeter inflows that are exhausted before they reach the Delta in the summer. A typical NDOI summer estimate for Yolo Bypass inflow is several hundred cfs; however, observed flow at Lisbon Weir is usually negative (i.e., away from the Delta) in summer, except for recent flow actions and experiments. Calaveras River inflow-a component of NDOI's eastside inflow-provides another example of a minor Delta inflow that is typically exhausted before it reaches the Delta in the summer but contributes to the NDOI summer estimate. We also note that the Dayflow Delta outflow water balance (https://water.ca.gov/Programs/Environmental-Services/Compliance-Monitoring-And-Assessment/ Dayflow-Data), while not strictly identical to the NDOI water balance, omits Sacramento Regional Sanitation District discharges to the Delta. Except for these discharges, the NDOI accounting issues noted here involve overestimation of surface water sources which, once corrected, would reduce outflow in summer months. By contrast, in this work we show evidence that outflow in late summer is higher than currently estimated in the NDOI water balance. These accounting issues are  not an alternative explanation for higher summer outflow; rather, the biases shown in this work would have to be even larger to compensate for such accounting issues.

METHODS
Following the seminal work of Owen and Nance (1962), we hypothesize that the NDOI-based estimate of freshwater flow to the San Francisco Estuary may be seasonally biased by its inherent conflation of in-Delta consumptive use and net channel depletions. We explore the occurrence of seasonal outflow bias through direct comparison with measured Delta outflow and through indirect comparison with the estuary's historical salinity record. We accomplished this latter comparison by comparing observed salinity (X2) with the isohaline position as predicted by a suite of flow-based empirical models (Jassby et al. 1995;Monismith et al. 2002;MacWilliams et al. 2015;Hutton et al. 2015;and Monismith 2017) that were consistently recalibrated by Rath et al. (2021, this volume). Regarding this latter comparison, we acknowledge that tides and wind influence the relationship between outflow and salinity intrusion, and could confound the assessment of seasonal bias. Our work considered the effect of tides by comparing results with an empirical X2 model that includes tidal inputs (Rath et al. 2017). However, wind effects on salinity could not be similarly considered because of the absence of a modeling framework that could be run over decadal time-scales. Methods used to undertake this work are summarized below.

Data
In this work, we used a 50-year record of Delta outflow and salinity that spanned WYs 1968 through 2017 as reported in a companion paper , this volume). The outflow record represents daily NDOI values approximated by the Dayflow Delta outflow water balance. We also used high-resolution (15-minute) outflow data measured by a four-station flow monitoring network maintained by USGS and located near the confluence of the Sacramento and San Joaquin rivers (Burau et al. 2016). We averaged the 15-minute gauge data, available from April 1996 through September 2017, to create a daily outflow record that is considerably shorter than the 5-decade NDOI record. We computed an outflow residual data set as the daily difference between gauged and NDOI outflow values. The salinity data (X2) represent daily average estimates of isohaline position as derived from surface salinity measurements in the estuary and as generated from model simulations. The observed X2 record is based on work reported by Hutton et al. (2015); the simulated X2 record is based on a suite of recalibrated empirical models reported by Rath et al. (2021, this volume). We computed empirical X2 model residual data sets as the daily difference between simulated and observed X2 values.

Statistical Analysis
We examined the outflow residual data set and the X2 model residual data sets for seasonal signals using a variety of statistical analyses. We tested both data sets using the Seasonal and Trend Decomposition by Loess (STL) method (Cleveland et al. 1990); this method fits smoothing splines of differing time-scales to extract long-term (trend), medium-term (seasonal), and short-term (remainder) signals from time-series data. We also examined the model residual data sets for seasonal signals using standard statistical metrics including box-and-whisker plots, the square of the correlation coefficient (R 2 ), the Nash-Sutcliffe efficiency (NSE) index, the mean residual, and the root mean square error (RMSE).
Outflow and X2 data were partitioned in several ways to conduct the statistical tests for the 50-year period of record, including bins that represented months and bins that represented different hydrologic ranges as measured by current day outflow and antecedent (previousday) isohaline position. The empirical X2 models tend to perform best under conditions that represent an inverse relationship between X2 and outflow. Conditions that are inconsistent with this relationship (e.g., a high outflow event after an extended period of high salinity intrusion) tend to be poorly replicated by the empirical X2 models, and reflect highly dynamic conditions not easily encapsulated in a single outflow time history. To differentiate empirical X2 model performance accordingly, we used a nine-bin data partition (see Table 1) following Rath et al. (2021, this volume).

Alternate Outflow Time Series Analyses
We imposed alternate outflow time-series on the recalibrated empirical X2 models reported in Rath et al. (2021, this volume) to draw further inferences about seasonal bias in the NDOI record. We used two distinct approaches, which are summarized below.
In the first method, we constructed an alternate outflow time-series (spanning April 1996 through September 2017) by assimilating the decomposed seasonal signal from the outflow residual data set described above. We hypothesized that a more accurate seasonal signal may be embedded in the USGS gauge data that is not captured by the water balance approach used in the NDOI estimate. Specifically, we constructed the alternate timeseries by adding the daily average of the STL seasonal component to the daily NDOI timeseries. We then recalibrated the empirical X2 models with this new outflow time-series using the same approach as in Rath et al. (2021, this volume) and compared their performance against those calibrated with the unmodified NDOI timeseries. Results from this first method are reported in Results (Seasonal Bias Revealed by Observed Outflow Data section).
In the second method, we used a data-driven statistical approach that assimilated optimized monthly bias adjustments into the original NDOI outflow time-series. Results from this second method are reported in Results (see "Seasonal Bias Revealed by Data-Driven Flow Adjustments").
For each empirical X2 model, we fit a Bayesian statistical model for 12 monthly NDOI adjustments and applied these 12 fitted parameters to each year in the data set. We constrained these monthly adjustments to approximately sum to zero; this constraint was implemented by adding a tight prior distribution on the sum of the 12 monthly adjustment values. Our rationale for imposing this constraint was that, while the annual outflow and NCD volumes may be reasonably accurate (or difficult to correct using inverse modeling as demonstrated by Rath et al. [2021, this volume]), the seasonal distribution of those volumes may be misallocated. The Bayesian method gives a probability distribution for each of the 12 monthly values instead of one best-fit value, but the average of the distribution can be used when a single value is preferred. The Bayesian method requires specification of a prior distribution, which is a constraint on the fitted value that is selected to be as strong (highly constraining) or as weak as necessary. Specification of prior distribution is advantageous when overfitting may be of concern, and when it is useful to incorporate prior knowledge. The method allows for explicit consideration of parameter uncertainty and for the ease of setting up the constrained fit in an elegant way. We obtained estimates of the posterior distributions associated with the monthly outflow adjustments with Markov Chain Monte Carlo (MCMC) sampling-a numerical, iterative procedure that gives samples from the probability distributions of the 12 monthly adjustments. Details on the data-driven bias adjustment methodology are provided in the supplemental information (Appendix A).
One could imagine recalibrating the salinity model parameters again (or even jointly optimizing both flow offsets and model parameters) to this new data set as was done in the first method (i.e., NDOI modified by USGS data). However, our approach was guided by an interest in investigating whether specific small, patterned perturbations of the underlying NDOI data can better explain the salinity measurements. In contrast, the first method was an independent observation-driven revision of the flow data. Updating both sets of parameters as part of the data-driven approach-salinity model and 12 monthly flow adjustments, iteratively or jointly-may have resulted in a superior calibration fit but would have hindered interpretability and risked issues associated with overfitting.
The monthly outflow adjustments derived from the data-driven statistical approach were calibrated from a subset of the 50-year record of Delta outflow and salinity. Specifically, calibration data were limited to three of the nine data bins highlighted in yellow in Table 1 (i.e., the bins most representative of an inverse relationship between X2 and outflow). Additionally, extremely low outflow events (< 350 cfs), extremely high outflow events (> 88,000 cfs), and days where observed or antecedent X2 was less than 45 km were excluded from the calibration set due to instabilities with some models' representation of very high or very low flows. This screening resulted in monthly fitting data subsets that ranged between 762 and 1,214 points, with the high-outflow months of December through April having the most excluded data. Once we determined the flow adjustments, we applied them to the full 5-decade NDOI time-series to evaluate how the empirical X2 models performed compared to predictions that used the original NDOI time-series.

RESULTS
Results from the statistical and alternate outflow time-series analyses are presented below. The various analyses uniformly reveal seasonal bias in the NDOI record.

Seasonal Bias Revealed by Observed Outflow Data
As mentioned above, Delta outflow estimates provided by the USGS gauge data and NDOI water balance are strongly correlated over the full annual range of flow. Unfortunately, discrepancies introduced by channel depletion uncertainty, accounting errors, gauging issues, and signal processing are at their greatest during summer months, when outflow is at its lowest and accuracy is most critically needed for modeling and analysis.
Some of the discrepancy between outflow estimates is tractable by tracing residuals over time. Figure 5 presents time-series of tidallyfiltered (Godin 24-25-25) USGS outflow estimates, NDOI water balance estimates (top panel), and associated residuals (bottom panel). Residuals exhibit a clear time trend, with USGS estimates low relative to NDOI in the spring and during large winter storms, and high (relative to NDOI) in the summer, with a zero crossing in July from negative to positive and back to negative in early January. The USGS outflow estimates exhibit a structural break in 2007 when the seasonal magnituFdes and mean shift significantly. Another more subtle change occurred in 2013. These changes may have been caused by re-calibrations of the Jersey Point station ratings near those times.
Given that flow and season are easily conflated, it is possible that the residual trend shown in the bottom panel of Figure 5 is produced by flowdependent bias rather than seasonal bias. A possible explanation for a flow-dependent bias is that flow station calibration periods are very brief compared to sub-tidal periods, and the fitting procedures for the stations favor fidelity to tidal frequencies. To test the robustness of this explanation, we fit a simple linear scaling to the flow data measured at the Jersey Point and Rio Vista stations, plus a combined scalar shift for the sum of all the stations using least squares to dates with NDOI below 40,000 cfs and with candidate structural breaks (changes of parameters) on the parameters on the dates in October 2007 and 2013 when station re-ratings occurred. Residuals between USGS and NDOI estimates become smaller after applying this correction to the USGS data (RMSE drops from 3,302 cfs to 2,540 cfs for flows under 40,000 cfs)-particularly in the early record-but the seasonal shape of the time-series plot did not change (see bottom panel of Figure 5). While this is apparently a significant effect, it does not adequately explain most of the observed residual pattern. Because such a re-scaling complicates the broader discussion, we did not apply it to the remainder of the work presented here.
To analyze the residual time-series' seasonal component, we applied the STL decomposition (Cleveland et al. 1990) to disaggregate the Delta outflow residual into long-term, seasonal, and high-frequency components. In pre-processing the data, we generally filled missing data by interpolation. In the case of Dutch Slough, we filled missing data by assuming zero net flow. More effectively, we applied a third-order SavGol filter with an 81-day window, so that fortnightly and monthly fluctuations sympathetic to the spring-neap cycle were suppressed. This essentially adds a preliminary low-pass filtering to the STL procedure, a somewhat sub-optimal way of achieving a low-pass step as in MSTL, the multiple band variant of STL.
Output from the STL analysis is shown in Figure 6; the trend, seasonal, and remainder components are shown in the top, middle and bottom panels, respectively. A close-up of the seasonal component is provided in Figure 7; this pattern is repeated for all years. The seasonal pattern consists of negative residuals (NDOI > USGSmeasured outflow) in winter and spring, positive residuals (NDOI < USGS-measured outflow) in summer and fall, and a zero crossing in July.   period that spans April 1996 through September 2017) when subjected to the baseline NDOI timeseries and an NDOI time series that was altered to reflect the seasonal signal associated with the USGS-measured outflow data. (This is the alternate outflow time series-developed using the first method as introduced in the "Alternate Outflow Time Series Analyses" section.) Model performance is presented as box-and-whisker plots for nine data subsets that follow the outflow and antecedent salinity bins defined in Table 1. Preferred model performance is characterized by a box centered near zero with narrow vertical range. In general, subjecting the empirical X2 models to the alternate NDOI time-series resulted in uniformly poorer performance across all models. However, the alternate NDOI time-series improved model performance in a few limited instances; see performance of the Jassby et al. (1995) and Monismith et al. (2002) models for the three mid-range outflow bins (5,000 cfs < Q < 30,000 cfs). Figure 9 presents box-and-whisker plots of recalibrated empirical X2 model residuals by Figure 8 This box-and-whisker chart compares model residuals from five recalibrated empirical X2 models (for 9 data bins spanning Apr 1996 through Sep 2017) when subjected to the baseline NDOI time series and an NDOI time series that was altered to reflect the seasonal signal associated with USGS-measured outflow.

Seasonal Bias Revealed by Observed X2 Data
month for the three most common (and least dynamic) hydrologic bins defined in Table 1. A pattern of seasonal bias that is consistent with the flow STL analysis is readily observed across all models for the low-range (Q < 5,000 cfs) and mid-range (5,000 cfs < Q < 30,000 cfs) flow bins.  Figure 10 shows data-driven seasonal flow adjustments by month for each empirical X2 model. The monthly estimates are presented as adjustments to gross channel depletion (GCD), where an increase in GCD is equivalent to a decrease in NDOI. The left panel presents the monthly adjustment in flow units of cfs, and the right panel presents the same value expressed as a percent of the overall monthly average NDOI. The plots show non-parametric density estimates of the posterior distribution (as estimated by the MCMC samples) for each combination of month and model. The "overall average" category pools all MCMC posterior samples together and treats that as one overall average distribution. Table 2 shows the mean of each distribution in tabular form; we used these mean values to generate an alternate outflow time-series (using the second method as introduced in the "Alternate Outflow Time Series Analyses" section).

Seasonal Bias Revealed by Data-Driven Flow Adjustments
The empirical X2 models used in this work are based on different assumptions and respond uniquely to flow variation (e.g., see the discussion of dynamic response in Rath et al. [2021, this volume]). Despite these differences, the models provide remarkably similar results for the direction (and, to a lesser degree, magnitude) of seasonal flow adjustments to optimally match salinity observations. The models generally prescribe the least change in flow during the spring months of April through June; they prescribe the greatest relative change in flow (a reduction in GCD or an increase in NDOI) during the summer-fall months of August and September. The models prescribe positive GCD adjustments during the winter months of Figure 10 This diagram shows posterior estimates of seasonal flow adjustments by month for each empirical X2 model. The monthly estimates are presented as adjustments to gross channel depletion (GCD), where an increase in GCD is equivalent to a decrease in NDOI. The left panel presents the monthly adjustment in flow units of cfs and the right panel presents the same value expressed as a percent of the overall monthly average NDOI. VOLUME 19,ISSUE 4,ARTICLE 4 November through February; while the winter adjustments are proportionally small relative to NDOI, they represent significant increases to GCD estimates, ranging from about 23% to 57% in November to about 76% to 128% in January, depending on the empirical model. Figure 11, through re-examination of the box-and-whisker analysis shown in Figure 9, provides a side-by-side comparison of empirical X2 model residuals with and without the datadriven outflow adjustments. In this figure, box widths are proportional to the square root of the number of data points. Thus, very narrow box widths indicate low sample sizes. The figure demonstrates that the adjusted flows generally produce superior model predictions with residuals nearly uniformly closer to being centered at zero. This result is not entirely surprising, given the additional degrees of freedom involved in fitting the flow adjustments. The most dramatic improvement in model performance is concentrated in the late summer and fall months. The seasonal outflow adjustments have a more muted effect on the subset of data represented by high outflow and low X2 conditions, where a very large flow adjustment is needed to effect changes in X2.

DISCUSSION
This work employed a 50-year observed record of salinity in the estuary along with gauged tidal flows and an ensemble of five empirical flowsalinity (X2) models to test whether a seasonal bias in Delta outflow estimates could be inferred. Over this time-period, there are minimal trends in the flow-salinity model residuals, i.e. the models appear to work well across the entire period (Rath et al. 2021, this volume). A specific hypothesis available at the inception of this work, based on limited field studies and island-oriented models (e.g., Owen and Nance 1962;Liang and Suits 2018a), is that NCD (a component of the NDOI water balance) does not equal DICU on daily or monthly time-scales, because soil moisture fluctuates throughout the year-depleting in dry months to provide a time-lagged irrigation supply and accreting during wet months. Thus, on any given island, peaks in consumptive use may not coincide with peaks in net channel depletion. This observation is somewhat supported by differences between appropriately filtered observations from USGS flow gages and NDOI, with measured outflow being greater than the water balance estimate in summer and fall months, and measured flows being less in winter and spring months. It is broadly understood that trends in gauged flow data can be explained by other sources of error (e.g., wind, tides) as described in CDWR (2016). Although empirical X2 model predictions were generally poorer when subjected to outflow adjustments that reflected these seasonal trends, the pattern is sufficiently striking in its magnitude and regularity to warrant further testing; this is currently in progress with a wellresolved 3-D circulation model. To evaluate the lagging hypothesis based on salinity data, the published empirical models were recalibrated with outflow and X2 data sets over a common period (WYs 2000(WYs -2009, and the recalibrated models were run with a seasonally adjusted NDOI to see if the adjustment resulted in improved salinity predictions. The resulting seasonal NDOI flow adjustments, while not fully explaining differences between observed and predicted X2, do indicate a consistent pattern across the models used. Assuming that deviations are associated with channel depletions, the seasonal flow adjustments suggest that a positive adjustment to GCD is needed in the winter months of December through February (i.e., NDOI estimates are biased high), and a negative adjustment to GCD is needed in the late summer months of August and September (i.e., NDOI estimates are biased low). The magnitude of the bias is about 1,000 cfs, which is a small percentage of typical winter outflow but a significant percentage of typical summer outflow. Although

Figure 11
This box-and-whisker chart compares empirical X2 model residuals with and without data-driven outflow adjustments (for three data bins spanning WYs 1968 through 2017). Box widths are proportional to the square root of the number of data points. this signal is weak relative to the errors in the quantification of flows, the finding is notable in that it is obtained from an assessment of five independently developed and independently recalibrated models using 5 decades of continuous salinity measurements.
Furthermore, this finding is consistent with the physical explanation, noted by Owens and Nance (1962), that channel depletions are not equivalent to consumptive use on daily or monthly time-scales because soil moisture fluctuates throughout the year to meet plant water needs. This finding is also consistent with an assertion in CDWR (2016) that, based on historical salinity and modeled outflow estimates, actual outflow may be somewhat higher during the irrigation season and lower during winter months. Importantly, this finding suggests a subtle timing difference that currently available NCD estimates do not include.
The biases reported here are also qualitatively consistent with the DCD model (Liang and Suits 2018a), the development of which has been influenced by salinity results obtained from a 1-D Delta hydrodynamic model (CDWR 2020) over decades of historical modeling. In particular, the DCD model features reduced depletions in summer and increased depletion in winter relative to NCD, although that adjustment is not as focused on particular months as it is here. In DCD, the winter channel depletions are explained as water table recovery, but the water table fluctuation is not explicitly modeled.
The inversion of salinity to seasonal bias outlined here is more dramatic than the associated improvement in the forward application of the empirical X2 models for salinity estimation. If outflow and salinity were perfectly related through the models, the open-ended NDOI  adjustments would be expected to improve salinity predictions up to the time resolution considered and the numerical tractability of the inversion. Instead, there are limits to what might be expected from the seasonal flow adjustments, particularly given that the models perform reasonably well in describing salinity with baseline flows. Model skill is generally limited by errors in salinity observations, X2 interpolations, salinity gate operations, and interactions with tides discussed by Rath et al. (2021, this volume).
All the empirical X2 models considered here are subject to similar limitations in their formulations. Chief among these is that they contain little treatment of the vastly heterogeneous zones in which the X2 isohaline typically resides: Carquinez Strait, Suisun Bay, the confluence of the Sacramento and San Joaquin rivers, and upstream channels are different in shape and depth, and therefore in the dispersion processes they are likely to sustain. Model performance is also limited by the exclusion of important salinity drivers (albeit secondary to Delta outflow) such as wind and tides that have characteristic seasonal cycles. Biases in minor Delta inflow accounting (e.g., Yolo Bypass inflow) are believed to be seasonal, and typical winter operation of the Suisun Marsh Salinity Control structure introduces some seasonal bias into the localized flow-salinity relationship.
We found that seasonal bias associated with the Rath et al. (2017) empirical X2 model-a model that accounts for spring-neap tidal variabilitywas less pronounced, particularly at fortnightly scales, but generally consistent with the other models at seasonal scales. This finding suggests the tidal driver is not responsible for the seasonal bias observed in the other five empirical models. We did not examine seasonal effects of wind on salinity predictions. Bever et al. (2018) identified statistically significant declines in wind speeds in northern San Francisco Bay and the Delta between 1995 and 2015, particularly from October through January. They used numerical modeling to argue that this trend may be related to observed decreases in turbidity over the same time-period. However, Bever et al. did not conduct a similar model analysis of the effect of wind on salinity; there are considerable hurdles to doing so beyond the level of a sensitivity study: (1) a multi-decadal scale analysis of wind effects would be required to separate its effect on salinity from confounding variables that share the same summer seasonality, and (2) the availability of high-quality wind records is limited before the year 2000.
Our finding of a seasonal bias in the reported NDOI over a robust 50-year salinity record builds on and quantifies current physical understanding of the island-channel interactions in the Delta and is especially meaningful during low flow periods. But we can only infer the bias estimates we present here; they cannot be directly linked to specific hydrologic processes on the Delta islands or even to outflow. Thus, this finding adds urgency to field and modeling studies that more directly focus on quantifying the timing and magnitude of Delta NCD. Resulting advancements in scientific understanding, including improved NCD estimates, would lead to fulfilling the state's regulatory and operational goals of developing more accurate NDOI estimates.
As described earlier, NDOI is derived from a water balance that includes several measured flow terms, each of which is associated with some uncertainty. Overestimation of surface water inflows that reach the Delta in the summer is a category of error that is so easy to identify that care must be taken to avoid a "bias of discovery" in NDOI accounting. Correcting these errors would reduce NDOI estimates in the summer. By contrast, evidence presented here based on measured outflow and salinity suggests just the opposite: that summer outflow is actually higher than current NDOI estimates. The dilemma is that-assuming the bias is related to NCD estimates-errors are related to subsurface flow estimates (or gauging error) that are not as directly observable or as easy to quantify as an omitted source, dry river bed, or ditch would be for surface-water components.
Future work could apply our methodology to explore the possible seasonal bias in another flow term that is of considerable importance in Delta water management: combined Old and Middle River (OMR) flow. Water balance estimates of this flow (Andrews et al. 2016) are expected to be similarly influenced by the timing of NCD, as is the Delta outflow estimate through the NDOI water balance. However, OMR flow is significantly influenced by the seasonal placement of temporary barriers in certain South Delta channels that may confound the relationship with island hydrologic exchanges. Despite this concern, the finding of a similar seasonal bias in OMR flow estimates would provide greater support for the concept identified here: that the full Delta island water balance significantly influences the timing of island water exchanges with Delta channels, and that evapotranspiration and precipitation estimates alone do not adequately represent it on daily and monthly time-scales.
This work addressed uncertainty in net Delta channel depletion estimates indirectly by means of tidal flow data, salinity data, and flow-salinity relationships aggregated at a large scale. The findings, however, point to the necessity of collecting more local-scale data at the island level to better quantify the water balance on individual islands, and to directly test the hypothesis of water storage on islands as a mechanism that modifies the timing of channel depletions in the Delta. Indeed, such measurement can only be done rigorously at the island scale and is inherently inferential at the scale of the entire Delta.
Current work has focused on evapotranspiration measurements at a fine spatial resolution across the Delta (Medellín-Azuara et al. 2018), although measurements of island diversions-required through the implementation of California's Senate Bill 88-have been limited in extent. Such island-scale measurements-coupled with (1) measurement of other island water balance terms (including base flows from the Sierra Foothills region on the eastern boundary of the Delta), (2) updates to Delta-wide consumptive use patterns associated with land use changes since the 1960s, and (3) ongoing refinements to Delta-wide modeling-will go a long way toward addressing the uncertainties of estimating Delta outflow and salinity during critical low-flow periods, when careful management of available water resources in California is most needed.