Coastal eutrophication is a prevalent threat to the healthy functioning of ecosystems globally. While degraded water quality can be detected by monitoring oxygen, nutrient concentrations, and algal abundance, establishing regulatory guidelines is complicated by a lack of baseline data (e.g., pre-Anthropocene). We use historical carbon and nitrogen isoscapes from sediment cores to reconstruct spatial and temporal changes in nutrient dynamics for a central California estuary, where development and agriculture dramatically enhanced nutrient inputs over the past century. We found strong contrasts between current sediment stable isotopes and those from the recent past, demonstrating shifts exceeding those in previously studied eutrophic estuaries and substantial increases in nutrient inputs. Comparisons of contemporary with historical isoscapes also revealed that nitrogen sources shifted from a marine-terrestrial gradient to amplified denitrification at the head and mouth of the estuary. Geospatial analysis of historical data suggests that an increase in fertilizer application – rather than population growth or increases in the extent of cultivated land – is chiefly responsible for increasing nutrient loads during the 20^th century. This study demonstrates the ability of isotopic and stoichiometric maps to provide important perspectives on long-term shifts and spatial patterns of nutrients that can be used to improve management of nutrient pollution.

To examine interrelationships between nitrogen pollution and anthropogenic sources over the past century, we parametrized a model of nitrogen inputs to the watershed. Our model was based on the Nitrogen Loading Model (NLM). We applied the NLM model to calculate watershed sources of nitrogen over time in decadal increments from 1930–2010. We compiled historical data on changes in human population from census data, atmospheric deposition, homes with wastewater treatment, the areal extent of cultivated and natural lands and impervious surface cover, and estimated changes in fertilizer application rates in the Elkhorn watershed (based on annual "Commercial Fertilizers" and "Fertilizing Materials" reports published by the California Department of Agriculture 1925–2012).

Eighty-five ~ three-meter-deep sediment cores were collected during 2010 from the vertices of a 200 m x 200 m grid superimposed over the tidal and never-diked portions of the estuary. Most of the sediment cores were collected using a Russian peat borer to minimize compaction; in a few locations, a vibracorer was necessary to penetrate sands. Six focal cores were selected for high-resolution analyses and were collected using a piston corer with polycarbonate liners to obtain intact core sections for scanning and archiving. Focal cores were split into 1-cm sections; the remaining cores were sectioned into 10 cm intervals for 0–50 cm depths, and into 25-cm intervals for 50–100 cm depths. Core splits were archived at the LacCore repository at the University of Minnesota. Chronologies were created using downcore profiles of ²¹⁰Pb, ¹³⁷Cs, and ²²⁶Ra measured with a low-energy germanium multichannel gamma spectrometer. Historical geochemical markers included Pb concentrations measured using ICP-AES following four-acid extractions, AMS radiocarbon dating of fossil peat, and magnetic susceptibility and imaging using a Geotek Multi-Sensor core logger. The maximum depth of radiocesium was assigned an age of 1953, radiocesium peaks were assigned an age of 1963, and total lead concentration peaks were assigned an age of 1974.

Lead-210, radiocesium, and radiocarbon dating were combined in an age-depth model using a Bayesian approach to construct chronologies for seven cores. The age model ²¹⁰Pb Plum in R version 4.0.5 uses the same statistical approach as the previous model Bacon, but incorporates radionuclide dating including parameters of deposition of ²¹⁰Pb, supported ²¹⁰Pb, and accretion rates. The Plum model was selected because it can account for incremental ²¹⁰Pb data over depth in the cores, as opposed to using the analytical approach of the continuous rate of supply model. Additionally, this model has been used previously for chronologies of estuarine sediments. Within Elkhorn Slough, sediment accumulation rates varied little from site to site over the past century and were similar to values reported previously; thus, to estimate ages for the 85 undated cores, we compiled a composite core chronology using the seven cores to represent mean age-date model for the entire estuary. This composite core chronology was then applied to the 85 undated cores, using the composite age-depth relationship to estimate dates for the depth segments utilized for isotopic and stoichiometric measurements. We report the mean year output of the model and 95% confidence intervals around the mean.

For the six high-resolution focal sites, cores were analyzed at 1-cm increments (for 0 to 50 cm depths) for stable carbon and nitrogen isotopic composition using a Finnegan Delta Plus continuous flow isotope ratio mass spectrometer using standard methods, and for carbon and nitrogen concentration using a Flash 1112 EA. For the 85 coarser resolution cores, sediments were analyzed for carbon and nitrogen abundance and stable isotope ratios using a Vario Cube elemental analyzer interfaced to an Isoprime 100 IRMS. Isotope ratios for carbon and nitrogen are reported in permille notation as:  where R is the abundance ratio of the less common (a) to more common isotope. The standard for nitrogen is atmospheric nitrogen gas; the standard for carbon is PeeDee Belemnite; by definition, standards have δ=0. Sediments were not pretreated to remove inorganic carbon, as acidification did not quantitatively shift ratios. Previous studies suggest little effect of diageneses on sediment δ¹⁵N ratios in coastal marine settings, but shifts of ~ -1.5‰ in δ¹³C ratios are expected and C/N ratios are thought to decrease over time. Furthermore, atmospheric δ¹³C ratios have declined by about -1.5‰ since 1850 associated with the Suess Effect – the release of lighter C from fossil fuel combustion.

Whole estuary isoscape and stoichioscape maps were produced using sedimentary stable isotope (δ¹³C and δ¹⁵N) and molar nutrient stoichiometric (C/N) ratios interpolated from the 85 core locations using ordinary kriging in ArcGIS version 10.2.2 (ESRI, Redlands, CA, USA) to the spatial extent of cored areas in Elkhorn Slough. Maps were created for six depth intervals dated using the composite chronology (ca. 1726–1839, 1839–1885, 1885–1951, 1951–1963,1963–1981, and 1981–2010). Different interpolation variogram models including spherical, circular, exponential, Gaussian, linear interpolation with linear drift, and linear with quadratic drift were tested. Leave-one-out cross validation of 15% of the points was used to choose the model which yielded the smallest root mean square error between predicted and actual values. To ensure that historical differences in interpolation maps were a function of data differences rather than variogram methodology, the spherical kriging method was used for all timepoints. We also applied data from monthly water quality sampling at a network of (~26) stations across Elkhorn Slough since 1988. Monthly nitrate data from the sites were averaged during the full year of 1995 and mapped using ordinary kriging for comparison to spatial patterns of the isoscape and stoichioscape maps.

Trends in isotopic and stoichiometric signatures since the 1850s were examined for the six high-resolution cores. Timeseries analysis of the high-resolution data investigated the statistical significance of trends during the period of increasing fertilizer application, as well as offsets in the signatures associated with the timing of marine inlet construction for the harbor. Statistically significant change points in the timeseries were determined using the Pettitt Test, a nonparametric test that identifies the year of a step change and assigns significance to the selection. Datasets of δ¹⁵N, δ¹³C, and C/N for each of the six high-resolution coring sites were separately tested for the period 1850–2010 (n = 45 time points each). Next, the timeseries were split at the significant step change points that were statistically identified, forming two datasets “before” and “after” the year of change. Trend analysis was performed using linear regression on the split datasets, to model the slope after the split as well as the difference of y-intercept at the step change year (Fig. S1 diagrams the slope and intercept of our statistical models). The difference of y-intercept at the step change year is interpreted as an offset in the timeseries, consistent with construction of the harbor inlet when the step occurred at the same time as the construction (1946 ± 10 years). The slope after this step change year is attributed to increasing fertilizer addition to the watershed from 1940–1980.

To compare sediment isotope results to dissolved nutrient concentrations, we compared water quality monitoring data to the high-resolution sediment cores during a 20-year period. Monthly water sampling of parameters (including salinity) were measured at the sites, and water samples were also collected into brown Nalgene bottles; stored on ice; filtered; and analyzed for nutrients, including nitrate (NO₃⁻), within 48 hours, or frozen for later analysis in accordance with standard methods. Three of the high-resolution sediment cores were collected at the same locations as water quality monitoring sites. For these three water quality sampling stations (Portero Road North, Kirby Park, and Hudsons Landing West), we compared annual mean water column dissolved NO₃⁻ (mM) and salinity (ppt) to sedimentary δ¹⁵N values during the same year from 1990–2010.

Excel, R, and a GIS software such as ArcGIS or QGIS.