Patterns, pace, and processes of water‐quality variability in a long‐studied estuary

Environmental time series have rich information content that is invaluable for measuring and understanding changes over time and guiding policies to manage change. I extracted information from measurements of 10 water‐quality constituents in upper San Francisco Bay from 1975 to 2016, one of the longest observational records in a U.S. estuary. Changes were detected at every time scale captured by monthly sampling. Long‐term trends included increased ammonium (+53%), nitrate + nitrate (+50%), silicate (+14%), Secchi depth (+42%), and decreased chlorophyll a (Chl a) (−74%) and suspended particulate matter (−45%). Changes at the decadal scale included abrupt shifts (Chl a, nitrate + nitrite) and oscillations between shorter trends of increase and decrease (Secchi depth, phosphate). Long‐term trends were not expressed equally across all seasons, and seasonal patterns of change varied across constituents. These examples illustrate key features of environmental variability at the land–sea interface: (1) water‐quality components change continually at time scales from months to decades; (2) patterns of seasonal, multiyear, and multidecadal change are complex and vary across constituents; (3) primary drivers of change are freshwater inflow, the master regulator of estuarine dynamics, and human activities such as river damming, water diversions, wastewater discharge, environmental policies, and species introductions; (4) extracting the full information content of time series requires multiple analyses, each revealing a different layer of insight into how changes develop over time; (5) water‐quality variability is nonstationary, so future changes cannot be forecast reliably; (6) repeated observation is an essential method of Earth system science with applications in the design and performance measures of environmental policies.

Repeated observations are an essential method of Earth system science. They detect changes over time, measure rates of change, and provide clues about their underlying causes. This is the era of the Anthropocene where humans dominate the planet (Vitousek et al. 1997), and the pace of change in some ecosystems is faster than we could have imagined a generation ago (Cloern et al. 2016). We now recognize that planet Earth has limits in the provisions required to sustain humanity, and for some those limits have been exceeded (Röckstrom et al. 2009). A grand challenge of science application is to use discoveries about past and ongoing changes measured in observational time series to understand, anticipate, and guide strategies for mitigation or adaptation to changes that will unfold in the future.
The iconic example of repeated observations as a method of Earth system science is the record of atmospheric CO 2 measurements at Mauna Loa initiated by Charles David Keeling in the late 1950s ( Fig. 1). This time series first documented a seasonal pattern of atmospheric CO 2 tied to the seasonal cycles of photosynthesis and respiration on land. After a decade of measurements, the record revealed a second pattern of increasing atmospheric CO 2 that has accelerated over time and is tied to fossil-fuel burning. This single time series has had monumental significance as evidence that humanity is changing atmospheric chemistry and, as a result, modifying the climate system (Melillo et al. 2014). This discovery has motivated international efforts to assess the pace and consequences of global climate change (IPCC 2014), and multinational agreements to reduce this threat to habitability of planet Earth by humans.
Changes detected in long-term observations have motivated many other kinds of policies and documented some impressive successes, such as reduced atmospheric nitrogen deposition in estuarine watersheds (Eshleman and Sabo 2016); recovery of fish and whale populations decimated by fishing (Richards and Rago 1999;Lyons 2016), brown pelicans extirpated by pesticides (Holm et al. 2003), and forest ecosystems damaged by acid precipitation (Likens et al. 1996); improved air quality in the UK after passage of its 1956 Clean Air Act (Singh et al. 2017); and returns of oxygen and fish communities to rivers after upgrading sewage treatment processes (Andrews and Rickard 1980). Every multidecadal observational record is valuable because of its rich information content. That information is essential for understanding the human dimension of Earth-system change, and it is a critical component of the evidence base required to guide environmental policy (Likens 2010). This article is about change in ecosystems where rivers and oceans meet, and its purpose is twofold. First, I use 42 yr of monthly measurements in upper San Francisco Bay to illustrate the diverse patterns, rates, and processes of variability characteristic of estuaries. Second, I use this example to show how the information content of long time series can be extracted from a suite of analyses, each revealing a different layer of information about how changes develop over time. This case study is intended to show how: (1) the dynamics of estuaries are driven by climate variability and diverse human activities; and (2) trend analysis, alone, can give misleading and incomplete understanding of those dynamics. Extraction of the full information content of long-time series requires use of multiple analysis types.

A multidecadal observational program where rivers and oceans meet
California's Central Valley and slopes of its bounding coastal range and Sierra Nevada Mountains form a 163,000-km 2 watershed. Runoff from that watershed flows into the Sacramento, San Joaquin and smaller rivers that converge to form an inland Delta, once a 14,000-km 2 tidal wetland system that has been transformed into a diked agricultural landscape (Atwater et al. 1979). Inflows to the Delta provide much of California's water supply to urban centers and agricultural production in the Delta and Central Valley. Pumped exports from the Delta to the State Water Project (SWP, see Fig. 2) and federal Central Valley Project (CVP) diverted a mean 18% of Delta Inflow during the period 1975-2016, resulting in reduced Outflow from the Delta to San Francisco Bay (Cloern and Jassby 2012). As a condition for permitting these exports, the State Water Quality Control Board established water-quality objectives and requirements for a comprehensive monitoring program "to provide information for water resource management in compliance with flow-related water-quality standards." Thus was launched the Environmental Monitoring Program (EMP) of the Interagency Ecological Program, one of the longest sustained programs of water-quality monitoring in a U.S. estuary.
The EMP includes boat-based collection of water samples at fixed sampling sites distributed across the Sacramento and San Joaquin Rivers, the interior Delta, and upper San Francisco Bay. I analyzed a 42-yr record of measurements from Sta. D8, situated in the channel of Suisun Bay downstream of the Sacramento-San Joaquin River confluence (Fig. 2). This site has a complete record of water-quality measurements made once or twice monthly since 1975, and includes constituents not sampled at all other sites, such as chlorophyll a (Chl a) and nutrient concentrations. This site is characteristic of transitional zones between land and sea where variability of water quality and biological communities is tied to seasonal and annual fluctuations of river inflow (Cloern et al. 2018). The period of record includes years of extreme (1977) and prolonged (2012-2015) drought, record-high annual inflow (1983), and extreme floods (e.g., March 1983, January 1997, February 1998. It therefore provides an opportunity to measure estuarine responses to a wide range of hydroclimatic forcings. It also was a period of population growth and deepening of the human footprint through activities such as wastewater disposal, river damming, water management, nutrient enrichment, and introductions of non-native species (Cloern and Jassby 2012). Measurements repeated over decades provide a record for identifying patterns of water-quality variability over time, and processes of that variability associated with the climate system and human activities.

Data sources and preprocessing
I accessed water-quality data for the period January 1975 through June 2015 at the EMP webpage http://www.water.ca. gov/iep/products/data.cfm (accessed 19 July 2018). Data for the period July 2015 through 2016 were provided by Jenna Rinde of the California Department of Water Resources (DWR, personal communication 18 October 2017). The record contains surface (0.9 m) measurements made at Sta. D8 (Fig. 2) on 621 dates. The number of missing entries for most constituents was < 15. Exceptions were Chl a (32 missing values), ammonium (62 missing values), and phosphate (101 missing values). I screened the data by examining values that could be classified as far outliers, i.e., more than three times the interquartile range above the 75 th percentile (Tukey 1977). If extreme values were consistent with other measurements (e.g., coincident far outliers of both dissolved oxygen

Cloern
Water-quality changes in a long-studied estuary S193 [DO] and Chl a), I retained them. Otherwise, I removed one value of Secchi depth > 120 cm, one value of nitrate + nitrite > 71 μM, two values of ammonium > 20 μM, and three values of suspended particulate matter (SPM) > 116 mg L −1 . From the screened data, I computed monthly mean measurements of Secchi depth, SPM, Chl a, DO, electrical conductivity, temperature, and concentrations of ammonium, nitrate + nitrite, phosphate, and silicate. Salinity was calculated from conductivity and temperature (Fofonoff and Millard 1983). I compiled monthly mean values of these same constituents from measurements made at the Sacramento-San Joaquin River confluence (Sta. D4, Fig. 2) to explore synchrony of water-quality variability in the rivers and downstream estuary.
Biological monitoring by DWR includes sampling and analysis of benthic invertebrate populations at Sta. D7 in the shallow subtidal region of Suisun Bay adjacent to channel Sta. D8 (Fig. 2). I accessed data online (http://www.water.ca.gov/bdma/meta/ benthic/data.cfm, accessed 02 January 2018) and computed monthly mean abundance of clams Mya arenaria and Potamocorbula amurensis. Their summed abundance was used as an index of potential phytoplankton grazing rate by benthic filter feeders. P. amurensis was introduced in late 1986 and has been abundant since; the marine clam M. arenaria colonized the upper estuary during the extreme drought of 1977 (Nichols 1985).
Daily flow data for the period 01 January 1975 through 30 September 2016 were accessed at the DWR website http://www.water.ca.gov/dayflow/output/Output.cfm on 21 September 2017. From this record, I computed mean monthly Inflow to the Delta, Export from the Delta to the SWP and CVP (see Fig. 2), and Outflow from the Delta to San Francisco Bay. Descriptions of sampling, analytical, and computational methods for these data are available at the DWR webpages above.
An important point source of nitrogen and phosphorus to the estuary is effluent from the Sacramento Regional Wastewater Treatment Facility discharged into the Sacramento River (Jassby 2008). Measurements of effluent flow and concentrations of total phosphorus (TP), Total Kehldal Nitrogen (TKN), and ammonium were made at frequencies that varied from weekly to monthly. These data were provided by Dr. Timothy Mussen of the Sacramento Regional County Sanitation District (personal communication, 04 October 2017). From these, I computed mean monthly loadings of TKN and TP (tons d −1 ) as indices of wastewater N and P sources. This treatment plant does not include a nitrification process, so most of the TKN loading was as ammonium nitrogen. These loadings are available beginning October 1986 for TP and January 1987 for TKN.

Methods and steps of data analysis Patterns and rates of change
The data records analyzed here include monthly time series of 10 water-quality constituents and six measures of hydroclimatic and human drivers of change: water Inflow to and Export from the Delta, Outflow from the Delta to San Francisco Bay, wastewater loadings of N and P, and abundance of clams as top-down regulators of phytoplankton biomass accumulation. All analyses were done in the R environment (R Development Core Team 2012) and most with package wql, a maintained version of now-archived package wq developed to explore multivariate environmental time series (Jassby and Cloem 2016). A common set of steps was used to identify key patterns of variability in each data record: 1. Plot the time series with a function to smooth monthly scale variability and aid visualization of multiyear patterns including trends. I used a locally weighted smoothing function (R function loess) with span 0.2. The span parameter prescribes the proportion of data used to compute means around each point in a series. As an exploratory step, experimentation with different spans can be insightful. I used span 0.2, producing smoothed values from measurements centered around 8.4-yr segments (0.2 span × 42 yr) of the full series. 2. Decompose the time series into annual and seasonal components (wql function decompTS) to (1) isolate and more clearly visualize the annual patterns, (2) identify a mean seasonal (monthly) pattern, and (3) compare the magnitudes of annual and seasonal variability. 3. Measure trends of change over time, using the Mann Kendall test (function mannKen) on annualized series, and the Seasonal Kendall test (function seaKen) on monthly series to remove any effects of seasonality. Nonparametric tests were used because they detect monotonic (not necessarily linear) trends, are resistant to the effects of extreme values, and are robust against departures of data from normality (Hirsch et al. 1982(Hirsch et al. , 1991. 4. Apply the Seasonal Kendall test with a 14-yr running window (function seaRoll) to detect oscillations-i.e., runs within the series having significant trends but of opposite signs. I chose a 14-yr window somewhat arbitrarily, but this does represent a third of the record length and allows series of sufficient duration to look for significant trends over periods shorter than the observational record. 5. Use the nonparametric Pettitt test (Pettitt 1979) (function pett) to determine if trends were the result of an abrupt shift rather than a steady period of increase or decrease. This analysis can be informative because the pattern of long-term change can give clues about its underlying cause. The Pettitt test finds, if present, the location of a single change point in a time series using the Mann-Whitney statistic to determine if samples before and after the change point come from the same distribution (Jassby and Cloern 2016). Most techniques of change-point detection were developed for independent observations. Their power of detection weakens and the likelihood of false-positive detections increases when the series have autocorrelation (Beaulieu et al. 2012). The monthly time series have high autocorrelation, but the series of annual mean values have no or weak autocorrelation so I applied the Pettitt test on annualized data. 6. Estimate the trends and their significance levels for each month, using function seasonTrend. This analysis can be insightful because some series with no overall trend do have significant trends of change for individual months, and because trend signs (direction) can vary across months.
Processes of change I used regression models to identify factors most strongly associated with annual variability of each water-quality constituent. The objective was not the best-fitting regression model, but rather a simple model that identifies a few environmental factors that explain a large fraction of that variability. I used ordinary linear regression implemented as function ols in the rms package (Harrell 2017). Nonlinear relationships of annualized data were explored using restrictive cubic splines. The model-building strategy was to first consider the full suite of factors that might drive change in a particular constituent. For example, long-term changes in Chl a might be driven by changes in N, P, or Si concentrations or inputs, temperature, salinity, outflow or the seasonal timing of outflow, turbidity (Secchi depth), or top-down control by clams. For each constituent, a full model was built allowing for nonlinear relationships with all factors considered. Factors that had no significant effect were removed, and factors with no significant nonlinear terms were recast as linear. The process of model simplification was iterative and continued to find the one or two factors that explained more than 50% of the annual deviance. Clam abundance was treated as a categorical variable and assigned a value of either 1 (abundance above the 20 th percentile = 270 m −2 ), or 0 (< 270 m −2 ). Model validation included steps to determine if residuals were normally distributed (using the Shapiro-Wilk test) and homogeneous (by testing for a linear relationship between model residuals and fitted values). When residuals departed from these expectations, I explored model improvements by log transformation of the predicted and/or predictor variables. When additive models produced little explanatory power, I explored models with interactions between predictor variables.
The regression models were used to measure and compare effects of individual predictor variables on annual variability of each water-quality constituent. Effects were recovered from the summary function of ols models produced in the rms package. The method measures effect size as the range of modelpredicted values along the interquartile range of the target predictor variable while holding other predictor variables at their median.

Patterns of variability, rates of change
Step 1: Visual inspection Measurements of atmospheric CO 2 at Mauna Loa ( Fig. 1) show a monotonic increase from January1958 through October 2017, with small-amplitude seasonal oscillations. The pattern is clear because it reflects only two processes: steadily increasing CO 2 emissions from fossil-fuel combustion (a large signal) and the annual cycle of net photosynthesis on land (a small signal).
The key information content of this time series can be extracted simply by visual inspection. However, the information content of water-quality time series is usually not so easily extracted because the patterns are more complex and varied. Ten examples from the San Francisco Bay data set are shown in Fig. 3, illustrating a diversity of pattern types. Some series have oscillatory patterns. Monthly salinity ranged from 0 to 21, reflecting the dynamic nature of mixing between freshwater and seawater in estuaries. The series has oscillations as periods of high salinity centered around 1977,1990, and the most recent years, and periods of low salinity in the early 1980s and early 1990s. Phosphate also has an oscillatory pattern, with an increase from the early 1980s to early 1990s followed by a decade of decrease. The temperature and DO series have large monthly variability around small or no long-term trends, illustrating that seasonal variability is more important than decadal variability for some constituents. The Chl a pattern illustrates a shift to low concentrations and damped seasonal variability in the late 1980s. Series of other constituents suggest trends of increase (Secchi depth, nitrate + nitrite, ammonium, silicate) or decrease (SPM). None of these trends are linear, or even monotonic, and some are difficult to visualize without the loess fits because of high-frequency monthly variability. The diversity and complexity of these patterns suggest that annual variability of water quality, unlike variability in atmospheric CO 2 , is caused by many processes and that for some constituents it is small relative to seasonal variability.
Step 2: Decomposition into annual and seasonal components Monitoring programs often sample at monthly or quarterly frequency, so records sustained over decades contain information about the relative magnitudes of annual vs. seasonal variability. Trends in annual variability are signals of long-term change, while seasonal variability can mask those signals and make them difficult to visualize or detect. The first few years of the Keeling CO 2 curve revealed what appeared, at the time, to be high-amplitude seasonal variability. As measurements continued for a decade, it became apparent that the seasonal component was small relative to the trend of increasing atmospheric CO 2 . Comparison of variability at the two scales provides a context for understanding long-term trends and their significance. I used a simple additive model to decompose monthly time series into four components using wql function decompTS: where c ij is the measured value for year i and month j; C is the overall mean of the series; y i is the annual effect, computed as the annual mean Y i divided by the overall mean C; m j is the mean monthly effect, adjusted for year-to-year variability; and ε ij is the residual (Cloern and Jassby 2010).
Results of this analysis revealed wide differences across constituents in the relative magnitudes of their annual and monthly variability. The annual component of temperature variability was small relative to the monthly component (Fig. 4A). The annual and monthly components of ammonium concentration were of comparable magnitude (Fig. 4B), and the annual component of SPM variability was larger than the seasonal component (Fig. 4C). Plots of the isolated annual component can be informative because they reveal details of trends not easily seen in time series plots, such as the anomalous high ammonium concentrations in 1977 (year of extreme drought) and 1994, and the nearly continuous run of positive deviations from the long-term mean beginning in 2000 (Fig. 4B). Plots of the seasonal component show the mean monthly pattern and identify constituents with either weak (SPM) or strong (ammonium) seasonality. I compared annual and seasonal components and their ratios for the 10 constituents measured in upper San Francisco Bay and the Keeling CO 2 series to illustrate the large seasonal variability of estuarine water quality relative to atmospheric CO 2 (Table 1).
Step 3: Overall trends I used two tests for monotonic trends, one on the series of annual mean values (Mann Kendall), and one on the monthly time series (Seasonal Kendall). Conclusions from both methods were identical: there were no trends of change in salinity, temperature, or phosphate over the period 1975-2016, but there were significant trends of decreasing SPM, Chl a, and DO, and increasing Secchi depth, nitrate + nitrite, ammonium, and silicate (Table 2). For these seven constituents with trends, I computed the percent change from the first to last decade of the observational record. Although the DO trend was significant, it was small (2% decrease). However, the > 40% decrease of SPM and increase of Secchi depth, 74% decrease of Chl a, and ≥ 50% increase of nitrate + nitrite and ammonium represent substantial changes in suspended sediments, light penetration, phytoplankton biomass, and inorganic nitrogen concentrations over the past four decades. These changes illustrate the varied and fast-paced changes often measured in estuarine-coastal ecosystems (Cloern et al. 2016). Hirsch et al. (1982) developed the Seasonal Kendall test, used above, as an exploratory step to detect and measure trends in water-quality time series. Environmental time series sustained over decades contain substantially more information than just trends, and further analyses can reveal insightful details of variability not captured by overall trend tests. In particular, trend tests provide no information about the course of change over time. Was it monotonic, abrupt, oscillatory, and was it expressed equally across all seasons? I used three additional analyses to extract further information from these time series.
Step 4: Oscillations The trends reported in Table 2 are from a test of monotonic trend over the full period of observation. However, inspection of those records (Fig. 3) shows that many patterns of variability are not monotonic. The wql function seaRoll detects oscillatory variability by applying the Seasonal Kendall test on user-specified windows as they move along the time series, 1 yr at a time. For example, with a 14-yr window, the positive Sen slope of Secchi depth at year 1977 (Fig. 5A) represents a significant trend over the period [1977][1978][1979][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990]. Four examples of  (Table 2) gives an incomplete description of how Secchi depth changed between 1975 and 2016. The time series of phosphate concentration (Fig. 5B) shows an early period of significant increase, a later period of significant decrease, and a final period of no trend. Thus, the rolling Seasonal Kendall test gives a very different depiction of how phosphate varied than the overall Seasonal Kendall test that indicated no significant change ( Table 2). As we will see below, trends over subsets of time series can give clues about processes of water-quality variability. The Chl a series has an overall negative trend of −0.04 mg m −3 yr −1 (Table 2). However, the rolling Seasonal Kendall test with 14-yr window shows that Chl a decreased significantly over the period 1975-1993, and then increased slowly after 1988 (Fig. 5C). Again, this analysis gives a more complete and insightful depiction of change than the overall trend test. This analysis reveals how the overall Chl a trend developed as a period of fast decrease (−0.2 mg m −3 yr −1 to −0.4 mg m −3 yr −1 ) followed by a period of slow increasedistinctly different from the assumed monotonic pattern of the trend tests used above. As a contrast to the oscillations of water-quality time series, rolling trends of the Mauna Loa CO 2 series show positive trends across all 14-yr segments of the full record (Fig. 5A), diagnostic of a monotonic trend.
Step 5: Abrupt changes One of the important ecological discoveries from multidecadal observations is that changes can occur abruptly as discrete shifts from one state to another (e.g., Anderson and Piatt Cloern Water-quality changes in a long-studied estuary S198 1999). The Pettitt test detected change points in every waterquality time series that had a significant overall trend (Table 3). For example, it detected a 1987 change point in Chl a when mean concentration decreased 5.9 mg m −3 (Table 3). This abrupt shift is well understood. It coincided with the explosive population growth of the non-native clam P. amurensis after it was introduced in late 1986, and has been attributed to grazing control by this filter-feeder (Alpine and Cloern 1992). Visual inspection of the Chl a time series (Fig. 3) suggested the possibility of an abrupt change, and the Pettitt test verified a shift and identified the year of change.
Further corroboration of a state change included the synchronous 1987 shift to higher nitrate + nitrite concentration (Table 3). Phytoplankton uptake is an important sink for inorganic nitrogen in this estuary (Peterson et al. 1985), and the abrupt nitrate + nitrite increase is explained by the synchronous decrease of phytoplankton biomass and N uptake.
Change points detected in series of other constituents, such as DO (Table 3), are not supported with corroborating evidence. Change-detection methods identify points along a time series where statistical properties of the measurement change significantly, but these points do not necessarily reflect abrupt shifts. For example, some tests identify a change point in the middle of series having a linear trend (Reeves et al. 2007). The midpoint of time series analyzed here is 1996, so change points detected around that year (e.g., DO, ammonium) should be interpreted with caution. The Pettitt test also detected a shift in the Keeling CO 2 series at the midpoint of that record (1988), but a discrete shift is not evident in the series (Fig. 1). Visual inspections of SPM, Secchi depth, and silicate series (Fig. 3) do not show the abrupt changes detected by the Pettitt test. Therefore, change-point analysis can reveal important information about trends that is not revealed by simple trend analysis. However, results of those tests should be interpreted with skepticism if the change point is near the series center, if it is not corroborated with other data, is not evident in visual inspection, and if the change has small magnitude.
Step 6: Trends by month Trend tests do not give information about how or if changes vary across seasons. The wql function seasonTrend extracts another level of information from time series by computing signs and magnitudes of trends for individual seasons or months. Six different patterns of seasonal change are illustrated in Fig. 6. The overall trend of decreasing SPM was a result of decreases every month except December and January (Fig. 6A), so there was no pronounced seasonal pattern in the SPM trend. Ammonium increased every month except July and August and, unlike SPM, it did have a pronounced seasonal pattern with largest increases during winter and smallest in summer (Fig. 6B). The overall trend of Chl a decrease was primarily a result of large decreases in summer, June- Cloern Water-quality changes in a long-studied estuary September (Fig. 6C). The trends of silicate were positive each month but they were largest in summer, particularly July and August (Fig. 6D). Salinity had no overall trend, but it increased in autumn, September-December (Fig. 6E). Water export had no overall trend, but there were significant trends of decrease during winter-spring and increase during summer-autumn (Fig. 6F). These examples illustrate that trends of change are not expressed equally across seasons, and that seasonal patterns of change vary among water-quality constituents. Those seasonal patterns can give clues about the processes of change. For example, the absence of seasonal variability in the trend of SPM decrease suggests a process of reduced sediment supply that operates year-round. The ammonium pattern suggests a process of increasing supply but rapid assimilation during summer, perhaps through fast nitrification. The Chl a pattern suggests a process of either reduced production or increased consumption during summer. The patterns of silicate, nitrate + nitrite, and ammonium increase could reflect reduced uptake after phytoplankton biomass decreased in 1987. The salinity increase during autumn could be related to the trend of increased summer-autumn Export. Each of these hypotheses that emerged from seasonal trend analyses can be explored with regression models. Development of those models can be guided by seasonal trend analyses that identify months when changes were large and dominated the long-term trends.

Drivers of change
Estuarine dynamics are strongly tied to climatic variability through seasonal and annual variability of precipitation, runoff, and river inflow. I used daily measurements of Inflow to the Delta to represent this hydroclimatic variability (Fig. 7A). Water management has changed the quantity and timing of river discharge to the estuary, and I used Export to the state and federal water projects as an index of human modification of hydrology (Fig. 7B). Summed abundances of M. arenaria and P. amurensis were used as an index of potential grazing rate by bivalve filter feeders (Fig. 7C). Seasonal, annual, and decadal variability of suspended sediment concentrations in estuaries are controlled by riverine inputs. I used SPM concentrations measured at the confluence of the Sacramento and San Joaquin Rivers (Sta. D4, Fig. 2) as an index of riverine sediment supply (Fig. 7D). The Sacramento Regional Wastewater Treatment Facility is a large point source of nitrogen and phosphorus to the Sacramento River. Monthly measurements of TKN and TP loading are available since 1987 (Fig. 7E,F). I used regression models to measure associations between these climatic and human processes and the water-quality changes described above. The strategy was to focus on changes during those seasons contributing most to long-term trends, guided by the monthly trend tests (Fig. 6).

Salinity
The salinity gradient is the defining feature of estuaries that sets them apart from the rivers and oceans they span (Cloern et al. 2018). It is structured by river input of freshwater that, in this estuary, is determined primarily by Inflow to the Delta (a climate-related process) and water management including upstream consumption and Export from the Delta (Cloern and Jassby 2012). Analyses above showed that salinity has increased during autumn (Fig. 6E), so I constructed annual series of mean September-December measurements to explore relationships between salinity, Inflow, and Export during that season. The best-fitting simple regression model described salinity as a nonlinear function of Inflow and a linear function of the ratio Export : Inflow, i.e., the fraction of Inflow diverted from the estuary (Fig. 8A). This regression model explained  Cloern Water-quality changes in a long-studied estuary 84% of the autumn salinity variability, with the largest effect (−5.2) from Inflow, and second largest effect (+2.0) from the ratio Export : Inflow (Table 4). Water Export increased during autumn months (Fig. 6E), so the salinity increase appears to be a response to changing water-management practices and, in particular, increases in the fraction of Inflow exported during autumn (Fig. 8A, bottom panel). The salinity example, and those that follow, illustrate how change in estuaries is driven by processes tied to the climate system and to human activities, and how their effect sizes can be compared.
Chl a Phytoplankton biomass in estuaries is regulated by riverdriven transport processes and the balance between growth and grazing rates. Chl a at Sta. D8 decreased abruptly in 1987 (Table 3), and the largest declines occurred in summer-May through September (Fig. 6C). The best-fitting model of summer Chl a had an interactive effect between Outflow and clam abundance. This interaction was expected because the relationship between Chl a and Outflow changed after the Potamocorbula introduction (Jassby 2008). Most (80%) of that variability was associated with Outflow and clam abundance treated as a categorical variable (Fig. 8B). No other variables, such as temperature, turbidity (Secchi depth), nutrient concentrations, or loadings had comparable effects on summer Chl a.
Summer phytoplankton blooms occurred regularly in this region of the estuary during the 1970s and 1980s. These were attributed to the accumulation of phytoplankton cells within Blue lines are partial-residual plots with 95% confidence intervals (gray band). Table 4. Summary of regression models used to identify processes of variability in the seven water-quality constituents having significant trends of change over the period 1975-2016. Models were built to explore annual variability for the season (months) when trends were largest (see Fig. 6). Effect sizes are differences in model predictions between the 75 th and 25 th percentile of the target predictor variable as the other is fixed at its median. 200-300 m 3 s −1 (Fig. 8B), primarily through its effect on residence time (Jassby 2008). A summer bloom did not develop in 1977, a year of extreme drought when salinity increased high enough that the upper estuary was colonized by marine benthic filter feeders including M. arenaria (Nichols 1985). The 1977 anomaly was the first evidence that bivalve grazing can suppress blooms in this estuary. Summer blooms disappeared completely after P. amurensis became established as a permanent resident in 1987. The regression model (Table 4) shows that the effect of clam presence was large (−13.1) compared to the outflow effect (−2.5). The nearly fourfold loss of phytoplankton biomass (Table 2) and production restructured biological communities and altered nutrient dynamics, establishing what has become an iconic example of abrupt disturbance by the introduction of a non-native species (Cloern and Jassby 2012). Regression models allow us to compare the effects of that human disturbance with the effect of varying freshwater inflow.

Silicate
The rivers flowing into San Francisco Bay have high (3 00 μM) silicate concentrations derived from dissolution of silicate minerals in the watershed. Spatial distributions in the estuary are determined by mixing of high-silicate river water with lower-silicate ocean water; concentrations decrease along the river-ocean salinity gradient (Cloern et al. 2018). Summer phytoplankton blooms that developed prior to the Potamocorbula introduction were dominated by diatoms, and silicate concentrations decreased as these blooms developed during summer (Peterson et al. 1985). After the summer blooms disappeared, diatom uptake decreased leading to a long-term trend of increasing silicate concentrations. The rate of silicate increase was largest in summer (June-October, Fig. 6D), the season of largest Chl a decrease.
A regression model verified the effects of Outflow and Chl a that, together, explained 65% of the mean June-October silicate variability (Fig. 8C). The Chl a effect was nonlinear and negative, with highest silicate concentrations when summer Chl a concentration was less than about 3 mg m −3 -i.e., during the post-Potamocorbula era. The Outflow effect was hyperbolic, with a steep slope of increasing silicate concentrations up to summer Outflow ≥ 300 m 3 s −1 (Fig. 8C) when salinity is low and the fraction of high-silicate river water is high. The patterns in Fig. 8C reflect control of summer silicate by two key processes-the riverine supply, and phytoplankton sink. Analyses of model results showed that the outflow effect of +37.9 was nearly double the phytoplankton (Chl a) effect of −21.3 (Table 4).

Ammonium
While the silicate source to San Francisco Bay is river inflow, the ammonium source is primarily municipal wastewater (Jassby 2008). Comparison of silicate and ammonium variability provides an opportunity to explore the differing patterns of nutrients delivered by runoff and point sources. The trends of increasing ammonium were significant for most months (Fig. 6B), so I explored models to explain variability of annual mean ammonium concentration. The best-fitting simple model explained 59% of ammonium variability with two processes-Outflow and wastewater inputs of TKN (Fig. 8D). The size effects of Outflow (−1.0) and TKN loading (+0.8) were comparable ( Table 4). The relationship with Outflow is nonlinear and negative, with highest ammonium concentrations during dry years and lowest during wet years (Fig. 8D). This is opposite the Outflow-silicate relationship because river inflow is a source of silicate but it dilutes wastewater-derived ammonium. The shape of nutrient-Outflow relationships can therefore give information about the relative importance of riverine and point-source inputs. Ammonium concentration in the estuary was significantly and linearly related to TKN loading (Fig. 8D), so the long-term trend of increasing ammonium concentration ( Fig. 4; Table 2) can be attributed to increasing wastewater loading over time (Fig. 7E), but not to changes in Outflow and dilution.

Ammonium at the monthly scale
Analyses of water-quality time series data can yield different insights into processes depending on the time scale considered. As an example, I explored regression models to identify the key drivers of ammonium variability at the monthly scale to compare against analyses at the annual scale. This was motivated by the strong monthly pattern of ammonium trends-largest increases in winter and smallest in summer (Fig. 6B), suggesting a possible influence of seasonal temperature variability. The data series included 346 monthly measurements beginning January 1987. The best fitting model described monthly ammonium variability as a function of three variables-Outflow and TKN loading (as above), and temperature. Monthly ammonium was a nonlinear and negative function of Outflow and a linear positive function of TKN loading (Fig. 9)-the same functional responses observed at the annual scale (Fig. 8D). The Outflow effect was −2.1 and the TKN-loading effect was +1.0. However, at the monthly time scale, temperature had the largest effect (−5.8) through a negative linear relationship. Thus, both annual and monthly variability of ammonium, including trends of increase over time, are associated with variability of TKN loading as a source and Outflow through its dilution effect. However, a stronger temperature effect was revealed by analyzing monthly data. This effect is not evident in analyses of the annual data because annual temperature fluctuations are small relative to the high-amplitude monthly variability (Fig. 4A).
The seasonal temperature effect on ammonium concentration in the estuary is presumably a response to the near-linear relationship between nitrification rate and temperature (e.g., Sudarno et al. 2011). Highest nitrification rates during the warm season would explain the negative relationship between ammonium and temperature (Fig. 9A), the seasonal pattern of lowest ammonium concentrations during summer (Fig. 4B), and the seasonal pattern of ammonium increase over time-largest in winter and near-zero in summer (Fig. 6B). Thus, the long-term trends of ammonium increase in the estuary can be attributed to increased wastewater inputs. But the expression of increased loading is not evident in summer (Fig. 6B), implying that wastewater ammonium is nitrified in the time of transit from the wastewater source to the sampling site in the estuary during the warm, low-flow season. Direct measurements of nitrate production and travel time in the lower Sacramento River support this conclusion (Kraus et al. 2017).

SPM and Secchi depth
The long-term decline of SPM in San Francisco Bay is well documented and understood as a response to a sequence of human activities that began in the 19 th century when hydraulic gold mining in the Sierra Nevada mountains and other watershed disturbances mobilized massive quantities of sediments and increased their supply to San Francisco Bay. The era of increased sediment supply peaked in the early 20 th century, and it was followed by an era of reduced sediment supply, resulting in part from dam construction on all major tributaries (Schoellhamer et al. 2013). The primary source of sediment to the estuary is from the Sacramento River that contributes most of Delta inflow. The trends of SPM increase from 1975 to 2016 were positive each month from February through November (Fig. 6A), so I analyzed SPM concentrations averaged over these months. The best-fitting regression model (adjusted R 2 = 0.67) described annual SPM as a linear function of SPM at the upstream Sta. D4 (an index of the river supply), and a nonlinear function of Outflow. The effects of Outflow (+11.8) and SPM at the river confluence (+9.6) were both positive and of similar magnitude. Therefore, Outflow plays an important role in the annual variability of SPM in the estuary, and the long-term trend of decrease results from the 50% reduction of sediment supply from the Sacramento River since the 1950s. That loss was a result of sediment retention behind dams and gradual removal by storms of the sediments deposited in the river system during the 19 th century (Schoellhamer et al. 2016).
San Francisco Bay is a turbid estuary, and that turbidity derives primarily from light scattering and absorption by the suspended fine-grain minerals delivered by rivers. Most (81%) of the variability of annual-mean Secchi depth is associated with annual-mean SPM concentration, so the trend of increasing Secchi depth is a consequence primarily of decreased SPM. Support for this conclusion is the nearly identical 42% increase of Secchi depth and 45% decrease of SPM over the study period (Table 2). A contributing, but smaller, effect was associated with the trend of decreasing Chl a.

Phosphate
Although phosphate concentration did not change over the long term (Table 2), there was substantial variability of phosphate concentration as an early period of increase during the 1980s, peak concentrations in the early 1990s (Fig. 3) followed by a period of significant decrease (Fig. 5B). Records of phosphorus loading from the Sacramento Regional Wastewater Treatment Facility are not available during the early era of phosphate increase. However, measurements since 1987 showed a steep decrease in wastewater TP loading in the early 1990s that paralleled the downstream phosphate decrease at Sta. D8. This suggests that wastewater effluent is a significant source of phosphorus, as well as ammonium nitrogen, to the upper estuary. The best-fitting regression model (Table 4) explained 73% of annual-mean phosphate concentration as a negative function of Outflow (effect size = −0.74) and a positive function of wastewater TP loading (effect size = +0.25), similar to model results for ammonium. Relative to the Outflow effect, the effect of wastewater TP loading was smaller than the effect of ammonium loading (Table 4). This suggests that the wastewater enrichment of riverine phosphate is smaller than its enrichment of riverine ammonium.
The parallel steep declines of phosphorus loading and concentrations in the estuary during the 1990s (Fig. 3) occurred after it became evident that phosphorus enrichment was degrading water quality of U.S. lakes and rivers. State bans on phosphorus-containing detergents grew in the 1970s and 1980s, and by 1994 the manufacture of phosphorus-based detergents ended (Litke 1999). Thus, the oscillating patterns of phosphate concentrations in San Francisco Bay (and other estuaries and lakes) include eras of increase associated with population growth (Jassby 2008) and decrease after policies were implemented to remove phosphorus from household detergents.

Conclusions
When the IEP Environmental Monitoring Program began in 1975, we could not imagine the changes it would capture over the succeeding four decades. Changes occurred at every time scale captured by monthly sampling. Long-term trends included increased ammonium (+53%), nitrate + nitrate (+50%), silicate (+14%), Secchi depth (+42%), and decreased Chl a (−74%) and SPM (−45%). Changes at the decadal scale included abrupt shifts (Chl a, nitrate + nitrite) and oscillations between shorter trends of increase and decrease (Secchi depth, phosphate). Every constituent varied from year to year, but the magnitude of annual variability ranged from small (temperature, DO) to large (SPM, Secchi depth, Chl a). Long-term trends were not expressed equally across all seasons, and patterns of change by month varied across constituents. For example, salinity increased in autumn, Chl a decreases were largest in summer, while ammonium increases were largest in winter. Thus, changes in estuarine water quality have a diversity of patterns, occur across a spectrum of time scales, and are more complex than the monotonic pattern of increasing atmospheric CO 2 . This complexity arises from the unique array of processes that drive variability in estuaries where forcings from rivers, oceans, the overlying atmosphere, and a deep human footprint converge (Cloern and Jassby 2012).
Changes measured in upper San Francisco Bay were strongly tied to seasonal and annual variability of flow because river discharge is the master regulator of estuarine dynamics through its influence on transport processes, stratification and mixing, salinity, inputs of sediments and nutrients, and dilution of point-source pollutants (Cloern et al. 2018). They were also strongly tied to human activities, including those in the upper watershed (river damming), in the tributary rivers (wastewater discharge, water diversions), within the estuary (species introductions through ship ballast discharge), and after implementation of policies (phosphorus removal from detergents). Regression models can be used to measure the relative importance of flow and human activities as drivers of water-quality change. They show that the human factor is large enough to be detected against a background of high climatic (flow) variability.
San Francisco Bay is a very different ecosystem today from the one we began to study in the 1970s. Changes described here are characteristic of the fast-paced changes occurring in estuaries globally (Cloern et al. 2016). We know with certainty that changes will continue into the future, but there is large uncertainty about how. Results presented here show that water-quality constituents are not stationary-their statistical properties (mean, variance, range) are changing. The challenge of forecasting and planning for future changes is daunting because statistical models have not been developed to forecast future states in a nonstationary world (Milly et al. 2008). Continuity of repeated observations, sustained indefinitely into the future, is therefore critical.