Ecosystem variability along the estuarine salinity gradient: Examples from long‐term study of San Francisco Bay

The salinity gradient of estuaries plays a unique and fundamental role in structuring spatial patterns of physical properties, biota, and biogeochemical processes. We use variability along the salinity gradient of San Francisco Bay to illustrate some lessons about the diversity of spatial structures in estuaries and their variability over time. Spatial patterns of dissolved constituents (e.g., silicate) can be linear or nonlinear, depending on the relative importance of river‐ocean mixing and internal sinks (diatom uptake). Particles have different spatial patterns because they accumulate in estuarine turbidity maxima formed by the combination of sinking and estuarine circulation. Some constituents have weak or no mean spatial structure along the salinity gradient, reflecting spatially distributed sources along the estuary (nitrate) or atmospheric exchanges that buffer spatial variability of ecosystem metabolism (dissolved oxygen). The density difference between freshwater and seawater establishes stratification in estuaries stronger than the thermal stratification of lakes and oceans. Stratification is strongest around the center of the salinity gradient and when river discharge is high. Spatial distributions of motile organisms are shaped by species‐specific adaptations to different salinity ranges (shrimp) and by behavioral responses to environmental variability (northern anchovy). Estuarine spatial patterns change over time scales of events (intrusions of upwelled ocean water), seasons (river inflow), years (annual weather anomalies), and between eras separated by ecosystem disturbances (a species introduction). Each of these lessons is a piece in the puzzle of how estuarine ecosystems are structured and how they differ from the river and ocean ecosystems they bridge.

Estuaries are transitional ecosystems between land and ocean. Their defining feature is the salinity gradient, a spatial pattern in the mixture of water from two distinctly different sources: seawater, and freshwater delivered as land runoff. The salinity gradient is a fundamental reason why the estuary is not just an ecotone-a transitional zone between two separate biomes-but rather an ecosystem class of its own, with unique species assemblages and transport processes. As a result, patterns of water properties and biota along the estuarine axis from freshwater to sea seldom represent the simple passive mixing of two water sources. Here, we examine the nature of and reasons for patterns along the salinity gradient of San Francisco Bay, ". . .a classic example of a coastal plain estuary in which terrestrial freshwater mixes with salt water entering the estuary from the ocean" (Schoellhamer et al. 2016).
Since the beginning of estuarine science, we have recognized the importance of the salinity gradient, using it as a frame of reference for understanding spatial structure and its variability in the form of property vs. salinity plots, or "mixing diagrams" (Stef ansson and Richards 1963;Ketchum 1967;Wollast and De Broeu 1971;Liss and Pointon 1973). The approach is informative and appealing because of its simplicity, exploiting differences in the chemical and physical properties of seawater and freshwater. When endmember concentrations vary more slowly than estuary flushing time and there are no intermediate transformations or sources, then spatial structure is determined primarily by mixing of the two end-members, and properties should vary linearly (conservatively) along the salinity gradient. Nonlinear patterns along the salinity gradient, on the other hand, often imply estuarine processes of transformation and nonconservative behavior (Boyle et al. 1974;Sharp et al. 1984).
Nonlinear patterns of conservative properties may also occur if end-member concentrations vary on time scales of estuary flushing (Loder and Reichard 1981), sometimes complicating the interpretation of property-salinity plots.
Property-salinity plots also have value independent of any role in deducing conservative behavior. They remain unchanged during simple movement of water parcels along the estuarine axis. Variability due to tidal movement is therefore decreased, allowing us to make more precise comparisons here among different seasons and eras. In addition, property-salinity plots facilitate comparisons among estuaries because they describe the estuarine axis in terms of a common salinity gradient rather than a multitude of estuarine sizes and shapes. Eyre (2000), for example, compared spatial patterns in nine Australian estuaries in this manner.
We illustrate the approach using measurements made in San Francisco Bay during pioneering studies of this estuary by the U.S. Geological Survey (Peterson et al. 1975a;Conomos et al. 1979). Seasonal sampling often revealed linear distributions of nutrients (nitrate 1 nitrite, silicate), phytoplankton biomass (chlorophyll a), and suspended particulate matter (SPM) along the salinity gradient during winter months (Fig. 1). During other seasons, however, the distributions of those constituents varied nonlinearly along the salinity gradient. For example, nutrient (N and Si) concentrations often had "concave" patterns with localized minima (Fig. 1A,B) while Chl a and SPM had "convex" patterns of localized maxima within the estuary (Fig. 1C,D). These kinds of pattern analyses were breakthroughs in the early development of estuarine science because they led to discoveries that biological processes (e.g., phytoplankton uptake) can strongly influence the spatial structure of reactive elements, and that physical processes associated with the salinity gradient concentrate phytoplankton and sediments in an estuarine turbidity maximum (ETM; Peterson et al. 1975b).
The examples in Fig. 1 illustrate that, as transitional ecosystems, estuaries are characterized by high spatial variability with spatial patterns that vary over time and differ among constituents. Thus, property-salinity plots have been instructive for interpreting or comparing results from individual sampling events. But they have not been combined for a unified and systematic depiction of how the salinity gradient structures spatial patterns of estuarine properties. Given the large variability of patterns illustrated in Fig. 1, we ask two fundamental questions about the functioning of estuaries as transitional systems between landscapes and seascapes: 1. What are the climatological spatial patterns of estuarine properties along the salinity gradient? (Here, we use "climatology" to refer to long-term means over the period of record.) 2. How and why do spatial patterns deviate from these climatological means?
We address these questions using long-term measurements of chemical, biological, and physical properties along the salinity gradient of San Francisco Bay as an example of a river-dominated tidal estuary. First, we describe the structure of the estuarine salinity gradient in this estuary and its seasonal and annual variability. Then, we identify and compare monthly mean distributions of six estuarine properties along the salinity gradient, selected to illustrate the diversity of spatial patterns that can develop in estuaries. Next, we show and discuss deviations from climatological patterns at the scale of events, years, and decades. We conclude with a discussion framed around these analyses to explain how estuaries are an ecosystem class distinct from the riverine and marine ecosystems they bridge.

San Francisco Bay as an estuarine ecosystem
The San Francisco Bay system includes two connected estuaries: South Bay is a marine lagoon situated in an urban landscape, and North Bay is the estuary of California's two largest rivers, the Sacramento and San Joaquin (Fig. 2). The rivers carry runoff produced in a 153,000-km 2 watershed bounded by coastal mountains and the Sierra Nevada Range and spanning the agricultural Central Valley. They deliver freshwater, sediments, nutrients, detritus, freshwater plankton, and agricultural chemicals to the estuary. The seaward boundary at the Golden Gate (Fig. 2) connects the estuary to shelf waters of the Pacific Ocean situated in the California Current System and strongly influenced by wind-driven coastal upwelling that seasonally alters oceanic water temperature, salinity, dissolved oxygen (DO), pH, nutrient concentrations, phytoplankton biomass, and primary and higher-level production (Raimonet and Cloern 2016). In this geographic setting, San Francisco Bay has been used as a model system to discover how ecosystems between land and ocean are structured, how they function, and how they change over time (Cloern and Jassby 2012).

The observational programs
The USGS began hydrographic and water quality studies of San Francisco Bay in 1969. Measurements have continued ever since, sustaining one of the world's longest observational programs in an estuary. We analyzed data collected in the northern estuary, defined here as the region between Sta. 18 near the Golden Gate and Sta. 657 in the lower Sacramento River (Fig. 2). Our primary analyses are of data collected since 1988 when we began using a Conductivity, Temperature, Depth sensor package (CTD) to measure highresolution vertical profiles of salinity and other constituents. Vertical profiles were obtained approximately monthly at 20 fixed stations along the salinity gradient. We also analyzed pre-1988 data to illustrate changes over time in the distributions of some properties. Sampling during that era was less frequent and included analyses of discrete water samples collected at 3-5 depths per station.
Methods have evolved over time as new sensor technologies and widely-accepted standard methods emerged (Cloern and Schraga 2016). In recent years, we measured vertical profiles with a Sea-Bird Electronics-9plus data acquisition system and sensors to measure: depth (Paroscientific Digiquartz pressure transducer), conductivity (SBE-4C conductivity sensor), temperature (SBE-3plus thermistor), SPM (Campbell Scientific OBS-3 optical backscatter sensor), Chl a (Turner Designs C7 fluorometer), photosynthetically active radiation (PAR; Li-Cor 192s quantum sensor), and DO (SBE-43 oxygen electrode). The optical backscatter sensor, fluorometer, and oxygen electrode were calibrated each sampling cruise with discrete measurements of ca. 20 samples collected in surface and near-bottom waters. SPM was measured gravimetrically as mass collected onto pre-weighed 0.4-lm pore size polycarbonate membrane filters, using a correction for salt retained on filters (Hager 1993). Chl a samples were collected onto GF/F filters, extracted in 90% acetone and analyzed using a Turner Designs fluorometer calibrated with Chl a standard. The acetone extracts were then acidified, and Chl a and phaeopigment concentrations calculated using equations of Jeffrey et al. (2005). Winkler reagents (Strickland and Parsons 1972) were added immediately to DO samples collected into Biological Oxygen Demand (BOD) bottles. Acidified DO samples were titrated with 0.01 N sodium thiosulfate dispensed by a Metrohm autotitrator (Titrino 798) using the potentiometric titration method of Gran eli and Gran eli (1991). Samples for dissolved inorganic nutrient analyses were filtered through 0.4-lm pore size polycarbonate filters under vacuum (less than 14 kPa), and the filtrates were frozen until analyzed by the USGS National Water Quality Laboratory. Concentrations were determined on an Aquakem 600 automated discrete analyzer (Thermo Scientific, U.S.A.) following Fishman and Friedman (1989) for nitrite, phosphate, and silicate, and Patton and Kryskalla (2003) for nitrate. All data and supporting metadata are available online at https://doi.org/10.5066/ F7TQ5ZPR.
From these measurements, we derived four additional estuarine properties. The light attenuation coefficient k (m 21 ) was calculated from vertical profiles of PAR measured with the quantum sensor. Stratification was computed from each vertical profile as the difference between sigma-t measured at 10 m and 1 m. We used this approach for consistency across sampling sites and because the pycnocline (median depth 4.5 m) was usually within the upper 10 m. Mixing diagrams shown below were essentially unchanged when stratification was based on density profiles of the entire water column. For comparison, we computed density stratification in Lake Tahoe and at station ALOHA in the N Pacific from profiles of temperature and salinity (assumed zero at Lake Tahoe), using function sw_den in R package marelac (Soetaert et al. 2016). Gross primary productivity (GPP, mg C m 22 d 21 ) was estimated with an empirical model derived from 14 C uptake assays to measure photosynthesis in San Francisco Bay (Cloern et al. 2007): GPP53:77 Á E Á C=k, where E is mean incident PAR (moles quanta m 22 d 21 ) for the month of sampling based on measurements with a LiCor 190s quantum sensor from 1979 through 2004, and C is Chl a (lg L 21 ).
We also used data collected by the California Department of Fish and Wildlife San Francisco Bay Study of the Interagency Ecological Program for the San Francisco Estuary. This program began in January 1980 and samples fish and macrocrustaceans at 42 sites using a midwater trawl for pelagic species and an otter trawl for demersal species. It also measures salinity and temperature with an SBE CTD. We examined spatial patterns along the salinity gradient for a pelagic marine fish, northern anchovy (Engraulis mordax), and three shrimp species: Crangon franciscorum, Crangon nigricauda, and Exopalaemon modestus (the latter introduced in the late 1990s). Abundances (as catch per unit effort, CPUE) of northern anchovy were measured through December 2014 and of shrimp through December 2013. Methods are described more fully by Baxter et al. (1999). These data were provided by Kathy Hieb, California Department of Fish and Wildlife (pers. comm., 18 August 2015).

Data analyses
We analyzed a data set that included 5429 vertical profiles of salinity and temperature, 5276 Chl a profiles, 4249 DO profiles, and 2064 measurements of near-surface dissolved inorganic nutrient concentrations from sampling between 28 February 1988 and 16 December 2015. From the vertical profiles, we used mean values in the upper 10 m, the depth of the shallowest station. Monthly climatology was identified by grouping all samples into 11 salinity bins (0-3, 3-6, . . ., 30-33) and then computing long-term mean values for each salinity bin and month. We did the same analyses for northern anchovy and shrimp abundances. Smooth lines in plots are local regression fits (Cleveland et al. 1992) using the loess function in R (R Core Team 2014) with span 5 0.7 and degree 5 1, i.e., locally linear fits. We use them only as a visual aid, not for any formal hypothesis testing.

Structure and its variability
The salinity gradient is a spatial pattern in the relative proportions of freshwater and seawater along the transition from river to ocean. We illustrate examples from San Francisco Bay for two hydrologic conditions-low freshwater inflow during the dry summer-autumn, and high inflow during a large winter storm (Fig. 3). These contrasting patterns illustrate how the structure of the salinity gradient varies with river inflow and is, therefore, tightly linked to variability of precipitation and runoff in the watershed. During the dry summer-autumn, the longitudinal salinity gradient extends along the full length of the estuary and salt penetrates > 80 km inland. During winter floods, the salinity gradient is compressed and displaced seaward.
The contrasting salinity distributions in Fig. 3 reflect changes in water circulation and mixing associated with variable river inflow. First, water density in San Francisco Bay is determined primarily by salinity (correlation coefficient 5 0.993), compared to lakes and oceans where density is controlled by temperature. The salinity range along San Francisco Bay (0-33) represents (at 158C) a density range of about 25 kg m 23 . This horizontal density gradient creates a pressure gradient that drives estuarine circulation-two-layer flow in which low-density surface water is advected seaward over a landward-flowing bottom layer (Monismith et al. 1996). Thirty to forty percent of the salt influx to San Francisco Bay is driven by estuarine circulation (Walters et al. 1985), so this is an important process of water movement and transport of solutes, particles, and organisms. The velocity of this estuarine transport process varies in direct proportion to the longitudinal salinity gradient, i.e., with river inflow (Geyer 2010). Second, estuarine circulation sets up stratification as a vertical salinity (density) gradient that can be stronger than the thermal stratification of lakes and oceans. As river inflow and strength of the longitudinal salinity gradient increase, so does the intensity of salinity stratification (Geyer 2010), as observed in San Francisco Bay (Fig. 3). Finally, river inflow also controls residence time. The water replacement time in San Francisco Bay is about 1 d during high-flow events, but about 2 months during the dry season (Walters et al. 1985). Physical dynamics of San Francisco Bay and other river-dominated estuaries, including structure of the salinity gradient, estuarine circulation, transport time scales, and stratification are driven by variability of freshwater inflow.

Three components of variability
The structure of the salinity gradient in San Francisco Bay can be parameterized with a single scalar, X2, defined as the distance (km) from the Golden Gate where bottom salinity is 2 (see Fig 3). The 2 isohaline was chosen because it demarcates the boundary between a fresh, vertically mixed region of net (tidally averaged) seaward flow from a salty, stratified region of estuarine circulation (Jassby et al. 1995). Seaward and landward bottom flows converge at X2, so its location can be a region of particle concentration including sediments and plankton (Fig. 1C,D). Several formulations have been developed to relate X2 to freshwater inflow Q (m 3 s 21 ), and we used the simplest based on the assumption of steady state (Monismith et al. 2002): X2 5 167ÁQ 20:141 . We used this equation to construct a 1988-2014 time series of monthly X2 using Q (Delta Outflow Index) computed by the California Department of Water Resources from measured river discharge and water exports and estimated local consumption (http:// www.water.ca.gov/dayflow/output/; accessed 17 June 2016).
The mean X2 over this period was 75.6 km but its position was highly variable, ranging from 48 km to 95 km (Fig.  4A). The series reveals variability at different time scales that we can compare. We decomposed the monthly X2 series into three components-annual, seasonal, and residual events-based on the additive model: X ij 5 X1y i 1m j 1r ij , where X ij is X2 for year i and month j; X is the long-term mean (75.6 km); y i is the annual effect; m j is the average seasonal (monthly) effect; and r ij is the residual series, which we sometimes refer to as the "events" component. The annual effect is simply the annual mean subtracted from the longterm mean. The monthly effect is given by the mean of that month's deviation from its corresponding annual mean. The events component is then obtained as the remainder. This simple additive decomposition mirrors the multiplicative decomposition described in more detail by Cloern and Jassby (2010). Both forms are motivated by the observation that many important events (e.g., persistent dry periods, species invasions) start or stop suddenly. Smoothing to extract the annualized term can disguise the timing of these events.
The annual component (Fig. 4B) tracks annual variability of precipitation. It reveals effects of multi-year droughts (1988)(1989)(1990)(1991)(1992)(2013)(2014)(2015) when annual mean X2 was displaced up to 10 km upstream from the long-term mean, and wet years (1995,1998,(2005)(2006) when X2 was displaced as much as 15 km seaward of the mean. Annual variability of precipitation is higher over California than other regions of the continental United States, and it is determined primarily by the number of large winter storms that provide much of the region's total annual precipitation (Dettinger 2016).
Multi-year droughts develop when large winter storms are absent; wet years have several large winter storms. Some of this annual variability is tied to the El Niño Southern Oscillation because extreme wet years, such as 1982-1983 and 1997-1998, were years of very strong El Niño, while the 1988-1989 drought occurred during a strong La Niña (http://ggweather.com/enso/oni.htm).
The average monthly pattern shows a 15-20 km seasonal excursion of X2 (Fig. 4C) that reflects California's Mediterranean climate of a cool wet winter and warm dry summer resulting from seasonal changes in the position and strength of the North Pacific High. This high pressure system has maximum intensity and northernmost position during summer, when it inhibits rainfall and blocks passage of storms formed in the Jet Stream (Iacobellis et al. 2016). It weakens and migrates to the southeast in winter, opening opportunities for storms formed over the North Pacific to pass over California.
The third component is the residual series ( Fig. 4D) representing event-scale deviations from the annual and seasonal components. These include large (> 20 km) downstream displacements of salt during exceptional flood events such as January and February 1997 (Figs. 3, 4) when X2 was 48 km and 53 km, respectively. These exceptional floods arise from landfall of atmospheric rivers-moisture-bearing air streams originating in the tropical or subtropical Pacific that can interact with storms along the Jet Stream to produce events of extremely high precipitation and flooding. These events are the source of 30-45% of California's annual precipitation (Iacobellis et al. 2016).
All components of the decomposition model are additive deviations from means with common units (km), so their importance can be compared directly. Standard deviations of each component were of comparable magnitude: 6.4 km, 6.4 km, and 4.7 km for the annual, monthly, and residual components (Fig. 4), corresponding to 40%, 39%, and 21% of the overall variance, respectively. Therefore, structure of the San Francisco Bay salinity gradient varies over three time scales, all of which must be considered to understand how spatial patterns of estuarine properties can be explained by spatially varying mixtures of fresh and seawater. Our strategy is to focus first on the seasonal component by describing monthly climatology (i.e., long-term means of spatial patterns by month), and then to explore annual and event-scale deviations from those mean patterns.

Monthly climatology
Contours of monthly mean salinity along the sampling transect reveal a clear annual cycle in the expansion and contraction of river and ocean influence on the estuary (Fig.  5). On average, river influence was strongest (isohalines displaced furthest seaward) in January-February, and ocean influence was strongest (isohalines displaced furthest landward) at the end of the dry season (August-October). The intervening months represent a transition from the wet to dry season as river inflow receded and isohalines migrated landward. The seasonal pattern of salinity distribution would look very different in estuaries situated in other hydroclimatic settings, such as the US Gulf or Atlantic coasts where rainfall occurs throughout the year and events of extreme high river inflow are associated with summer thunderstorms and tropical cyclones (Peierls et al. 2012) that develop during California's dry season.
The longitudinal salinity gradient provides a frame of reference for examining spatial patterns of estuarine properties along a seasonally-oscillating gradient of river-ocean influence (Fig. 5). It is particularly advantageous compared to the geographic frame of reference in tidal estuaries, such as San Francisco Bay, where water parcels are in continuous motion and displaced about 10 km up and down estuary twice daily by the semi-diurnal tide. We use this frame of reference to determine how estuarine properties vary along the mixing lines between freshwater and seawater, independent of geographic location. Climatological patterns along the salinity gradient

Diverse patterns
As an exploratory step, we examined contour plots of 1988-2015 mean values by month and salinity (Fig. 6). This revealed that some (not all) properties did vary along the salinity gradient, but the patterns took many different forms; we illustrate four. Silicate is an example of a property with highest concentration at the river end and lowest at the ocean end, regardless of season (Fig. 6A). Primary productivity followed the inverse pattern, consistently highest at the ocean end (Fig. 6B). Temperature varied along the salinity gradient, but the range of spatial variability for a given month was much smaller than the range of seasonal variability at a particular salinity (Fig. 6C). Stratification peaked at intermediate salinities in mid-estuary, and was strongest during the wet season (Fig. 6D). This diversity of patterns suggests that spatial-temporal distributions of estuarine properties can be structured by a variety of processes, some generating sharp spatial gradients and others stronger seasonal variability.
Next, we explored climatological patterns of six different properties to look for clues about processes that generate this diversity of spatial structure. We simplify the presentation by comparing mean spatial patterns for months representing three primary seasons and a range of mean X2 positions (Fig. 3): cool-wet (January), transitional spring (April), to warm-dry (August). Although nonlinear conservative mixing lines can be estimated when the sampling frequency is high enough (Eyre 2000), that is not an option here with our monthly historical data. Instead, we rely on existing studies of the estuary to support the presentation of underlying mechanisms rather than the shapes of property-salinity plots alone.

Solutes and particles Silicate
Weathering of silica-containing minerals in the watershed yields runoff with very high concentrations of dissolved silicate, on the order of 250-300 lM in the Sacramento River (Peterson et al. 1985). This is more than ten times greater than silicate concentrations in the adjacent coastal ocean, so San Francisco Bay is set up to have a large spatial gradient of silicate. The shape of that gradient can be modified by within-estuary processes such as diatom uptake (Peterson et al. 1975a) or dissolution of biogenic silica in sediments and release to the overlying water column (Grenz et al. 2000). Mean patterns of silicate for 1988-2015, however, showed linear relationships with salinity for every month (Fig. 7A). This is the classical pattern of conservative behavior, implying that within-estuary processes are either balanced or are slow relative to transit time through the estuary (Peterson et al. 1978). Lesson: the silicate example from San Francisco Bay exemplifies dissolved constituents having high concentration in a river source, low concentration in the ocean, and spatial distributions shaped by mixing between the two water sources. This conservative pattern is very different from the nonconservative patterns measured in the 1970s and 1980s (Fig. 1B), and we discuss this change in Si biogeochemistry later.

SPM
The Sacramento and San Joaquin rivers deliver an average 5700 tons of SPM, mostly fine-grained mineral particles, to the estuary annually (Schoellhamer et al. 2016). Most (82%) of that sediment influx comes in winter, and in particular with the first large discharge pulses (first flush) that carry extremely high SPM concentrations. Much of this sediment mass is deposited downstream in the seaward estuary and later, as freshwater inflow recedes, is transported up-estuary by the combination of wind-and tidally-driven suspension and density-driven circulation (Schoellhamer 2001).
These annual cycles of changing sediment supply and transport produce large seasonal changes in the spatial distribution of SPM that were not seen for silicate. The mean winter (January) pattern is a gradient of decreasing SPM concentration between the river (mean $ 80 mg L 21 ) and ocean ends ($ 20 mg L 21 ). As sediments are carried seaward  by high river flow they settle, producing a concave pattern in winter (Fig. 7B) diagnostic of an internal sink, in this case deposition. The mean spring (April) pattern is convex, traditionally interpreted as an internal source. However, in this case, the localized maximum SPM concentration results from particle accumulation in an ETM that forms in the lowsalinity zone of this estuary (Fig. 7B). ETMs are distinctive features of estuaries and their formation in the low-salinity zone of San Francisco Bay during spring is a result of several processes (Schoellhamer 2001). First, as river inflow recedes the density-driven landward bottom flows transport sediment further into the estuary. Second, winds intensify in spring to suspend the erodible sediment pool supplied by winter runoff, so SPM concentrations increase. Third, the low-salinity zone is positioned near a sill that blocks estuarine circulation and causes convergence to produce the ETM. During summer-autumn, the sediment supply, wind intensity, and strength of estuarine circulation all diminish, so SPM concentrations are reduced throughout the estuary (Fig. 7B).
Lessons: the SPM example exemplifies a constituent delivered by river flow that decreases along the salinity gradient but exhibits patterns that vary seasonally. The contrasting spatial patterns of silicate (linear, constant) and SPM (nonlinear, variable) illustrate the effect of processes that transport particles but not solutes-sinking, deposition, wind and tidal resuspension, and accumulation in ETMs. The SPM example also illustrates the effect of bottom topography on particle transport and spatial patterns in San Francisco Bay and other estuaries (Schoellhamer 2001). Location of the spring SPM maximum is set partly by the geographic location of a sill that controls upstream flow. Therefore, spatial variability of bottom topography can be as important as the salinity gradient in structuring spatial patterns of some estuarine properties. This includes an optical property, the light attenuation coefficient k (m 21 ), which is highly correlated with and follows the same seasonal-spatial patterns as SPM. Spatial structure of k has an important ecological consequence that we describe below.

Physical properties Temperature
Water temperature in San Francisco Bay is determined by heat fluxes across its interfaces with the coastal ocean, rivers, and overlying atmosphere (Monismith et al. 2009). We know, for example, that temperature variability in San Francisco Bay is strongly correlated with temperature variability in the adjacent coastal ocean, but the correlation weakens with distance into the estuary (Raimonet and Cloern 2016). We also know that temperature in the estuary responds rapidly to weather anomalies such as heat waves . However, the relative importance of heat fluxes by water exchange and atmospheric exchange is not known.
Monthly mean water temperature in San Francisco Bay followed linear or near-linear distributions along the salinity gradient each season (Fig. 7C), indicating that river-ocean mixing might play a role. The sign of the temperature gradient shifted seasonally from positive-seaward during winter to positive-landward during summer. This reflects the large heat capacity of the Pacific Ocean that dampens its seasonal temperature variability relative to that over land, leading to seasonal reversals of the atmospheric temperature gradient between ocean (warmer in winter, cooler in summer) and land (Iacobellis et al. 2016). The linear temperature patterns could thus also be explained by the along-estuary gradient of surface heat flux associated with the river-to-ocean gradients of atmospheric temperature and cloudiness.
Lessons: the temperature example illustrates how differing phenology (seasonal phasing, amplitude) of variability at river and ocean boundaries can alter spatial gradients within estuaries. The temperature patterns also illustrate that some estuarine properties have greater seasonal variability than spatial variability. This implies that their dynamics are driven primarily by processes other than river-ocean mixing or within-estuary transformations. We would expect this to be true for other properties whose variability is tied to atmospheric forcing, such as incident solar radiation, wind mixing, and gas exchange.

Stratification
Thermal stratification of lakes and oceans is a central feature of their physical dynamics because it dampens turbulent mixing and slows exchanges between surface and deep waters. The maximum density difference from thermal stratification, even in deep lakes and oceans, is on the order of 1-4 kg m 23 . We show examples from Lake Tahoe, where a 13.58C temperature gradient in the upper 500 m represented a density difference of 1.4 kg m 23 , and from the North Pacific at station ALOHA where a weak salinity gradient and 22.48C temperature gradient in the upper 1000 m represented a density difference of 4.2 kg m 23 (Fig. 8A,B). Water density in estuaries is controlled primarily by salinity, and the density difference between freshwater and seawater creates the potential for stronger density stratification than seen in lakes and oceans. We show an example from Sta. 18 during a winter flood, when a 1.58C temperature difference and a 20.6 salinity difference corresponded to a density difference of 15.7 kg m 23 over 39 m (Fig. 8C).
The climatological patterns in San Francisco Bay always showed weakest stratification (< 1 kg m 23 density difference between depths of 10 m and 1 m) near the river and ocean boundaries and strongest stratification within the estuary. Spatial structures were smooth, convex patterns of variability with peak stratification in the salinity range 10-25 (Fig. 7D). The amplitude of the spatial pattern varied seasonally, with strongest stratification during winter and weakest during summer. This seasonal pattern is well understood because stratification of estuaries varies in proportion to freshwater inflow as a source of buoyancy (Geyer 2010). Several factors generate the spatial pattern of strongest stratification at intermediate salinities in San Francisco Bay. Salinity stratification requires a supply of both freshwater and seawater, and the maximum potential stratification would occur when freshwater overlays a seawater layer of equal thickness so the depth-averaged salinity is $ 17, i.e., within the estuary. In addition, the longitudinal salinity gradient that sets up vertical stratification in San Francisco Bay is strongest in the midsalinity range (Monismith et al. 2002). This combination of factors creates a particular spatial-seasonal pattern that was unique among the properties we investigated: strongest stratification at intermediate salinities during winter (Fig. 6).
Lessons: a distinctive feature of estuaries is the potential for strong vertical stratification arising from the salinity difference between fresh and seawater. Stratification in estuaries varies with the seasonal pattern of freshwater inflow, distinct from the seasonal patterns of thermal stratification in lakes and oceans tied to the seasonal cycle of surface heat flux. Stratification is not uniform along the salinity gradient but, instead, is strongest where mean salinity is around the average between freshwater and seawater. The applicability of these lessons to other estuaries probably depends on their bottom topography (through its control of circulation) and hydroclimatic setting.

Biota and biological processes Primary productivity
Estuaries are often described as high-productivity ecosystems because they and their connected marshes and mudflats have a diverse and abundant community of algal and vascular-plant producers. However, the phytoplankton component of annual primary production ranges widely, from < 10 g C m 22 to > 1000 g C m 22 . Many estuaries are enriched in nutrients. For example, 1988-2015 mean concentrations of phosphate (2.4 lM), dissolved inorganic nitrogen (32.4 lM), and silicate (178 lM) in northern San Francisco Bay exceed levels that limit algal growth. As a result, primary productivity in this and other estuaries is often regulated by high turbidity and associated light limitation of photosynthesis (Cloern 1987).
Monthly mean GPP in San Francisco Bay was always higher at the ocean end than at the river end (Fig. 7E). The spatial patterns showed a monotonic increase of productivity along the salinity gradient. The amplitude of that spatial gradient varied seasonally, ranging from $ 20-100 mg C m 22 d 21 during winter to $ 200-1000 mg C m 22 d 21 during spring and summer. These productivity estimates were derived from three quantities: monthly mean solar radiation, Chl a, and turbidity indexed as the light attenuation coefficient k.
San Francisco Bay is a turbid estuary because of its high SPM concentrations-the Spearman rank correlation coefficient between climatological mean k and SPM was 0.93 (n 5 132)-and so the spatial patterns of turbidity (k) mirror the spatial patterns of SPM (Fig. 7B). The general pattern of increasing primary productivity from river to ocean is thus partly explained by the pattern of decreasing k and increasing photic depth toward the ocean. For example, mean August k was 2.52 m 21 in the lowest salinity bin and 0.83 m 21 at highest salinity. This difference corresponds to an estuarine gradient of deepening photic depth from 1.8 m . The station ALOHA data were obtained from Hawaii Ocean Time-series data report 1 (Chiswell et al. 1990). The San Francisco Bay data were obtained from the USGS San Francisco Bay water-quality website (http://sfbay.wr.usgs.gov/access/wqdata/; accessed 11 August 2016).
to 5.6 m (assuming photic depth is the depth of 1% surface irradiance). Chl a concentrations also generally increased along the salinity gradient and, for every month, mean values were highest at the ocean boundary. Therefore, spatial structure of primary productivity in this estuary (Fig. 7E) arises from the compounding effects of increasing phytoplankton biomass and light availability along the salinity gradient.
Lessons: the primary productivity example illustrates sharp spatial gradients that can develop in estuaries. During summer, mean primary productivity varied five-fold along a 100-km transect. Much of that variability is associated with spatial gradients of river-derived sediments and their attenuation of light, a feature not present in the open ocean or beyond the littoral zone of lakes. The general pattern of increasing productivity from river to ocean occurs in many other estuaries, fjords, and bays influenced by land runoff including the Ems-Dollard, Corpus Christi Bay, Howe Sound, Delaware Bay, Bristol Channel, and Chesapeake Bay (Cloern 1987).

Biological communities
Aquatic microbes, plants, and animals have evolved physiological and life-cycle adaptations to salinity variability. Few organisms can osmoregulate or complete their life cycles along the full gradient from fresh to seawater, but many are adapted to a salinity range that is taxon-specific. Therefore, biological communities are structured by the estuarine salinity gradient. We show an example as mean abundances of three shrimp species along the salinity gradient during August (Fig. 9). C. franciscorum (California bay shrimp) and C. nigricauda (blacktail bay shrimp) are marine species that migrate into estuaries as juveniles or reproductive adults. C. nigricauda has an optimal salinity range of 18-32 (Unger 1994) and its abundance in the estuary was highest at salinity > 20. C. franciscorum has an optimal salinity range of 2-22 (Unger 1994) and its population was centered in the salinity range 5-20. E. modestus (Siberian prawn) is a freshwater species found only in the upper estuary and tidal rivers (Fig. 9). Therefore, three species of one community partition estuarine habitats by their adaptations to different salinity ranges. Habitat partitioning also occurs within species because early-stage juveniles of the Crangon species occupy low salinities of their habitat range and then migrate seaward to occupy higher salinities as they mature (Hatfield 1985).
Other biological communities are similarly structured by the salinity gradient. The fish community of San Francisco Bay includes 137 species that group into guilds based on their salinity ranges (Feyrer et al. 2015), including: marine species found only at high salinity (e.g., yellowtail rockfish Sebastes flavidus, Dover sole Microstomus pacificus, white seabass Atractoscion nobilis); freshwater species confined to low salinity (e.g., juvenile Western brook lamprey Lampetra richardsoni, spotted bass Micropterus punctulatus, tule perch Hysterocarpus traskii); and generalists found across a wide salinity range (e.g., adult American shad Alosa sapidissima, starry flounder Platichthys stellatus, longfin smelt Spirinchus thaleichthys). The phytoplankton community also includes freshwater taxa found only at salinity < 5 (Anabaena affinis, Chlorella vulgaris, Cryptomonas ovata, Fragilaria crotonensis), marine taxa found only at salinity > 25 (Ceratium furca, Biddulphia pulchella, Pyramimonas parkeaea), and others found across a wide salinity range (1-25) such as Gyrosigma acuminatum, Plagioselmis prolonga, and Cyclotella striata (Cloern and Dufford 2005). Salinity is the primary control on the distribution of marsh plants in San Francisco Bay (Watson and Byrne 2009) and other estuaries (Crain et al. 2004). Spatial patterns of microbial communities, such as ammoniaoxidizing bacteria and archaea, are also structured by the salinity gradient (Mosier and Francis 2008).
Lessons: these examples illustrate that the estuarine salinity gradient is a habitat gradient occupied by organisms having species-specific adaptions to salinity variability. The salinity ranges of some species are narrow, so biological communities have the same sharp spatial gradients of other properties, changing over distances of a few tens of kilometers. The geographic location of motile estuarine species (plankton, crustaceans, fish) tracks the salinity gradient as it oscillates with tidal currents and variable river inflow (Fig.  5).

Deviations from the climatological patterns
Climatological patterns are instructive because they reveal mean spatial structures along estuarine salinity gradients and provide clues about their underlying processes. They also provide a benchmark for detecting and measuring variability of estuarine spatial patterns over time. We show three examples as deviations from the 1988-2015 climatological patterns at time scales of decades, years, and events.

An ecosystem state change
The northern San Francisco Bay ecosystem was transformed after the 1986 introduction of Potamocorbula amurensis, a small clam indigenous to Asian rivers and estuaries that, unlike marine mollusks, survives across the full salinity range of San Francisco Bay. Its population exploded, reaching mean abundance of 4000 m 22 within a year of its introduction (Nichols et al. 1990). At that density, the rate of P. amurensis filter feeding exceeds the phytoplankton growth rate (Kimmerer and Thompson 2014), and the most immediate effect of this species introduction was a five-fold reduction of phytoplankton biomass and primary production in low-salinity regions of the estuary previously uninhabited by bivalves (Alpine and Cloern 1992). This large-scale, unintended experiment provides an opportunity to learn how the spatial structures of estuarine properties can be modified by an ecosystem-scale perturbation. To do this, we compared mean distributions of four properties over the eras before and after the 1986 P. amurensis introduction.
Early studies of San Francisco Bay discovered a summer phytoplankton bloom that developed each year over a period of months as a consequence of slow flushing, high productivity in shallow bays, and biomass accumulation in the ETM (Cloern et al. 1983). In that pre-clam era, a broad phytoplankton biomass peak developed in the salinity range 3-14 where Chl a concentration approached or exceeded 10 lg L 21 . That summer pattern disappeared in the post-clam era, and biomass has been uniformly low (mean Chl a 5 2.6 lg L 21 ) along the full salinity gradient (Fig. 10A). The summer blooms were dominated by diatoms, and mean silicate concentrations followed the classical concave pattern resulting from silicate uptake within the estuary (Fig. 10B). Using deviations along the river-ocean mixing line, Kimmerer (2005) estimated that silicate uptake rate averaged 2.6 lmol Si L 21 d 21 during the pre-clam era but only 0.4 lmol Si L 21 d 21 after the summer blooms disappeared. As a result, the silicate spatial pattern now follows the conservative mixing line between rivers and ocean (Fig. 10B). The spatial pattern of nitrate also changed after 1986. Summer nitrate concentrations during the pre-clam era followed a pattern of increase from river to ocean, with minimum concentrations aligned with the Chl a maximum (Fig. 10C). Nitrate concentrations in the post-clam era increased by 10-20 lM along the entire salinity range, and the spatial pattern shifted to a mid-estuary maximum and lowest concentrations at the river and ocean boundaries (Fig. 10C). This new spatial pattern reflects two processes of change: (1) reduced phytoplankton uptake after the clam invasion, and (2) increased nitrate inputs from population growth in the surrounding landscape, including discharges from 24 municipal wastewater treatment plants distributed along northern San Francisco Bay (Novick and Senn 2014).
Loss of the summer bloom also had major ecological effects because phytoplankton primary production is the largest source of organic matter supporting metazoan food webs in the San Francisco Bay-Delta (Sobczak et al. 2005). The reduction of mean summer Chl a concentrations from ! 10 lg L 21 to < 3 lg L 21 (Fig. 10A) represented a shift to a state of chronic food limitation of consumers such as calanoid copepods and cladocerans (Mueller-Solger et al. 2002;Kimmerer et al. 2005). As a result, summer abundances of key components of the zooplankton community-the rotifer Synchaeta bicornis, copepods Eurytemora affinis and Acartia spp., and mysid shrimp Neomysis mercedis-all decreased abruptly after the summer blooms disappeared in 1987 (Kimmerer 2002).
Higher-level consumers can adapt to changes in their food supply by altering their spatial distribution within estuaries. We show an example, comparing pre-and post-clam distributions of northern anchovy, the biomass-dominant planktivore in San Francisco Bay. Prior to 1987, the summer distribution of northern anchovy was centered in the salinity range from about 10-23 (Fig. 10D). After 1986, the population shifted seaward so that maximum abundance is in the highest salinity bins. Kimmerer (2006) first reported this shift, noting that summer abundance of northern anchovy decreased 94% in low salinity habitats after the clam introduction. He interpreted the population decline and altered spatial distribution as behavioral accommodation of northern anchovy to the reduced supply of phytoplankton and zooplankton food resources. Other kinds of animal behaviors, such as vertical migrations of copepods in oscillating tidal flows, can retain smaller organisms and control their spatial distributions within estuaries (Kimmerer et al. 1998).
Lessons: this example illustrates that spatial patterns of estuarine organisms and biogeochemical processes can be restructured by ecosystem disturbances. In this case, a disturbance (increased grazing) altered an internal estuarine process (nutrient uptake), leading to system-wide changes as elevated concentrations and altered spatial patterns of silicate and nitrate along the entire salinity gradient. Disturbance effects captured by long-term observations can reveal processes that shape spatial patterns. The post-clam spatial pattern of nitrate has lowest concentrations near the river and ocean sources, reflecting inputs from wastewater treatment plants distributed along the estuary. Therefore, spatial patterns in estuaries can be shaped by sources other than riverine and oceanic inputs. The altered spatial distributions of northern anchovy illustrate a very different process: organisms can adjust the spatial structure of their populations through behavioral adaptations to changing environments, such as altered food supplies.

A 2-yr weather anomaly
Our second example illustrates deviations from mean patterns caused by annual weather anomalies. Global mean temperature reached record highs in 2014 and then in 2015 when mean temperatures over land and over oceans exceeded long term means in the northern hemisphere by 1.338C and 0.748C, respectively (NOAA National Centers for Environmental Information 2016). During these two years, a "highly anomalous weather pattern" developed in the NE Pacific Ocean as a strong high pressure ridge persisted and an unprecedented mass ("the blob") of warm water expanded along the west coast of North America (Bond et al. 2015). Manifestations of this event, "the largest marine heat wave ever recorded" in the NE Pacific (Di Lorenzo and Mantua 2016), were expressed inside San Francisco Bay. We measured record-high water temperatures in either 2014 or 2015 at every sampling location for each month, and the 2014-2015 mean temperatures exceeded the climatological means in each salinity bin by 0.5-38C (Fig. 11). Positive temperature anomalies occurred along the full salinity gradient, suggesting that warming of the estuary was a response to warming of both coastal waters and air temperature as effects of the NE Pacific weather anomaly propagated landward (Bond et al. 2015). Lesson: this example illustrates how estuaries respond to global-and regional-scale weather anomalies and that those responses can be system-wide and expressed along the entire salinity gradient.

A weather event
Deviations from mean spatial patterns can also occur as responses to short-term events, and we illustrate an event of anomalously low DO. DO concentrations in the upper 10 m of San Francisco Bay averaged 91% saturation over the period 1988-2015, and mean DO varied weakly along the salinity gradient (Fig. 12A). Under-saturation of DO is a characteristic of estuaries as net-heterotrophic ecosystems (Caffrey 2004). But the absence of spatial structure in DO is unexpected, given the strong spatial pattern of primary productivity (Fig. 7E). Maintenance of mean DO around 90% saturation indicates an overall equilibrium between mean rates of photosynthesis, respiration, and atmospheric exchange. Any effects resulting from spatial variability of metabolism are apparently buffered by atmospheric exchange in this estuary. However, large  deviations around this climatological pattern can occur, and we show an example from 09 May 2006 when mean DO in the upper 10 m decreased significantly from river to ocean and was only 65% saturation near the ocean boundary (Fig.  12A). Water temperature was also anomalously low during this event, particularly near the ocean boundary (Fig. 12B). The combination of low DO and low temperature suggests that this estuarine anomaly resulted from an intrusion of deep ocean water brought to the surface by an event of winddriven coastal upwelling (Raimonet and Cloern 2016). Lessons: the DO example illustrates that processes, such as atmospheric gas exchange, can buffer variability of estuarine processes to dampen spatial structure that might otherwise develop. It also illustrates how event-scale ocean variability can propagate into estuaries to alter their spatial patterns.

Conclusions
A central problem of ecology is to understand how spatial patterns of habitats, biota, and biogeochemical processes are structured. The transitional zone between rivers and oceans is a gradient of salinity that develops as land runoff encounters and mixes with seawater. The salinity gradient provides a framework for characterizing and understanding spatial patterns in estuaries and their variability over time. Lessons illustrated with long-term measurements along the salinity gradient of San Francisco Bay provide a foundation for developing and testing fundamental principles of how estuaries are structured and how they function as ecosystems. These lessons include: 1. The structure of the salinity gradient varies with river inflow, so spatial and temporal patterns of estuarine ecosystem variability are strongly tied to the hydroclimatic setting. 2. Dissolved constituents having different concentrations in rivers and oceans (e.g., silicate) can vary linearly along the salinity gradient, a pattern shaped by passive mixing between freshwater and seawater. 3. Dissolved constituents can also have nonlinear patterns (e.g., silicate during diatom blooms) shaped by the combination of river-ocean mixing and within-estuary processes such as phytoplankton uptake. However, nonlinear patterns can also arise from departures from steady state, so their interpretation requires measurements of variability at estuarine boundaries. 4. Particles accumulate in estuarine turbidity maxima formed by the combination of particle sinking and estuarine circulation, so they can have different spatial patterns than dissolved constituents. 5. The location of estuarine turbidity maxima can be shaped by bottom topography and its control of densitydriven currents and particle transport.
6. Some constituents have weak or no mean spatial structure along the salinity gradient, reflecting spatially distributed sources along the estuary (nitrate) or atmospheric exchanges that buffer spatial variability of ecosystem metabolism (DO). 7. Turbid, nutrient-rich estuaries commonly have a spatial pattern of increasing primary productivity along the salinity gradient as river-derived sediments sink and the photic zone deepens. 8. The density difference between freshwater and seawater creates the potential for density stratification in estuaries stronger than the thermal stratification of lakes and oceans. Stratification in San Francisco Bay has a consistent spatial-seasonal pattern, being strongest around the center of the salinity gradient and when river discharge is high. 9. Spatial distributions of motile organisms are shaped by species-specific adaptations to salinity variability (shrimp) and by behavioral responses to other components of environmental variability such as food supply (northern anchovy). 10. Estuarine spatial patterns can change over time scales of events (intrusions of upwelled ocean water), seasons (river inflow), and years (annual weather anomalies), and between eras separated by ecosystem disturbances (a species introduction).
The classical definition of an estuary is "a semi-enclosed coastal body of water which has a free connection with the open sea and within which seawater is measurably diluted with freshwater derived from land drainage" (Cameron and Pritchard 1963). Collectively, the results presented here illustrate how estuaries are more than mixing zones between fresh and seawater. Rather than just ecotones or ecoclines along the land-sea transition, estuaries are a separate class of ecosystem having complex and diverse spatial patterns, sharp spatial gradients, unique structuring processes, and biogeochemical patterns and biological communities distinct from those in lakes, rivers, and oceans.