Advanced searchCitation search
Skip main navigation
Close Drawer MenuOpen Drawer Menu
Previous article
Next article
Open access
Technical Papers
Dec 18, 2015

Nine Decades of Salinity Observations in the San Francisco Bay and Delta: Modeling and Trend Evaluations

Authors: Paul H. Hutton, John S. Rath, Limin Chen, Michael J. Ungs, and Sujoy B. Roy, A.M.ASCE [email protected]Author Affiliations
Publication: Journal of Water Resources Planning and Management
Volume 142, Issue 3
https://doi.org/10.1061/(ASCE)WR.1943-5452.0000617
PDF

Abstract

The position of the low salinity zone in the San Francisco Bay Delta—given its correlation with the abundance of several estuarine species—is used for water management in a system that supplies water to more than 20 million people and contains one of the most diverse ecosystems on the Pacific coast. This work consolidates legacy and modern salinity data to develop a reasonably complete daily record spanning nine decades. The position of the low salinity zone, which is effectively characterized by an empirical model that was developed to support data cleaning and filling, reveals statistically significant trends consistent with increasing water demands and introduction of upstream reservoirs, e.g., increasing salinity trends in wet months and decreasing salinity trends in dry months. Reservoir effects are particularly apparent in drier years, with greater seasonal variability in the early part of the record before major reservoirs operated in the watershed. These data provide a basis for further analysis of how and why the position of the estuary’s low salinity zone has changed over time.

Introduction

Freshwater inflows have a direct influence on the salinity structure in estuaries. In the San Francisco Bay, the salinity structure has been related to the health of estuarine species in the Suisun Bay and the western Delta (Fig. 1). In particular, the location or position of two parts per thousand (ppt) bottom salinity—hereafter referred to as X2—has been correlated with the abundance of several species (Jassby et al. 1995). Using data collected over different time periods, the low salinity zone in general and the X2 position in particular have been associated with the greatest abundance of pelagic organisms, including the protected longfin smelt (Spirinchus thaleichthys) (Jassby et al. 1995) and the Delta smelt (Hypomesus transpacificus) (Feyrer et al. 2011). The X2 position has also been associated with the abundance of undesirable species such as the invasive Asian clam (Corbula amurensis). The relationship between the low salinity zone and the responses of individual species are a topic of continued research interest (Feyrer et al. 2007, 2011; Kimmerer et al. 2009; Moyle et al. 2010), and the broader science underlying the driving mechanisms between water quality, habitat quality, and species abundance continues to evolve (Reed et al. 2014).
Fig. 1. Key salinity stations are identified in Suisun Bay and the western Delta; salinity data from these and other locations were used to develop a long-term record of X2, the position of two parts per thousand bottom salinity in the estuary; the X2 position is reported as the distance in kilometers from Golden Gate along the axis of the estuary, following the original definition of the term (Jassby et al. 1995); the Sacramento and San Joaquin River branches are identified on the map
The position of the X2 isohaline (defined as the distance from Golden Gate in kilometers, Fig. 1) during the months of February through June is currently used as the basis of flow management in San Francisco Bay (CSWRCB 1999). Estuarine flows can be managed through upstream reservoir releases and exports of water from the Delta. The recent biological opinion on Delta smelt (USFWS 2008) regulates X2 position in fall months (September through November) following wet and above normal water years. Much of the published literature on X2 and its relationship to various biological indicators is based on data collected over limited periods, typically spanning the mid-1960s to the present.
Although X2 is defined in terms of bottom salinity, much of the published analysis is based on surface salinity measurements, including the seminal work on X2 (Jassby et al. 1995). Use of surface salinity as a surrogate for bottom salinity is largely motivated by the abundance of surface salinity measurements throughout the estuary, due in part to historical precedent and to the operational challenges of maintaining salinity sensors at depth. However, the estuary is known to be vertically stratified, with increasing stratification at greater river flows (Monismith et al. 2002). Stratification has been addressed by using a constant factor to relate the bottom salinity to surface salinity, i.e., two ppt bottom salinity is assumed to correspond to 1.76 ppt surface salinity (Jassby et al. 1995). Current regulations assume two ppt bottom salinity corresponds to 2.64  mS/cm surface specific conductance (CSWRCB 1999).
Given the importance of the low salinity zone for estuarine species, and of X2 in the management of water in the estuary, the present analysis builds on past work by extending the readily available surface salinity data. The earliest salinity data incorporated in this work are based on technical reports published by the California Department of Public Works (CDPW) and its successor agency, the Department of Water Resources (CDWR), beginning in the 1920s. This work also extends previously published salinity trend evaluations in the Bay-Delta, which have focused on more limited time periods or on station-specific salinity rather than isohaline position (Fox et al. 1991; Shellenbarger and Schoellhamer 2011; Enright and Culberson 2009; Moyle et al. 2010). Although the data used here do not represent predevelopment conditions such as those obtained through analysis of paleoclimatic signals (Stahle et al. 2001), they do represent a wide range of hydrologic conditions and watershed development activities, including reservoir construction, water exports, and land use changes (Fox et al. 1990).
The major objectives of this work were to (1) develop a cleaned database for salinity across Suisun Bay and the western Delta for the longest observational record possible and compute isohaline positions at each point in time, (2) develop and calibrate an empirical salinity model that could be used for additional diagnostic evaluation of the data, and (3) evaluate changes in the isohaline position over the 91-year period of record from water year 1922 to 2012. Water years in California begin on October 1 of the preceding calendar year. By extending the starting point from the mid-1960s to the early 1920s, the readily available data set now incorporates a period of record prior to the construction of major water storage and diversion projects (completed between 1944 and 1967) and a period of severe drought that occurred between 1928 and 1934.
Not surprising given the extensive period of record, the data compiled in this work are not noise-free and error-free, arising in part from variations in sampling and analytical methodology. A significant effort was expended to clean the data to remove values that appeared to be inconsistent with other values. These data were then used to develop daily salinity estimates at each station, and gaps were filled through interpolation and comparison with neighboring stations. The empirical salinity model, based on a formulation accounting only for flow inputs, was calibrated using these cleaned data. Finally, statistical analyses were performed on the individual station salinity and interpolated isohaline positions to detect changes over time and across different water year classes.

Methods

Salinity Data Sources and Cleaning

The data incorporated in this work include historical grab sample data and modern conductivity sensor data. The historical grab sample data record, hereafter referred to as the Bulletin 23 data record, is based on legacy reports spanning the period October 1921 to June 1971 (CDWR 2015b), and also summarized in Roy et al. (2014). Scanned paper copies of these reports were used to develop an electronic database of salinity throughout the Delta and portions of San Francisco Bay. An important salinity data set that predates the Bulletin 23 data (but was not employed in this study) is based on records by the California Hawaiian Sugar Refining Corporation (C&H). C&H, which obtained most of its fresh water supply in the early twentieth century by transporting water to its refinery in Crockett, maintained a record on the distance its barges traveled to obtain fresh water (typically less than 50  mg/l chloride) and the quality of water obtained (CDWR 2015b; Lund et al. 2007). Although the C&H records are of great historical interest and demonstrate the seasonal variability in the salinity field prior to extensive upstream development, the nominal isohaline position of 50  mg/L chloride was not reported with commensurate tidal cycle information and was too low (i.e., too fresh) to accurately characterize the general salinity gradient.
Modern databases were queried from several sources to supplement the Bulletin 23 data, including: (1) the California Data Exchange Center (CDEC 2013), (2) the interagency ecological program, and (3) U.S. EPA’s STOrage and RETrieval (STORET). These modern data, hereafter referred to collectively as CDEC data, were further supplemented by U.S. Geological Survey data (USGS 2013) to represent high outflow periods when the low salinity zone extended far downstream into the San Francisco Bay. The combined data gathering effort resulted in a master database containing salinity records from October 1921 to September 2012, i.e., water years 1922–2012. The locations of key salinity stations used in this analysis are shown in Fig. 1. Additional information on data sources, stations, and time periods are provided in the “Supplemental Data” section, Tables S1–S3.
The raw data contained errors associated in part with variations in sampling and analytical methodology. The Bulletin 23 data report salinity as chloride concentrations. The CDEC data report salinity as specific conductance, or electrical conductivity (EC) standardized to 25°C. All data were converted to specific conductance in units of milliSiemens (mS/cm) using regression relationships developed from colocated chloride and specific conductance data in the estuary (Denton 2015).
The CDEC data are collected by continuous EC sensors that report every 15 min, and daily averages were computed directly using these subdaily values. Given that the original subdaily data were frequently unavailable, averages were computed over 24 h rather than a tidal day (25 h). Monismith et al. (2002) reported that the errors associated with this approximation were very slight. The Bulletin 23 data were collected nominally every four days at higher high tide or low high tide. Because estuarine salinity can vary significantly over the course of a day, these grab sample data were converted to approximately equivalent daily averages using simulation output from a hydrodynamic and water quality transport model, Delta Simulation Model version 2 (DSM2), a linked-node model that is widely used for studying Delta flow, stage, and water quality (CDWR 2015a). This tidal correction was successfully validated by comparing the resulting daily average estimates with co-occurring CDEC data (Roy et al. 2014). Enright and Culberson (2009), when confronted with the same problem, tidally corrected Bulletin 23 grab sample chloride data through linear correlations with co-occurring CDEC specific conductance data to produce long-term salinity time series records for three stations in the estuary.
Additional data cleaning and filling was performed by comparing daily average specific conductance at pairs of stations and assuming that under conditions of moderately high salinity reflecting strong ocean influence, salinity decreases monotonically downstream to upstream. When data at a pair of stations are inconsistent with this behavior, i.e., an eastern (upstream) station has a higher salinity than a western (downstream) station, a procedure was required to determine which of the two salinity values was erroneous, acknowledging there is no a priori way to make this determination. To perform the data cleaning step, data were correlated with nearby station data using least-squares regressions. Measured values that differed greatly (by more than four standard errors) or too frequently (by more than two standard errors multiple times) from regression predictions were removed from the data set. This step is considered an approximate way to remove potentially erroneous values from the data set, and it is possible that some true data values are excluded in the process. However, because this analysis is not focused on the behavior of extreme values, this approach is unlikely to affect the conclusions.
The method used to calculate isohaline position, discussed in the next section, requires a reasonably complete salinity record. Missing values were filled based on the salinity data of nearby stations using the correlations discussed in the previous paragraph. Filling missing downstream station values from upstream station data was found to be particularly challenging when upstream conditions were fresh, as downstream salinity can vary across orders of magnitude for the same low (or fresh) upstream salinity. After this step was completed, any remaining short gaps (up to eight days) were filled through linear interpolation. When there was an overlap of the Bulletin 23 and CDEC data (i.e., the 1964 to 1971 period), the latter were used in preference.

Isohaline Calculations

Isohaline position was calculated through interpolation of the cleaned and filled salinity record. Theoretically, different interpolation approaches may be used to calculate X2 position. The longitudinal salinity gradient changes with flow and with distance along the estuary (among other factors); thus the estimated isohaline position is somewhat dependent on the interpolation approach and stations used. The daily X2 position was estimated assuming a log-linear relationship between surface salinity and distance, interpolating across two stations that bound a specific conductance of 2.64  mS/cm (which under current regulations is assumed to correspond to a bottom salinity of 2 ppt). In a limited number of cases a weighting approach over additional stations was used if the data exhibited nonmonotonic behavior near a salinity value of interest. If the bounding stations were further apart than 25 km on any given day, the X2 position was not estimated because of uncertainty about interpolation accuracy. This condition resulted in X2 position not being estimated for 3.2% of the days over the study period that had one or more salinity data points. This interpolation method was used to calculate unique isohaline positions along the Sacramento and San Joaquin River branches upstream of their confluence (Fig. 1).
Monthly X2 values were estimated from the daily interpolated isohaline values. Monthly X2 position was defined as the mean value of all nonmissing daily X2 values for months where atleast 14 daily values were computed. Using similar methods, additional surface salinity isohalines (e.g., six ppt surface salinity isohaline) were estimated on daily and monthly time steps to more fully characterize the estuary’s low salinity zone.

Modeling Approach

Denton (1993) developed an approach to estimate salinity at fixed locations in the estuary based on a modification of the steady-state solution of the tidally-averaged advection-dispersion equation for salinity transport in a one-dimensional estuary. His empirical approach utilizes boundary conditions representative of the downstream ocean and upstream riverine environments and a concept called antecedent outflow, representing flow time-history in the estuary. The equation can be represented as
S(t)=(So−Sb)×exp[−α×G(t)]+Sb
(1)
where S(t) = salinity at a given location; So and Sb = downstream (i.e., ocean) and upstream (i.e., riverine) salinity boundaries, respectively; α = empirically determined location-specific constant (units  of  flow−1); and G(t) = measure of the antecedent outflow. Antecedent outflow is defined by the following routing function, similar to one proposed by Harder (1977):
∂G∂t=[Q(t)−G(t)]×G(t)β
(2)
where Q = Delta outflow; and β = empirically determined constant (units  of  flow·time). Denton (1993) observed that the term β/G is a time constant governing the rate at which G approaches steady state. These equations can be calibrated to predict site-specific salinity.
In reference to an autoregressive empirical model for calculating the X2 position proposed by Jassby et al. (1995), Monismith et al. (2002) argues on theoretical grounds that power-law relationships with flow are preferable over logarithmic relationships and proposed an autoregressive X2 function of the following form:
X2(t)=ω1×Q(t)ω2+ω3×X2(t−1)
(3)
where ω1, ω2, and ω3 = empirically determined constants.
Jassby et al. (1995) observed that the entire mean salinity field can be predicted if the X2 position is known, i.e., the salinity field is self-similar and can be predicted as a function of the longitudinal distance from Golden Gate (X) when normalized by X2. Thus, salinity as a function of X/X2 is relatively uniform for a wide range of flows. Following this observation, an integration of the Eulerian modeling approach of Denton (1993)—focused on a fixed station—and the Lagrangian modeling approach of Monismith et al. (2002)—focused on a fixed salinity—was performed to develop a tool for diagnostic applications in the salinity data cleaning and filling process. The resulting empirical model, which is capable of estimating salinity at variable locations and X2 and other isohaline positions, is termed the Delta Salinity Gradient (DSG) model. Formulation of the DSG model is described briefly in the remainder of this section. Details on model formulation are provided elsewhere (Hutton 2014).
The steady state solution to Eq. (3) can be derived by setting X2(t)=X2(t−1) to obtain X2¯=ω1/(1−ω3)×Q¯ω2, where X2¯ and Q¯ denote steady state conditions. Substituting antecedent flow G(t) for steady state flow Q¯ gives an approximation to the unsteady response of X2 to flow variations if G(t) does not vary too rapidly. This substitution of antecedent flow is similar in concept to, and motivated by, Denton’s (1993) derivation of the empirical Eq. (1), in which he proposed using the G flow instead of Q¯ in a steady state analytical solution of salinity transport. Reparametrizing the constants as Φ1=ω1/(1−ω3) and Φ2=ω2 gives a new empirical relationship between X2(t) and G(t) where Φ1 and Φ2 are independently calibrated to X2 from observed data
X2(t)=Φ1×G(t)Φ2
(4)
This empirical formulation, in contrast to those proposed by Jassby et al. (1995) and Monismith et al. (2002), is capable of estimating X2 during the early period of record when daily (and even monthly) Delta outflows frequently turned negative. Redefining the location-specific constant α as a function of X and scaling distance to the X2 isohaline (S=2.64  mS/cm) results in the following relationship:
S(X,t)=(So−Sb)×exp{τ×[XX2(t)]−(1/Φ2)}+Sb
(5)
where τ=ln[(2.64−Sb)/(So−Sb)]; and salinity is reported as specific conductance in units of mS/cm. Eq. (5) implicitly assumes that the estuary’s salinity structure is self-similar under all flow conditions. However, Monismith et al. (2002) showed that the structure changes under high flow conditions. To address this response to flow, the downstream boundary condition So is assumed to vary with X2 as a sigmoidal function:
So(t)=S^+(2.64−S^)×exp[−γ×X2(t)δ]
(6)
where S^ = ocean salinity (≈53  mS/cm); and γ and δ = empirically determined constants. Eq. (5) can be used to determine salinity at any longitudinal distance from Golden Gate given X2 position and Φ2 and assuming a reasonable value for Sb. If appropriate salinity observations are unavailable, X2 can be estimated from antecedent outflow using Eq. (4). Eq. (5) can be rearranged to predict surface salinity isohaline positions as a function of X2:
X(S,t)=X2(t)×{ln[S−SbSo(t)−Sb]τ}−Φ2
(7)

Statistical Analyses

Sen’s nonparametric estimate of slope (Gilbert 1987 and references therein) was used to perform a trend analysis of the monthly X2 estimates over the entire period of record and two additional intervals: water years 1922 to 1967 and 1968 to 2012. These intervals were selected to coincide with Enright and Culberson’s (2010) pre–water project periods and post–water project periods. The significance of the breakpoint between periods is that, although the Central Valley and State Water Projects began pumping water from the Delta in 1940 and 1967, respectively, they did not begin year-round pumping operations until 1968 when the San Luis Reservoir was completed to store water south of the Delta. The Sen slope is the median of all slopes between all possible unique pairs of individual data points in the time period being analyzed. If there are n time points or periods of time, then there are a total of n(n−1)/2 possible pairs of time points one could use to calculate a slope, and Sen’s slope is the median of these values. The method is robust and fairly insensitive to the presence of a small fraction of outliers, nondetect, or extreme data values; thus, trend estimates based on Sen slope are not biased by the occurrence of drought in the early part of the record.
The Mann-Kendall test was performed on the Sen slope at the 95% confidence level (Gilbert 1987 and references therein), with results identified as either an upward trend, a downward trend or no trend. The trend slope was computed using the median value of the Sen slope. Nonzero slopes may or may not be found to be statistically significant using the Mann-Kendall test. The nonparametric Wilcoxon Rank Sum test was used for the comparison of isohaline values for specific water year classes (i.e., wet, above normal, below normal, dry, and critically dry). The Mann-Kendall test was also performed on monthly average specific conductance values over the entire period of record at five locations. Details of the implementation of the statistical procedures are presented in the “Supplemental Data” section, Appendix S1.

Results

Cleaned and Filled Salinity Data

Summary statistics for the resulting cleaned and filled daily average surface specific conductance data based on the Bulletin 23 grab samples are shown in Table 1. Although the filling process provides a fairly complete record for key stations downstream of approximately 100 km, substantial gaps remain in upstream station records that were used exclusively to characterize extreme drought conditions in the 1920s and 1930s. Statistics for these salinity stations are not provided in Table 1. Similar statistics for the CDEC data are presented in Table 2. The cleaned and filled CDEC data show a more complete record than the Bulletin 23 data across all stations. These data are provided electronically in the Supplemental Data. Given our goal to interpolate X2 position and other isohaline positions in the low salinity zone, the available data provide an adequate basis for the calculation.
Table 1. Bulletin 23 Data Summary
Station nameDistance from golden gate (km)Approximate period of recordCleaned and filled data completenessSpecific conductance percentiles (mS/cm)
Count (days)Missing (%)10%50%90%
Bay stations
 Point Orient19.8February 1926–June 197113,97616202729
 Point Davis40.6February 1926–June 197114,802117.02027
 Crockett44.6February 1926–June 197113,685175.71926
 Benicia52.3February 1926–June 197113,706173.61624
 Martinez52.6February 1926–June 197113,760171.91322
 Bulls Head Point54.7February 1926–August 19579,007221.91525
 West Suisun59.5October 1921–June 197113,410260.79.721
 Bay Point64.2October 1921–December 196811,174350.59.322
 Port Chicago66.0October 1921–June 197114,745190.38.321
 O and A Ferry74.8October 1921–June 197115,522140.31.714
Lower Sacramento River stations
 Collinsville81.8October 1921–June 197115,751130.20.48.9
 Emmaton92.9October 1921–June 197115,185160.20.32.0
 Three Mile Slough Bridge96.6October 1921–June 197115,178160.20.30.8
 Rio Vista102.2September 1922–June 197114,408190.10.20.3
Lower San Joaquin River stations
 Antioch88.4October 1921–June 197115,451150.20.46.0
 Antioch Bridge93.7October 1921–June 197114,760190.30.31.8
 Jersey Point98.8October 1921–June 197115,380150.20.31.6
 False River101.2October 1921–June 197113,883240.20.20.8
 Oulton Point108.1September 1952–June 19715,395210.10.20.4
 San Andreas Landing113.1September 1952–June 19715,395210.10.20.3

Note: Statistics for the resulting cleaned and filled daily average Bulletin 23 specific conductance data are shown for key locations by river branch.

Table 2. CDEC Data Summary
Station nameDistance from golden gate (km)Approximate period of recordCleaned and filled data completenessSpecific conductance range (mS/cm)
Count (days)Missing (%)10%50%90%
Bay stations
 Point San Pablo22January 1965–September 201216,8393253844
 Carquinez45.5January 1965–September 201217,0102122736
 Martinez54September 1995–September 20126,03332.71826
 Martinez (USBR)55January 1965–April 199611,34512.21627
 Port Chicago64January 1965–September 201217,38900.39.320
 Mallard Island75July 1964–September 201217,50510.23.012
Lower Sacramento River stations
 Collinsville81July 1964–September 201216,98540.21.07.7
 Emmaton92July 1964–September 201217,42010.10.22.2
 Rio Vista101July 1964–September 201217,42010.10.20.3
Lower San Joaquin River stations
 Pittsburg77January 1965–September 201217,40500.22.310
 Antioch85.8July 1964–September 201217,56100.20.74.8
 Blind Point92.9July 1964–September 201217,54000.20.42.4
 Jersey Point95.8July 1964–September 201217,38810.20.31.7
 Three Mile Slough at San Joaquin River100.4July 1964–September 201217,32020.10.31.1
 San Andreas Landing109.2July 1964–September 201217,52600.10.20.3

Note: Statistics for the resulting cleaned and filled daily average CDEC specific conductance data are shown for key locations by river branch. CDEC = California Data Exchange Center; USBR = U.S. Bureau of Reclamation.

Interpolated and Model-Predicted X2 Position

Daily and monthly X2 positions were estimated for the period of record using the previously described approach. Daily X2 position was also predicted from the DSG model for the same period following the succeeding procedure:
•
Antecedent outflow was calculated from Eq. (2) assuming a nominal value for β of 475 cfs-years and assuming daily Delta outflows (Q) from the DAYFLOW model (CDWR 2015b). As detailed elsewhere (Hutton 2014), daily outflows prior to October 1929 were estimated from monthly outflow volumes and daily inflow volumes (CDWR 2015b). The same Delta flow time series is used for calibrating the model for both the Sacramento and San Joaquin River branches, and the channel-specific responses are embedded in the fitted model parameters.
•
Interpolated daily X2 values for the Sacramento River branch, spanning Water Years 2000 through 2009, were used to calibrate model parameters Φ1 and Φ2 from Eq. (4) through least-squares minimization. Best fit parameter values Φ1=456±3.93 (mean±1  SE) and Φ2=−0.193±0.001 resulted after data points representing extremely high outflow events (X2<38  km) were removed from the analysis. The coefficient of determination R2=0.92 and the standard error of estimate is 3.2 km. Our parameter estimates, when estimated using antecedent outflow in comparable units (m3/s), are similar to those reported by Gross et al. (2009) for various steady fit models. Differences in parameter estimates are attributed primarily to the use of a different calibration period. The autoregressive X2 function proposed by Monismith et al. (2002) [Eq. (3)] was calibrated to the same data set using a nonlinear least-squares fitting procedure and resulted in a coefficient of determination R2=0.89 and a 3.8 km standard error of estimate, when applied with modeled data, i.e., at each timestep, the antecedent X2 in the equation was based on the modeled value.
•
Best fit parameter values were also calibrated for the San Joaquin River branch, resulting in Φ1=502±4.63 and Φ2=−0.203±0.001 with R2=0.92 and a 3.6 km standard error of estimate. X2 values along the San Joaquin River branch are typically higher (i.e., further upstream) than those along the Sacramento River branch, due in large part to smaller freshwater inflow volumes from the San Joaquin River available to repel salinity. Eq. (3) was calibrated to the same data set, as reported previously, and resulted in a coefficient of determination R2=0.89 and a 4.1 km standard error of estimate. For both river channels, therefore, the DSG fits are slightly better than that obtained from Eq. (3).
A time series of the daily X2 position along the Sacramento River branch is shown in Figs. 2(a–d) over the full 91-year period of record. The time series reveals a wide range in daily X2 position from approximately 20 km to greater than 100 km. At the lower extreme, the X2 falls in a broad region of the estuary (San Pablo Bay), where the one-dimensional approach may be limiting and there may be significant lateral gradients in salinity. The X2 position is generally more upstream (i.e., higher) in dry and critically dry years, corresponding to sustained periods of low Delta outflow. The trace in Fig. 2(a), representing a period before Shasta Dam and other large upstream reservoirs were constructed, is visually distinct from the remaining time series. X2 values calculated from Eq. (4) are superimposed on the interpolated X2 values for comparison. The DSG model fits the time series reasonably well, with some exceptions in the pre-Shasta period corresponding to extreme salinity incursions during major drought periods that were well beyond the model’s calibration range. The generally slow rate of change in the salinity field and the use of an antecedent outflow term appear to justify the steady state approximation under most nonextreme flow conditions. The DSG model predictions show some seasonal bias when compared with interpolated X2 values (Roy et al. 2014). This bias is hypothesized to be related to inaccuracies associated with estimating net water use by agriculture in the Delta, particularly during low flow periods when this water use is a significant fraction of the Delta outflow water balance.
Fig. 2. Time series of interpolated and DSG-predicted daily X2 values on the Sacramento River branch (water years 1922–2012); the X2 position is generally more upstream (i.e., higher) in dry and critically dry years, corresponding to sustained periods of low Delta outflow; (a) represents a period prior to large reservoir construction in the Sacramento Valley, is visually distinct from the remaining time series; (a and b) represent the preproject period; (c and d) represent the postproject period
The interpolated monthly X2 time series was evaluated by grouping individual values into water year classes. Fig. 3 shows the monthly X2 position for the Sacramento River branch averaged by water year class for the previously defined preproject and postproject periods. The difference between preproject X2 and postproject X2 is greatest in critically dry years and diminishes with wetter conditions (plot panels from left to right). The postproject period exhibits a dramatically reduced X2 range, relative to the preproject period, during dry and critically dry water years. This reduced range is characterized by higher values in winter and lower values in summer. Water project operations, which typically store runoff in the winter and spring and release storage in summer months to maintain in-basin water quality standards, clearly have a strong influence on the estuary’s intra-annual salinity pattern except under wet hydrologic conditions. However, the differences between preproject and postproject conditions shown in Fig. 3 cannot be fully attributed to operations of the central valley and state water projects. Intensified upstream agricultural and urban water use and associated water projects for in-basin and out-of-basin water uses, and changes in estuarine geometry, mean sea level and watershed snowmelt patterns have also contributed to changes in salinity patterns. A similar figure for the San Joaquin River branch is shown in the “Supplemental Data” section, Fig. S1.
Fig. 3. Average monthly X2 position is shown by water year class on the Sacramento River branch under preproject (water years 1922–1967) and postproject (water years 1968–2012) conditions, with lines connecting the seasonal medians; symbols show individual year values (red = preproject and blue = postproject), with lines connecting the seasonal medians; in all but wet years, the postproject X2 position tends to be further downstream (i.e., lower) in summer months and further upstream (i.e., higher) in other months; the X2 position in October and November is generally more closely associated with the previous water year; thus, the x-axis spans the months December through November

Other Model Predictions

Eqs. (5) and (6) of the DSG model were applied to predict daily specific conductance at Collinsville over the six-year drought period Water Years 1928–1934 using a subset of the model-predicted X2 time series illustrated in Fig. 2. Collinsville (X=81  km) was selected to illustrate the model’s predictive capability as this station plays a critical role in X2 management during spring and fall. To conduct the illustrative simulation, the following model constants were assumed: Φ2=−0.193; Sb=0.2  mS/cm; γ=2.29×10−4; and δ=1.83. Fig. 4(a) compares the DSG-predicted time series with the cleaned and filled specific conductance data. The time series is also compared with predictions from a site-specific calibration of Eq. (1) reported by Denton and Sullivan (1993). The DSG model effectively represents the observed salinity variation at Collinsville over two orders of magnitude, although the extreme event in 1931 is overpredicted. Furthermore, the DSG model provides salinity estimates comparable to those provided by the site-specific empirical model.
Fig. 4. These figures illustrate predictive capability of the DSG model during a six-year drought in the early part of the record (water years 1928–1934): (a) shows a time series of observed and DSG-predicted salinity [Eqs. (5) and (6)] for a representative station, Collinsville (X=81  km), following the data cleaning and filling procedures described in the text; predictions from a site-specific empirical model (Denton and Sullivan 1993) are provided for comparison; (b) shows a time series of observed and DSG-predicted surface isohalines [Eq. (7)] that bound the estuary’s low salinity zone (1–6 ppt)
To further illustrate the utility of the DSG model, Eq. (7) was applied with the same model constants to predict low salinity zone position (bounded by surface salinities of one to six ppt) for water years 1928–1934. Fig. 4(b) compares the DSG-predicted time series with the interpolated isohaline data. Again, the data provide a reasonable validation of the DSG model except for the extreme event in 1931.

Isohaline Position Trend Analysis

Statistical analyses were performed on the interpolated X2 values to characterize behavior over time and in response to different hydrologic conditions. Results from the trend analysis are shown in Tables 3 and 4 with the analysis focusing on the Sacramento River branch. Similar analyses for the San Joaquin River branch are shown in the “Supplemental Data” section (Tables S4 and S5). Key results for the Sacramento River branch are summarized as follows:
•
The monthly trend evaluation for the entire period of record (1922–2012) shows statistically significant increases in X2 from November through June. Statistically significant decreases in X2 occur in August and September.
•
Over the preproject period (1922–1967), there is no significant change in X2 from January through July and a statistically significant decrease in X2 from August through December. The trend directions are identical for both river branches.
•
Over the postproject period (1968–2012), there is a nearly inverse response in trends, with a statistically significant increase in X2 from September through December. Again, the trend directions are identical for both river branches.
Table 3. Mann-Kendall Test Results: Sacramento River Branch X2
MonthFull period 1922–2012Preproject period 1922–1967Postproject period 1968–2012
Sample sizeSen’s slope median (km/year)Test decisionSample sizeSen’s slope median (km/year)Test decisionSample sizeSen’s slope median (km/year)Test decision
January830.10Up39−0.11No trend440.23No trend
February820.09Up39−0.03No trend430.10No trend
March830.08Up400.09No trend430.02No trend
April820.15Up390.19No trend430.01No trend
May850.13Up400.13No trend45−0.18No trend
June850.10Up400.01No trend45−0.08No trend
July87−0.04No trend42−0.04No trend45−0.06No trend
August86−0.13Down41−0.2Down450.06No trend
September88−0.12Down43−0.43Down450.20Up
October880.00No trend43−0.31Down450.28Up
November860.11Up41−0.2Down450.37Up
December850.13Up40−0.19Down450.37Up

Note: Over the entire period of record, the test shows (1) statistically significant increases in X2 from November through June and (2) statistically significant decreases in X2 for August and September. Over the preproject period (1922–1967), there is no significant change in X2 from January through July and a statistically significant decrease in X2 from August through December. Over the postproject period (1968–2012), there is a nearly inverse response in trends, with a statistically significant increase in X2 from September through December. Results are reported as an upward trend (up), a downward trend (down), or no trend.

Table 4. Wilcoxon Rank Sum Test Results: Sacramento River Branch X2
MonthYear type
CriticalDryBelow normalAbove normalWet
January7.6 (no trend)9.1 (up)6.4 (no trend)−1.1 (no trend)3.7 (no trend)
February16.3 (up)10.4 (up)4.4 (up)0.7 (no trend)5.8 (up)
March14.2 (up)7.4 (up)3.2 (no trend)0.9 (no trend)1.7 (no trend)
April13.9 (up)9.1 (up)9.0 (up)7.2 (no trend)4.5 (no trend)
May12.8 (up)9.1 (up)12.9 (up)2.4 (no trend)6.3 (up)
June2.0 (no trend)5.5 (no trend)15.5 (up)1.0 (no trend)6.8 (up)
July−10.6 (up)−3.2 (down)6.1 (up)−7.3 (down)1.7 (no trend)
August−12.0 (down)−7.2 (down)−2.8 (no trend)−9.8 (down)−3.1 (down)
September−15.6 (down)−3.8 (no trend)−0.8 (no trend)−6.6 (down)−5.4 (down)
October1.9 (no trend)0.6 (no trend)−5.5 (no trend)6.1 (no trend)0.0 (no trend)
November9.4 (up)6.9 (no trend)2.7 (no trend)15.4 (no trend)3.8 (no trend)
December9.2 (up)8.6 (up)4.2 (no trend)11.8 (no trend)3.0 (no trend)

Note: In general, dry and critical year postproject X2 was statistically significantly higher (i.e., upstream) in December through May and lower (i.e., downstream) in August and September. Wet year postproject X2 was statistically significantly higher in May and June and lower in August and September. Results are reported as a nonparametric estimate of the median of the difference (km) between a postproject X2 and a preproject X2, and significance is reported as an upward trend (up), a downward trend (down), or no trend.

The nonparametric Wilcoxon Rank Sum test was used for interperiod comparison of X2 position by month and water year class. The results of the Wilcoxon Rank Sum test, assuming a 95% confidence level, are summarized for the Sacramento River branch in Table 4. In general, postproject X2 positions during dry and critically dry water years were statistically significantly higher (i.e., upstream) in December through May and lower (i.e., downstream) in August and September. Although the test shows fewer statistically significant trends under wetter conditions, the trend of lower August and September X2 during the postproject X2 held. These statistical tests add more detail to the visual patterns displayed in Figs. 2 and 3.
To further evaluate the isohaline trend analysis results, the Mann-Kendall test was performed on observed and DSG-predicted monthly average specific conductance values over the entire period of record at five locations: Martinez, Port Chicago, Mallard Island (represented by the O&A Ferry location in the Bulletin 23 data), Collinsville, and Emmaton. Location-specific trends are compared with X2 trends (both derived from DSG predictions) in Table 5. Although observed data trends generally matched predicted data trends, the latter are presented to avoid bias that may be introduced by gaps in the observed salinity record. When a trend was detected in both the salinity and X2 time series, the trends are uniformly consistent. When a trend was not detected in the X2 time series, the salinity trends are generally consistent, with exceptions in January (Mallard Island and Collinsville) and March (Mallard Island).
Table 5. Mann-Kendall Test Results for DSG-Predicted Salinity and X2 at Selected Locations for the Entire Period of Record (1922–2012)
MonthSalinity trend: 1922–2012X2 Trend: 1922–2012
MartinezPort ChicagoMallard IslandCollinsvilleEmmaton
JanuaryNo trendNo trendUpUpNo trendNo trend
FebruaryUpUpNo trendNo trendNo trendUp
MarchNo trendNo trendUpNo trendNo trendNo trend
AprilUpUpUpUpNo trendUp
MayUpUpUpUpNo trendUp
JuneUpUpUpUpNo trendUp
JulyNo trendNo trendNo trendNo trendNo trendNo trend
AugustDownDownDownDownDownDown
SeptemberDownDownDownDownDownDown
OctoberNo trendNo trendNo trendNo trendNo trendNo trend
NovemberNo trendNo trendNo trendNo trendNo trendNo trend
DecemberUpUpUpUpUpUp

Note: Detected salinity and X2 trends are generally consistent. When a trend was detected in both the salinity and X2 time series, the trends are uniformly consistent. When a trend was not detected in the X2 time series, the salinity trends are generally consistent, with exceptions in January (Mallard Island and Collinsville) and March (Mallard Island). Results are reported as an upward trend (up), a downward trend (down), or no trend.

Discussion

Although the underlying data presented in this work were available in different documents or electronic sources, the cleaning of the raw data and integration into a single data set of daily average salinity in San Francisco Bay provides a unique perspective on the changes that have occurred over the past nine decades. This period has seen unprecedented anthropogenic change (e.g., land-use, water diversions, and reservoir construction) and significant hydrologic variability, including major floods and multiyear droughts. Additional drivers over the 20th and early 21st centuries include sea level rise and shifts in precipitation, snow accumulation, and runoff patterns. Understanding salinity behavior in this region is of general significance because of the ecological importance of the San Francisco estuary on the Pacific coast and because of the economic significance of the water withdrawals from the Delta that are the single largest source of California’s water supply. These data allow direct examination of the salinity responses to historical events and provide a basis for (1) relating salinity conditions in the current severe California drought to similar conditions that occurred in the past, and (2) refining existing models and exploring future responses in the combined human-hydrologic system, as society adapts to changing natural dynamics and environmental requirements (embodied in the new science of sociohydrology, Sivapalan et al. 2012). Improved understanding of processes affecting the salinity in the western Delta will enhance future management of the upstream reservoirs, withdrawals, and estuarine habitat quality. Key observations from the data evaluation follow.
The construction of upstream water storage and increased in-basin and out-of-basin water use has affected the isohaline positions in different ways, depending on season and water year class. For example, X2 position exhibits less intraannual variability in the postproject period than it did in the preproject period. Postproject X2 position is typically further upstream (i.e., higher) in wet months (February through May) of dry and critically dry years and further downstream (i.e., lower) in the dry months of August and September. This reduction in dry year variability is a straightforward result of reservoirs being operated to store water in wet periods and to release water during dry periods, thus damping the variation in Delta salinity. At the other hydrologic extreme, in wet years, flows are sufficiently high that reservoir operations have less effect on the Delta salinity gradient, resulting in great similarity between preproject and postproject X2 position.
The monthly trend evaluation for the entire period of record shows statistically significant increases in X2 position from November through June and statistically significant decreases in August and September. When the preproject and postproject periods are evaluated separately, important differences emerge. The preproject period is characterized by a statistically significant decreasing trend in X2 position from August through December, reflecting project objectives to maintain freshwater conditions in the Delta during the irrigation season and to evacuate reservoir storage in the fall for winter flood control operations. The postproject period is characterized by a statistically significant increase in X2 position from September through December, reflecting increasing in-basin use and Delta exports. These observations make clear the value of utilizing data from the entire period of record to assess changes in the salinity regime of the estuary. Much of the published literature on X2 and its relationship to various biological indicators is based on data collected over limited periods, typically spanning the mid-1960s to the present. Although it is recognized that such analyses are limited by lack of available biological data prior to the 1960s, conclusions drawn from this partial time interval should be evaluated in light of the more comprehensive description of the estuary’s salinity regime provided in this paper.
Salinity trends, as measured by specific conductance at fixed locations, are broadly consistent with detected trends in X2 position and the conceptual model of increasing salinity with decreasing freshwater flows and with greater proximity to Golden Gate. However, salinity response to flow trends is not uniform along the estuary. Flow trends in high flow months are more likely to translate into detectable salinity trends at downstream (higher salinity) locations, and flow trends in low flow months are more likely to translate into detectable salinity trends at upstream (lower salinity) locations. For example, detection of statistically significant long-term salinity trends was limited to three months at Emmaton [an upstream location (Table 5)], compared with ten months of statistically significant long-term X2 trends (interpolated—Table 3, trends for 1922–2012). Antecedent outflows are often sufficiently high that, at upstream locations such as Emmaton, salinity is not sensitive to modest changes in outflow, i.e., ∂S/∂G is small. The foregoing observation demonstrates the limitations of using a single location for evaluating salinity trends in the estuary and argues for the use of a Lagrangian approach, i.e., evaluating isohaline trends derived from multiple stations.
The X2 time series reported here integrates the effects of multiple drivers, some of which act over decades, and thus affirms the importance of considering longer-term records in defining baselines or targets for defining environmental goals and assessing changes. The periods and statistical analyses presented here are illustrative, and alternative periods or seasons could be considered to examine the response of the system to specific drivers that have the potential to effect isohaline position in the estuary. The data integration presented through this work serves as a foundation for the continuing analysis of salinity behavior in the San Francisco Bay and Delta, anticipating continued interest in the health of the Delta ecosystem in response to anthropogenic and other stressors. The findings presented in this paper are influenced by the data and the cleaning procedure employed, all of which are made available electronically (refer to the “Supplemental Data” section). Future work will consider alternative modeling approaches and statistical analyses to expand on the evaluation of how and why salinity trends in the San Francisco Bay and the Delta have changed over time.

Supplemental Data

Supplemental data files containing Fig. S1, Tables S1–S5, Appendix S1, and electronic data files of salinity and interpolated X2 value are available online in the ASCE Library (www.ascelibrary.org).

Supplemental Materials

File (supplemental_data_wr.1943-5452.0000617_hutton1.pdf)
File (supplemental_data_wr.1943-5452.0000617_hutton2.xlsx)
File (supplemental_data_wr.1943-5452.0000617_hutton3.xlsx)

Acknowledgments

This work was made possible through funding from the San Luis and Delta Mendota Water Authority and the State Water Contractors. We thank the California Department of Water Resources (specifically Tara Smith, Eli Ateljevich, and Joey Zhou) for supporting this study through developing an independent data cleaning procedure for the CDEC salinity data and conducting a Delta hydrodynamic and water quality transport simulation to allow computation of daily ratios between higher high tide salinity and daily average salinity. We also thank three anonymous reviewers for their constructive comments.

References

CDEC (California Data Exchange Center). (2013). 〈http://cdec.water.ca.gov〉 (Jul. 2013).
Google Scholar
CDWR (California Dept. of Water Resources). (2015a). “DSM2 documentation.” 〈http://baydeltaoffice.water.ca.gov/modeling/deltamodeling/models/dsm2v6/documentation.cfm〉 (Mar. 9, 2015).
Google Scholar
CDWR (California Department of Water Resources). (2015b). “DWR Bulletins and Publications.” 〈http://www.water.ca.gov/waterdatalibrary/docs/historic/bulletins.cfm〉 (Oct. 22, 2015).
Google Scholar
CSWRCB (California State Water Resources Control Board). (1999). “Water right decision 1641.” 〈http://www.swrcb.ca.gov/waterrights/board_decisions/adopted_orders/decisions/d1600_d1649/wrd1641_1999dec29.pdf〉 (Jun. 11, 2014).
Google Scholar
Denton, R. A. (1993). “Accounting for antecedent conditions in seawater intrusion modeling—Applications for the San Francisco Bay-Delta.” Proc., 1993 Hydraulic Division National Conf., H. W. Shen, ed., ASCE, San Francisco, 821–826.
Google Scholar
Denton, R. A. (2015). “Delta salinity constituent analysis.” State Water Project Contractors Authority.
Google Scholar
Denton, R. A., and Sullivan, G. D. (1993). “Antecedent flow-salinity relations: Application to Delta planning models.” Contra Costa Water District (CCWD) Technical Memorandum.
Google Scholar
Enright, C., and Culberson, S. D. (2009). “Salinity trends, variability, and control in the northern reach of the San Francisco Estuary.” San Francisco Estuary Watershed Sci., 7(2), 1–28.
Crossref
Google Scholar
Feyrer, F., Newman, K., Nobriga, M., and Sommer, T. (2011). “Modeling the effects of future outflow on the abiotic habitat of an imperiled estuarine fish.” Estuaries Coasts, 34(1), 120–128.
Crossref
Google Scholar
Feyrer, F., Nobriga, M. L., and Sommer, T. R. (2007). “Multidecadal trends for three declining fish species: Habitat patterns and mechanisms in the San Francisco Estuary, California, USA.” Can. J. Fish. Aquat. Sci., 64(4), 723–734.
Crossref
Google Scholar
Fox, J. P., Mongan, T. R., and Miller, W. J. (1990). “Trends in freshwater inflow to San Francisco Bay from the Sacramento-San Joaquin Delta.” J. Am. Water Resour. Assoc., 26(1), 101–116.
Crossref
Google Scholar
Fox, J. P., Mongan, T. R., and Miller, W. J. (1991). “Long-term annual and seasonal trends in surface salinity of San Francisco Bay.” J. Hydrol., 122(1), 93–117.
Crossref
Google Scholar
Gilbert, R. (1987). Statistical methods for environmental pollution monitoring, Van Nostrand Reinhold, New York, 320.
Google Scholar
Gross, E. S., MacWilliams, M. L., and Kimmerer, W. J. (2009). “Three-dimensional modeling of tidal hydrodynamics in the San Francisco Estuary.” San Francisco Estuary Watershed Sci., 7(2), 1–37.
Crossref
Google Scholar
Harder, J. A. (1977). “Predicting estuarine salinity from river inflows.” J. Hydraul. Div., 103(HY8), 877–888.
Google Scholar
Hutton, P. (2014). “Delta salinity gradient (DSG) model, version 1.0 Model Documentation, Metropolitan Water District of Southern California.” Technical Memorandum.
Google Scholar
Jassby, A. D., et al. (1995). “Isohaline position as a habitat indicator for estuarine populations.” Ecol. Appl., 5(1), 272–289.
Crossref
Google Scholar
Kimmerer, W. J., Gross, E. S., and MacWilliams, M. L. (2009). “Is the response of estuarine nekton to freshwater flow in the San Francisco Estuary explained by variation in habitat volume?” Estuaries Coasts, 32(2), 375–389.
Crossref
Google Scholar
Lund, J., Hanak, E., Fleenor, W., Howitt, R., Mount, J., and Moyle, P. (2007). Envisioning futures for the Sacramento-San Joaquin Delta, Public Policy Institute of California, San Francisco.
Google Scholar
Monismith, S. G., Kimmerer, W., Burau, J. R., and Stacey, M. T. (2002). “Structure and flow-induced variability of the subtidal salinity field in northern San Francisco Bay.” J. Phys. Oceanogr., 32(11), 3003–3019.
Crossref
Google Scholar
Moyle, P. B., Lund, J. R., Bennett, W. A., and Fleenor, W. E. (2010). “Habitat variability and complexity in the upper San Francisco Estuary.” San Francisco Estuary and Watershed Science, 8(3), 1–24.
Crossref
Google Scholar
Reed, D., et al. (2014). “Workshop on Delta Outflows and related stressors: Panel summary report, Delta Science Program.” 〈http://deltacouncil.ca.gov/sites/default/files/documents/files/Delta-Outflows-Report-Final-2014-05-05.pdf〉 (Oct. 13, 2015).
Google Scholar
Roy, S. B., Rath, J., Chen, L., Ungs, M. J., and Guerrero, M. (2014). “Salinity trends in Suisun Bay and the Western Delta (October 1921–September 2012).” San Luis and Delta Mendota Water Authority and State Water Contractors.
Google Scholar
Shellenbarger, G. G., and Schoellhamer, D. H. (2011). “Continuous salinity and temperature data from San Francisco Estuary, 1982–2002: Trends and the salinity-freshwater inflow relationship.” J. Coastal Res., 27(6), 1191–1201.
Crossref
Google Scholar
Sivapalan, M., Savenije, H. H., and Blöschl, G. (2012). “Socio-hydrology: A new science of people and water.” Hydrol. Processes, 26(8), 1270–1276.
Crossref
Google Scholar
Stahle, D. W., Therrell, M. D., Cleaveland, M. K., Cayan, D. R., Dettinger, M. D., and Knowles, N. (2001). “Ancient blue oaks reveal human impact on San Francisco Bay salinity.” Eos Trans. Am. Geophys. Union, 82(12), 141–145.
Google Scholar
USFWS (U.S. Fish and Wildlife Service). (2008). “Formal endangered species act consultation on the proposed coordinated operations of the central valley project (CVP) and state water project (SWP).” 〈http://www.fws.gov/sfbaydelta/documents/SWP-CVP_OPs_BO_12-15_final_OCR.pdf〉 (Feb. 5, 2013).
Google Scholar
USGS (U.S. Geological Survey). (2013). 〈http://ca.water.usgs.gov/projects/baydelta/〉 (Jun. 8, 2013).
Google Scholar

Information & Authors

Information

Published In

Go to Journal of Water Resources Planning and Management
Journal of Water Resources Planning and Management
Volume 142 • Issue 3 • March 2016

Copyright

This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/.

History

Received: Oct 1, 2014
Accepted: Sep 3, 2015
Published online: Dec 18, 2015
Published in print: Mar 1, 2016
Discussion open until: May 18, 2016

Authors

Affiliations

Paul H. Hutton
Principal Engineer, Bay Delta Initiatives, Metropolitan Water District of Southern California, 1121 L St., Suite 900, Sacramento, CA 95814-3974.
View all articles by this author
John S. Rath
Environmental Scientist, Tetra Tech Research and Development, 3746 Mount Diablo Blvd., Suite 300,Lafayette, CA.
View all articles by this author
Limin Chen
Principal Engineer, Tetra Tech Research and Development, 3746 Mount Diablo Blvd., Suite 300, Lafayette, CA 94549.
View all articles by this author
Michael J. Ungs
Principal Scientist, Tetra Tech Research and Development, 3746 Mount Diablo Blvd., Suite 300, Lafayette, CA 94549.
View all articles by this author
Sujoy B. Roy, A.M.ASCE [email protected]
Director, Tetra Tech Research and Development, 3746 Mount Diablo Blvd., Suite 300, Lafayette, CA 94549 (corresponding author). E-mail: [email protected]
View all articles by this author

Metrics & Citations

Metrics

Citations

Download citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by

View Options

Media

Figures

Other

Tables

Share

Share

Copy the content Link

Share with email

Email a colleague

Share