Hindcasting estuary ecological states using sediment cores, modelled historic nutrient loads, and a Bayesian network

1 Introduction

Estuaries are among earth’s most anthropogenically impacted ecosystems and threatened habitats (Kennish, 2002; Barbier et al., 2011). Globally, estuarine ecological degradation has accelerated over the past 150 to 300 years (Lotze et al., 2006) as a result of anthropogenic overexploitation, habitat transformation and pollution (Jackson et al., 2001; Lotze and Milewski, 2004; Lotze et al., 2005), particularly from increased diffuse-source catchment inputs of sediment and nutrients (Conley, 1999; Le Moal et al., 2019; Vieillard et al., 2020). The effects of these contaminant load increases include eutrophication (phytoplankton and macroalgal blooms; Barr et al., 2020), and increased sedimentation and organic content in sediments (Engelsen et al., 2008), with attendant ecological disruption including increased light attenuation, hypoxia and sulfide toxicity (Howarth and Marino, 2006; Fowler et al., 2013). These stressors can result in loss of important macrofaunal communities and habitats including seagrass beds and their ecosystem services (Green et al., 2014; Pratt et al., 2014; Sutula et al., 2014).

Estuaries in Aotearoa New Zealand have been increasingly impacted by catchment land use changes causing excessive nutrient and sediment inputs in freshwater flows (Lohrer et al., 2006; Morrison et al., 2009; Handley et al., 2020), especially since European colonization (Zeldis et al., 2010; Snelder et al., 2018). Increased sedimentation due to sediment losses from land has been attributed to deforestation for agriculture (Hume et al., 1992), and agricultural land uses and urbanization have increased nutrient inputs to freshwater (Larned et al., 2016) and estuaries (Dudley et al., 2020). Major changes in landcover in New Zealand resulted from large-scale catchment deforestation and conversion to pastoral agriculture during the late 19^th and early 20^th centuries. Technological advances, such as fertilizer top-dressing, increased agricultural production in the mid-20^th century, and increases in irrigation, stock densities and fertilizer application to pasture after the mid-20^th century (Ministry for the Environment, 2007; Howard-Williams et al., 2010; Ministry for the Environment and Stats NZ, 2019) allowing the intensification of agricultural land use and an increase in the area of urban land cover and human populations have resulted in increases of estuarine contaminant loads.

To manage these impacts, there are current New Zealand National Planning Framework (NPF) requirements that local governments develop policy statements (e.g. the National Policy Statement for Freshwater Management; Ministry for the Environment, 2023) that consider receiving waters (including estuaries) when setting freshwater and marine environmental limits and targets for restoration. Toward this, there is demand for management tools that quantify relationships between nutrient and sediment loads and estuary ecological condition (Townsend and Lohrer, 2015; Plew et al., 2020). Recent work in New Zealand has delivered a series of tools to link contaminant loads to estuarine trophic state (Plew et al., 2018, 2020; Zeldis and Plew, 2022), although these have only been applied in the context of modern estuarine conditions. Furthermore, quantitative monitoring data describing the progression and impact of nutrient/sediment stressors in New Zealand’s estuaries was usually begun only within the last two decades, often after degradation of estuary condition had commenced or was well advanced (Hewitt et al., 2014; Dudley et al., 2020). This has made it difficult to ascertain historical estuary condition and trajectories of change.

Sediments deposited within estuaries contain a wealth of information about the history and time course of changes of estuary biogeochemistry and ecology (Savage et al., 2010; Wigand et al., 2014). Here, we provide an approach to determine historic estuarine condition by combining estimates of historical nitrogen (N) loading to an estuary with sediment core biogeochemical properties and an ecological model of responses to sedimentation and nutrient loading. We hypothesize that these tools will provide an estuary-specific time series of ecological response to contaminant load changes and allow examination of estuarine condition tipping points (e.g., sediment and/or nutrient inputs at which seagrass or macrobenthic habitats and health were lost). Such knowledge could be used to determine realistic restoration objectives and their associated contaminant input levels, thereby providing management targets for restoration and catchment limit settings.

To test this hypothesis, we relate hindcasts of estuarine nutrient loading derived from catchment land use information to sediment accumulation rates (SAR) and geochemical and biological attributes in dated strata from sediment cores in a eutrophic estuary: New River Estuary (NRE) in Southland, New Zealand (Stevens et al., 2022). This estuary has a history of increasing nutrient and sediment loading corresponding with land clearance for sheep and beef farming from the early 20^th century, with intensive dairy production in the late 20^th century. It has also received urban-sourced wastewater nutrient loading from the adjacent Invercargill city, all of which have led to ecological change (Dudley et al., 2017; Stevens et al., 2022; Roberts et al., 2022a). We use contents of mud, photosynthetic pigments, carbon (C) and N isotopes, and bulk C and N from dated sediment core layers to describe estuary geochemical and biological conditions through time, back to 1847. We then estimate historic and present-day nutrient loading to the NRE by applying archival GIS-based land use layers within a nutrient runoff model (Elliott et al., 2016), and estuary nutrient concentrations arising from those loads using an estuary mixing model (Plew et al., 2018). Finally, these inputs are combined with other known drivers of estuary condition within a Bayesian network model (Zeldis and Plew, 2022) to hindcast values of a range of estuary ecological health indicators and an overall estuary health score, from 1847 to the present day.

To date, Bayesian modelling has only been applied to make present-day descriptions of ecological condition for New Zealand estuaries. However, by combining sediment core results with historic nutrient load and concentration estimates, the modelling can show ‘from where we have come’ in terms of estuary ecological condition over time compared to an earlier, pre-impacted state. Examining those historic states (including the pre-European colonization less impacted state), perspectives may be gained on when important changes in ecological condition occurred, along with co-occurring land use, and nutrient and sediment load forcings. By providing views of past ecological states and their drivers, we can inform efforts to set management targets and accompanying load limits for estuary remediation.

2 Methods

2.1 Sediment core collection

Replicate sediment cores (n = 3) were collected from two depositional sites (Upper: 46°25’51.7”S 168°19’35.5”E, Central: 46°26’28.4”S 168°19’39.8”E, WGS84) in the Waihopai Arm of the NRE, Southland, Aotearoa New Zealand (Figure 1), a region of the estuary where the most significant sedimentation has occurred with observed eutrophic changes (increased macroalgal cover and reduced dissolved oxygen). Sediment cores were collected using an Aquatic Research Large Bore Universal Sediment Corer fitted with 1.2 m x 102 mm OD (95 mm ID) polycarbonate core barrels with sharpened ends. Sediment cores of 1.0 - 1.2 m were collected by driving the cores into the sediment using a 3 kg percussion hammer attached to the top of a polyethylene and stainless-steel core head unit. The head unit is fitted with a one-way check valve that seals automatically, creating a vacuum when retrieving the core. The sediment around the core was dug out and removed using spades to a depth of 600 - 700 mm to facilitate easier retrieval. Polyethylene end caps and EVA inserts were fitted after retrieval and taped in place with duct tape, prior to shipping. Each core was allocated to be subsampled for layer dating, physical and biogeochemical analysis (particle size, phytopigment analysis, organic matter content, carbon content, nitrogen content, δ¹³C, δ¹⁵N, and sediment metal content), or stored (Table 1).

Figure 1

Figure 1 Location of the Upper and Central sediment coring sites within (A) Aotearoa/New Zealand showing the location of the Waihopai Estuary and its upstream catchment within Southland, (B) the catchment river network, with nutrient load monitoring sites and point sources, and (C) the lower New River Estuary and its main tributaries, the Waihopai River and the Ōreti River, and location of the Invercargill Wastewater Treatment Plant (WWTP) discharge at Clifton.

Table 1

Table 1 Details of the sediment cores collected from the two locations (Upper, Central) in the Waihopai Estuary, core length retrieved (cm), estimated compression during collection, and the analyses undertaken.

2.2 Geochronology

To determine the age of sediment preserved in the deposited layers, the cores were sampled in 1-cm thick increments (40 – 60 g dry sediment sample) every 5 cm downcore from one replicate core per site (UW1C, n =13, CW3C, n = 13, Table 1) down core and dated using gamma spectroscopy. The activity of excess ²¹⁰Pb and ¹³⁷Cs (number of disintegrations per second, Bq kg^-1) was counted at the ESR National Radiation Laboratory for 23 hours using a Canberra Model BE5030 hyper-pure germanium detector. Excess ²¹⁰Pb was determined from the ²²⁶Ra (t_½ 1622 y) activity and uncertainty determined following Swales et al. (2021). Excess ²¹⁰Pb was used to determine time-averaged SAR using a regression analysis of natural log-transformed data (Swales et al., 2021). The maximum depth of ¹³⁷Cs was used to verify the excess ²¹⁰Pb dating estimates. ¹³⁷Cs deposition was first detected in 1953 in New Zealand and therefore provides an approximate baseline for appearance in the sediment layer, with peak atmospheric deposition in 1963/64 (Matthews, 1989; Swales et al., 2021). Following initial dating, further refinement was provided by additional analyses of layers subsampled from a second replicate (UW1A n = 3, CW3B n = 2) to confirm the ²¹⁰Pb SAR inflection points in the regression analysis. However, the SAR estimates carry a certain amount of uncertainty. Fitting the unsupported ²¹⁰Pb data underpinning the SAR calculations involves some subjectivity when selecting the most likely break points in their trajectory. When estimated ages of adjacent core slices are widely spread, this can introduce error in the estimation of dates of the break points, which may have contributed to the differences found in the estimated break point dates between the Upper (1983) and Central (2006) cores. Also, naturally occurring frequent deposition and erosion of the surface sediment layers (such as by high runoff or storm events) mean that long-term sediment accumulation captures only a small fraction of the sediment deposited at any one time (Bentley et al., 2014). Additionally, bioturbation of sediments, resulting in mixing of the sediment surface with deeper layers, results in disruption of clear strata formation until sediment is buried below the maximum mixing depth (estimated to be approximately 15 cm at another similar New Zealand estuarine site (Bentley et al., 2014). Below this horizon, we assumed that C and N were no longer remineralised and lost across the sediment water interface, preserving the C and N contents and their ratios within the deeper sediment strata (Canfield, 1994). We also assumed that these deeper sediments were mixed during the period when they were within the bioturbation horizon, which may disrupt sediment core strata. Finally, core compression during collection adds uncertainty to age estimations of the sediment layers. We assumed the compression occurred predominantly in the upper layers of the core that had lower bulk density, meaning recent age estimations were more affected, while the age determination of the older layers was more accurate.

2.3 Physical and biogeochemical analyses

Sediment cores were sliced into 1 cm samples based on the calculated age model of the layers to obtain a time series of core sediment properties. For particle size analysis, a subsample of each slice (10 g) was taken, lyophilised and ultrasounded for 4 minutes to break up flocculants. The mean particle size (µm) and sediment mud content (%) were determined via volume particle diameter using an EyeTech time-of-transit (obscuration time) instrument (AmbiValue B.V., Dussen, Netherlands). A subsample (10 g) was taken for phytopigment analysis for chlorophyll-a and phaeophytin-a content (µg g^-1) following Sartory (1985). Pigments were extracted in 95% ethanol and determined spectrophotometrically. A preservation index (the ratio of chlorophyll-a:phaeophytin-a) was determined as a proxy for organic matter preservation (Le Guitton et al., 2015). Chlorophyll-a degrades more slowly in low oxygen conditions (Sun et al., 1991, 1993), therefore a greater preservation index can indicate either high deposition rates of chlorophyll-a or hypoxic/anoxic conditions of the overlying water (Rabalais et al., 2004). Organic matter content (%) was determined on a subsample (25 g) using loss on ignition. Contents of carbon (C, %), nitrogen (N, %) and stable isotopes δ¹³C and δ¹⁵N (‰) were determined on subsamples (5 g) using a DELTA V Plus continuous flow isotope ratio mass spectrometer linked to a Flash 2000 elemental analyser using a MAS 200 R autosampler (Thermo-Fisher Scientific, Bremen, Germany) at the NIWA Environmental and Ecological Stable Isotope Facility in Wellington, New Zealand. For sediment metal concentrations, subsamples (50 g) were air dried at 35^°C, acid digested (nitric/hydrochloric acid: US EPA 200.2) and assayed using ICP-MS mass spectrometer.

2.4 Modelling historic nutrient load to the New River Estuary

Point sources (e.g., discharges from waste-water treatment plants either directly into the estuary or into tributaries of the estuary) and diffuse sources (run-off from land) both contribute substantial nutrient loads to the NRE.

2.4.1 Point sources of nutrients

Loads of the four major point sources within the catchment (i.e., the Invercargill City Council Wastewater Treatment Plant (ICC WWTP), and three meat processing plants: Prime Range Meats, Alliance Lorneville, and Alliance Makarewa) were determined for all years with available data by summing daily loads derived from multiplication of total nitrogen (TN) concentration data and flow data. This method is reliant on accurate reporting of flow on all days when discharges occur; however, unreported flows have been observed on some dates between 2012 and 2013 (Robertson and Stevens, 2013a). Because it is not possible to accurately quantify unreported losses through time, we calculate loads using reported discharges and note that the loads from point sources may be higher than those stated here.

Wastewater loads from the ICC WWTP to NRE in years for which no data were available were calculated by assuming a linear relationship between wastewater TN and TP discharges and human population within the Invercargill urban area. This relationship was calculated based on population and ICC WWTP discharge records between 2004 and 2014. This period spanned the earliest date that daily consent data were available (2004) until the diversion of the Prime Range Meats meat processing plant flow to the ICC WWTP in 2015. Population data were sourced from Statistics New Zealand, and populations for non-census years were calculated by linear interpolations between census years.

Discharges from meat processing within the catchment in years for which no data were available were calculated from relationships between total meat processing plant TN and TP discharges and stock units (number of livestock) within the Southland Region. This relationship was calculated from discharge and stock data between 1999 and 2011. Annual Southland stock data were taken from Statistics New Zealand records.

2.4.2 Diffuse (non-point) sources of nutrients

Diffuse source nitrogen loads to the NRE between 1840 and 2020 were estimated using the Catchment Land Use for Environmental Sustainability (CLUES) model (Elliott et al., 2016). The CLUES model calculates annual riverine loads of nitrogen (as TN), based on land use within sub-catchments covered by the area of one or more of 19 land use classes representing key rural enterprises, native and exotic forest, and urban areas.

The land use together with animal stock numbers, physical catchment characteristics, rainfall and mean flows were used in a regression model to predict the total loads. Land use input data for the pre-European colonisation and post-Māori arrival period for the year 1840 were calculated from the vegetation layers of Leathwick et al. (2012), adapted from the work of Leathwick (2001) and Leathwick et al. (2003) which predicted potential vegetation cover composition from regressions relating the distributions of major canopy trees to environmental variables. This potential vegetation cover was adjusted to account for clearance of lowland indigenous forest prior to European arrival, using the 1840 grassland extent map of Ausseil et al. (2011). Land use for 1935 was derived by digitising the map of grasslands surveyed by Hilgendorf (1935). The vegetation cover from the first edition of New Zealand Land Resource Inventory was used as a 1978 Land use map (Landcare Research Ltd., 2018). Land use maps for the years between 2000 and 2020 were from Land Cover Database maps (LCDB), ranging from LCDB1 to LCDB5 (Landcare Research Ltd., 2013), intersected with Agribase land use boundaries (AsureQuality, 2008). Due to the lack of Agribase data for the year 2012, the 2012 CLUES model run was based on the 2008 land use/agribase layer but with updated coefficients using 2012 stock intensity data. Slope data for the catchment was acquired from the Land Environments of New Zealand classification (LENZ, Leathwick et al., 2002). Generated land use maps from 1840 to 2020 are shown in Figure 2. Agribase and LENZ were used to split the LCDB classification of grassland into pastoral land uses for different stock types that are characterised by different contaminant yields (e.g., dairy, deer, other animals and lowland intensive, hill county and high-country sheep and beef farming). Other land use classes such as cropland, orchards and perennial crops were further subdivided into specific classes, such as apples, kiwifruit and grapes, or for crops such as onions, potatoes, or maize. The rest of the LCDB classes were mapped directly without subdividing (e.g., native and exotic forests, scrub, tussock, urban areas, rivers and other water bodies).

Figure 2

Figure 2 Derived land use maps of the New River Estuary catchment at seven timepoints from the pre-European (1840) period to 2020.

The CLUES model uses two metrics to adjust nutrient loads generated on each land use class: a proportional adjustment for stock density, and an adjustment for land use intensity, both relative to the 2008 baseline used in model development (Elliott et al., 2016). Stock numbers (used to calculate stock density) were obtained from Agribase (AsureQuality, 2008) and summed accordingly for combined classes (e.g., sheep and beef). The maps were clipped to the New River catchment boundary and added to the CLUES GIS application directly and run as different scenarios to predict the loads. In the absence of stock density data prior to the year 2000, stock density for ‘sheep and beef’, ‘dairy’, and ‘deer’ land use classes were approximated using changes in stock numbers (or stock unit equivalent in the case of sheep and beef) in the Southland region (excluding Fiordland and the islands which have negligible stock), relative to the 2008 baseline. Adjustments for land use intensity for sheep and beef, and deer land classes were made using annual excretal N loads for these stocks taken from the New Zealand Greenhouse Gas Inventory 1990-2019 (Ministry for the Environment, 2021). Excretal N loads for the years prior to 1990 use the 1990 values. Adjustments for the dairy land class were made using milk solids per cow statistics based on relationships between milk solids output per cow and nutrient inputs per cow to dairy farmland (Luo and Ledgard, 2021). Area of land use by type, livestock density, and estimated nitrogen load are shown in Figure 3.

Figure 3

Figure 3 New River Estuary catchment (A) land use, (B) livestock density, and (C) estimated nitrogen load input.

2.4.3 Calibration of CLUES model results

Recent observed nutrient data from the six monitoring sites in rivers draining into the NRE were used to calibrate CLUES model results. These include daily flow and monthly concentration data from Makarewa River at Wallacetown (2007-2021), Ōreti River at Lumsden Bridge (2006-2020), Ōreti River at Three Kings (2007-2020), Ōreti River at Wallacetown (2006-2020), Otepuni Creek at Nith Street (2007-2020), and Waihopai River upstream of Queens Drive (2007-2020). Daily and annual TN concentrations at these locations were calculated by interpolating between monthly values for each site. Daily loads were then calculated by multiplying daily flow, nutrient concentration values and time steps. A linear regression was made between observed annual loads (t y^-1) and CLUES model results for years with land use maps and observed data. A separate regression relationship was made for sites on the Ōreti and Makarewa Rivers and two other sites at Waihopai River and Otepuni Creek. There was a reasonable agreement between observed and modelled TN with R² values of 0.88 and 0.97 for combined TN load from Ōreti and Makarewa Rivers sites and Waihopai and Otepuni sites, respectively. The uncertainty associated with the regression relationships between observed annual TN load and model estimates was assessed using the standard error of the regression. This measure provides insight into the variability or margin of error in the predictions made by the model. The standard error between observed annual TN load and CLUES model estimates were 304 t y^-1 for the Ōreti and Makarewa River reaches, and 31 t yr^-1 for the Waihopai and Otepuni. These high regression values and relatively low error estimates, compared with the average of observed TN values underscore the effectiveness of the model in accurately predicting TN levels in the study area. Coefficients from these regression relationships were then used to calibrate CLUES model results for years without observed data. To calculate annual NRE TN load, calibrated CLUES results for all terminal river reaches which drain into the NRE were summed.

A separate linear regression relationship was developed to estimate TN loads for years between 1840 and 2020 without available land use maps, using dairy stock excreta N load for the Southland region (for the years 2000-2020) and total stock excreta N (for the years 1840-1999) as the predictors. Stock excreta data from Ministry for the Environment, and stock numbers from Statistics New Zealand were used for this calculation (Ministry for the Environment, 2021).

To investigate the sensitivity of the CLUES model estimates to changes in stock density values, the data for the year 1935 was run in the CLUES model using 0.5 and 2 times the original stock density values, but with the same land use intensities. By reducing the stock density values to 0.5 of the original values, the calibrated TN load was reduced by 8.2% (123 t y^-1). When doubling the stock density values, the TN load increased by 9.6% (143 t y^-1). We note that the difference between modelled TN loads caused by these changes in stock density was less than the uncertainty in the CLUES model estimates.

2.5 Dilution modelling

Total nitrogen concentration in NRE was calculated as ‘potential TN’, which is defined as the concentration in the estuary the absence of uptake by algae or changes due to non-conservative processes such as denitrification (Plew et al., 2018). As such, potential TN represents the nitrogen available for primary production. Past and present estuary potential TN concentrations were estimated using estuary hydrological dilution models developed within ETI Tool 1 (Plew et al., 2018, 2020), which enables users to assess susceptibility of estuaries to eutrophication based on their nutrient loads and their flushing/dilution characteristics. Potential TN concentrations were calculated based on the relative proportions of riverine and marine water within the estuary, using salinity to tune a modified tidal prism model (Plew et al., 2018). For NRE, the tidal prism model was tuned using volume-averaged salinity extracted from a calibrated three-dimensional hydrodynamic model (Measures, 2016). The potential TN concentrations represent a time- (annual) and volume-averaged concentration for the entire estuary.

2.6 Ecological state modelling

Past and present estuary ecological states were estimated using Tool 3 of the online Estuary Trophic Index (ETI) toolset (Zeldis et al., 2017). ETI Tool 3 is a Bayesian Belief Network (BBN, Zeldis and Plew, 2022) that uses knowledge of the ecological connections between drivers of estuary trophic condition (e.g., estuary type, nutrient and sediment loads, flushing rate, etc.) and responses of estuary health indicators to calculate scores for the individual indicators, and a combined overall estuary health (ETI) score. Indicators are distinguished as ‘primary’, being macroalgae and phytoplankton that have demonstrably increasing responses to nutrient enrichment, and ‘secondary’, being ecological attributes that show symptoms of estuary health impairment (e.g., high sediment carbon, reduced sediment oxygen, impacted seagrass and macrobenthic communities, in response to nutrient loads and sedimentation. The ETI BBN combines the primary and secondary indicator scores into the final ETI score (Figure 4), where the score of the most impaired primary indicator is combined with the average of the secondary indicator scores and is discretised into 10 bands (Zeldis and Plew, 2022). ‘Intermediate calculation’ nodes do not directly input to the scores but are used to calculate some of the indicators (Figure 4).

Figure 4

Figure 4 Schematic of ETI Tool 3 BBN. Information on ‘drivers’ (blue nodes) are input by the user and the BBN calculates values of primary and secondary trophic indicators (yellow and pink nodes). Primary and secondary indicator values are used to produce ETI primary and secondary indicator scores, which are then combined to give the final ETI ‘performance’ score and band (red nodes). Intermediate calculation nodes (grey nodes) are used to estimate linked indicator values. Sediment accumulation rate, % mud, and % TOC were set using the sediment core results (Table 3). Because these are the downstream ‘children’ of the sediment load driver node (Figure 4), that node was not set. All other driver nodes were set (Table 4), including hindcasted potential TN concentrations obtained from ETI Tool 1 (Table 3). Adapted from Zeldis and Plew (2022) under CC-BY 4.0.

In NZ water quality policy statements (e.g., the National Policy Statement for Freshwater Management; NPS-FM; Ministry for the Environment, 2023), freshwater quality health variables are scored using a four-level banding (A-D). Therefore, the ETI score can be recast as four bands (A-D) if needed, to make comparison with the 4-state banding used within the NPS-FM. The estuary health conditions associated with the 4-state bandings of the ETI final score are described in Table 2.

Table 2

Table 2 Estuary health conditions associated with the bandings (A – D) and ETI scores 0 to 1 of the ETI Tool 3 BBN (Zeldis and Plew, 2022).

In default applications of the ETI Tool 3 BBN, all inputs are provided solely by values of the driver nodes (Figure 4). However, the BBN can also be informed by inputting values of indicator or intermediate calculation nodes, if these are known a priori, to improve its predictions (Uusitalo, 2007; Zeldis and Plew, 2022). Accordingly, in this application, intermediate calculation node values for sediment accumulation rate and % mud, and secondary node values for sediment % TOC (all derived from the sediment core results: Table 3), were input a priori to the BBN model runs. Because these are the downstream ‘children’ of the sediment load driver node (Figure 4), we did not set the value of the sediment load node in the BBN model runs. All other driver nodes were input, including hindcasted potential TN concentrations obtained from ETI Tool 1 dilution modelling (Table 4).

Table 3

Table 3 Values used to set potential TN (mg m^-3), % mud, sediment accumulation rate (SAR) (mm y^-1) and % TOC nodes in the BBN, for years which had values for % mud and/or % TOC derived from the Central and Upper cores.

Table 4

Table 4 Values used to set driver nodes in the BBN, for all years and both cores.

We provide resulting scores for four eutrophication and sedimentation associated indicators: Macroalgal ecological quality rating (EQR) is provided with thresholds of 0.8, 0.6 and 0.4 between A–B (minimal-moderate), B–C (moderate-high) and C–D (high–very high) bands of macroalgal eutrophication respectively (Robertson et al., 2016; Plew et al., 2020); Apparent Redox Potential Discontinuity (aRPD) depth, the depth of the boundary between oxic near-surface sediment and the underlying suboxic or anoxic sediment, where shallower scores are related to reduced benthic substrate oxygen penetration and habitat quality; Seagrass health driven by SAR and macroalgal growth (EQR) where larger scores indicate greater deleterious impacts; and a macrobenthos NZ AZTI Marine Benthic Index (AMBI, Borja et al., 2000) ranging from ‘Normal’ (band B: 1.2 – 3.3) to ‘Transitional to pollution’ (band C: 3.3 – 4.3), to ‘Polluted’ (band D: > 4.3).

The relationships between nodes of the ETI Tool 3 BBN are detailed in Zeldis and Plew (2022) and its supplementary information file. These relationships and modifications to the BBN since publication of Zeldis and Plew (2022) are provided in Supplementary Datasheet 1, Section 1 with modifications also in the ‘Change Log’ file available in the online ETI Tool 3.

3 Results

3.1 Sediment accumulation rate

Sediments at the Upper site were dated back to 1847 and at the Central site to 1874. At the Upper site, the SAR from 1847 (91 cm below surface) to 1983 (60 cm below surface) was 2.4 mm y^-1 and from 1983 to the present day it was 16.7 mm y^-1 (Figure 5; Supplementary Datasheet 1, Figure S1). At the Central site, SAR from 1874 (91 cm below surface) to 2006 (40 cm below surface) was 3.8 mm y^-1 and from 2006 to the present-day SAR was 33.6 mm y^-1 (Figure 5; Supplementary Datasheet 1, Figure S1).

Figure 5

Figure 5 Log-linear regression graphical interpretations of unsupported ²¹⁰Pb profiles with sediment core depth for (A) the Upper (orange) and (B) the Central (blue) sites, and (C) calculated sediment accumulation from the modelled log-liner regression fitted ages at the Upper (from 1847) and the Central (from 1874) locations in the New River Estuary.

3.2 Sediment physical and biogeochemical properties

Sediment organic matter and C content increased from the oldest dated layers to the present. Sediment OM ranged from 1.1% to 6.7% at the Upper site and from 1.4% to 8.2% at the Central site (Figure 6A) and % C ranged from 0.32% to 2.95% at the Upper site and from 0.16% to 3.12% at the Central site (Figure 6B).

Figure 6

Figure 6 Sediment core physical and biogeochemical properties at the Upper (orange) and Central (blue) sites in the New River Estuary: (A) organic matter content (%), (B) carbon (C) content (%), (C) δ¹³C (‰), (D) mud content (%), (E) nitrogen (N) content, (F) δ¹⁵N (‰), (G) C:N molar ratio, (H) chlorophyll-a content (µg g-1), (I) organic matter preservation index (chlorophyll-a:phaeophytin-a).

δ¹³C values, corrected for the Suess effect (the decrease in atmospheric δ¹³C caused by CO₂ emissions of fossil fuel depleted in ¹³C; Verburg, 2007), decreased from -24.47‰ to -26.64‰ in the Upper core between 1847-1983 (Figure 6C). In the Central core there was little variation (range from -25.27‰ to -26.16‰) and no clear trend between 1874 and 1985 (Figure 6C). However, sediment laid down in recent decades (i.e., post-1985) showed a steady increase in δ¹³C values at both locations up to -23.67‰ at the Upper site and -23.38‰ at the Central site. The conversion of the δ¹³C data to correct for the Suess effect was crucial: without it δ¹³C in the first two decades of this century was lower than in the early 20th century, suggesting a decrease in dissolved carbon demand and lower primary productivity this century compared with a century ago. This would not be consistent with the multi-pronged evidence for eutrophication from other sediment characteristics.

Mud content had steep increases in both the Upper and Central cores from the older sediments (1880 and 1916 at the Upper and Central sites respectively) to the present (Figure 6D, Table 4). While the increase in mud content at the Central site was somewhat variable prior to the 1980’s, its subsequent persistent 100% values (at both core sites) indicated formation of a very muddy environment in recent decades.

Sediment core N content increased from 0.03% to 0.34% at the Upper site and from 0.01% to 0.36% at the Central site (Figure 6E) from oldest dated layers sampled to the present day. From ~ 1850 to 1990, δ¹⁵N values were within the range of marine-derived N sources in southern New Zealand of ~ 5-6‰ (Cornelisen et al., 2007; Durante et al., 2021) and slightly below the 6.6 – 8.8‰ range of δ¹⁵N values given for macroalgal (Ulva spp) tissue collected from unimpacted, exposed coastal sites around New Zealand Barr et al., 2013; Figure 6F). Increased δ¹⁵N values, above the marine derived 5-6‰ range, are apparent from ~ 1990 and both cores show particularly sharp increases in δ¹⁵N values in sediments deposited since the early 2000’s. C:N ratios (Figure 6G) largely followed inverse trends to those seen in δ¹³C and δ¹⁵N values, with steady decreases from ~11 to 9 from the early 2000’s in both cores.

Sediment chlorophyll-a concentrations were relatively stable in both cores over the 20th century (Figure 6H), but from ~2000 onwards, sharp increases occurred at both core sites. The preservation index (chlorophyll-a:phaeophytin-a ratio) was relatively invariant until the early 2000’s, when it underwent sharp increases at both sites (Figure 6I).

3.3 Total N load and potential TN time series

CLUES-modelled TN loads to the estuary from freshwater and point sources increased from 947 t y^-1 in 1840, up to a maximum of 5715 t y^-1 in 2007 (Figure 7A). Hindcasted estuary potential TN concentration (Figure 7B) derived from dilution modelling increased ~4-fold from 1840 to the present, consistent with the increasing total load to the estuary. During the more recent period of rapid increase from 2000-2005 to 2015-2019 (5-year averages to reduce interannual variability), potential TN increased from 449 to 713 mg m^-3. There were very high correlations between the potential TN and N in the sediments at corresponding times in their respective time series (Figure 7C: Upper core R² = 0.949, TN = 110 + 1751 (% N); Central core R² = 0.926, TN = 112 + 1916 (% N)).

Figure 7

Figure 7 (A) CLUES modelled (white circles) and interpolated values for TN load (t y^-1), (B) potential N (mg m^-3) in the New River Estuary from 1847 to 2019 at the Upper (orange) and Central (blue sites), and the (C) relationship between sediment core nitrogen content (%) in the dated sediment core layers and the calculated potential N (mg m^-3). Upper core R² = 0.949 (y = 1751.40x - 110.20), Central core R² = 0.926 (y = 1916.36x - 112.00).

3.4 Ecological state modelling

NRE ecological state over the time series was described by BBN indicator and ETI scores (Figure 8) that resulted when the BBN was informed by the driver settings and observed sediment core values in Tables 3, 4. Hindcasted macroalgal EQR was stable at band A (healthy: greater than or near 0.8) for much of the time series prior to the 1950’s (Figure 8A). It worsened through band B (0.6 to < 0.8) to band C (0.4 to < 0.6) by the late 1980’s, coincident with gradually increasing potential TN (Figure 7A), but then declined more steeply starting in the late 1990’s in response to sharply increasing potential TN, reaching band D (< 0.4) by ~2010.

Figure 8

Figure 8 Responses of trophic indicator scores (A) macroalgae ecological quality rating (EQR), (B) sediment aRPD, (C) seagrass impact, (D) macrobenthic NZ AMBI, and (E) ETI score, with associated bandings (for (A–D) c.f. Table 2).

BBN-predicted aRPD (Figure 8B) was relatively stable prior to the 1930’s at band A (> 4 cm depth), after which it started shallowing initially in response to increasing % TOC (Figure 6B; Table 3) and then declining macroalgal EQR (Figure 8A) after the 1950’s. It reached band D (< 1 cm depth) in the early 2000’s, driven both by increasing % TOC and steeply declining macroalgae EQR. Overall, the shallowing aRPD indicated that the depth of the boundary between oxic near-surface sediment and the underlying suboxic/anoxic sediment was likely to have shallowed to very near the surface of the estuary bed.

Seagrass decline (Figure 8C) was predicted to be minor to moderate (bands A-B: scores < 0.5) until ~1983 for the Upper core area and ~ 2006 for the Central core area, after which they underwent severe decline to band D (> 0.8) at both sites. These dates corresponded with inflection points in SAR at the respective sites (Figure 5C), although the transition was less well defined for the Upper core because of the lack of SAR results between 1966 and 1983 (Table 3). For both cores, D-band status corresponded with steep declines in macroalgal EQR (Figure 8A).

Macrobenthic NZ AMBI index (Figure 8D) indicated ‘normal’ (≤ 3.3) status until 1966 and 1953 for the Upper and Central cores respectively, after which they assumed ‘transitional to pollution’ (≤ 4.3) status, eventually reaching ‘polluted’ status (> 4.3) in 1998 at both sites. Their worsening scores occurred in response to increasing % TOC and % mud (Figures 6B, D) and finally, in response to steeply worsening macroalgae EQR (Figure 8A).

Overall NRE estuary health, indexed by the ETI score (Figure 8E), was in band A (≤ 0.25: Table 2) until the early 20^th century. It degraded slowly through band B (> 0.25 to ≤ 0.50) from the 1930’s to the 1980’s, but then worsened rapidly through band C (> 0.50 to ≤ 0.75), reaching band D (> 0.75) by 2007. The degradation arose as an interaction of steadily increasing SAR, % mud, % TOC and potential TN, but its acceleration in the late 1990’s was associated with rapidly worsening macroalgae EQR state, driven by the sudden increase of potential TN.

4 Discussion

In this study, we reconstructed the historic trophic condition of the NRE. We used historical land use information to calculate nutrient loading to the estuary and determine estuary-wide time series of water column TN concentrations. We compared TN concentration time series to biogeochemical indices in dated sediment core layers at two locations in the estuary. Both TN concentration model results and sediment core chemistry indicated catchment and estuary N loads increased greatly from the year 2000 to the present day, and further modelling using a BBN (Zeldis and Plew, 2022) provided a time series of ecological state values. Our results indicate that the currently eutrophic NRE (Robertson et al., 2017; Forrest et al., 2022; Roberts et al., 2022a), whilst historically muddy, was not necessarily eutrophic in the past with a much lower SAR. While overall, the NRE showed increasing degradation, it showed differing responses to differing land use profiles through time (i.e., sheep and beef, then dairy), firstly with increased muddiness, and then macroalgal eutrophication.

4.1 Temporal trajectories of sediment and nutrient loading

The sediment core results from the upper NRE showed clear sign of increased sedimentation with increased muddiness and % TOC commencing ca. early-20^th century. SAR increased steadily at relatively low rates until inflection points appeared in the late 20^th and early 21^st centuries, with sharply increased rates thereafter. These rates concur with sedimentation estimates from sediment cores collected in 2007 by Robertson and Stevens (2007) and in 2017 by Brown (2019) in depositional zones in the NRE. Sedimentation rate in the Upper Waihopai Arm was estimated to have increased from 3 mm y^-1 from 1906 to 1967 to 28 mm y^-1 from 2001 to 2007 by Robertson and Stevens (2007), from 13.3 mm y^-1 from 1915 to 1935 to 22.4 mm y^-1 from 2009 to 2017 by Brown (2019), and from 5.9 mm y^-1 from 1956 to 2017 to 17.5 mm y^-1 from 2007 to 2017 in the Lower Waihopai Arm by Brown (2019). Sediment plate monitoring conducted in in Lower Waihopai Arm estimated sedimentation rates of 30-40 mm y^-1 (including both fine sediment an macroalgae entrained within the sediments) from 2007 to 2020 (Roberts and Ward, 2021) and 18.1-23.2 mm y^-1 from 2007 to 2022 (Forrest et al., 2022).

Increases in several core biogeochemical indicators (% OM, % C and % N) since the 1800’s, and at a more rapid rate from mid-20^th century, indicated gradually accelerating sedimentation and eutrophication in the upper NRE as catchment stock numbers and human population increased (AsureQuality, 2008; Grant, 2023). Additionally, and in contrast to % OM, % C, and % N, between about 1980 and the early 21^st century, core sediments showed a sudden switch from little to no change in sediment δ¹³C, δ¹⁵N, chlorophyll-a, and C:N ratios since the 1800’s to steeply increasing sediment δ¹³C and δ¹⁵N values, large increases in chlorophyll-a content, and decreasing C:N ratios. C:N ratio decreases in organic material in mixohaline systems are indicative of increased contributions of marine-derived algal material relative to contributions of vascular plant material, and/or an increase in N availability to primary producers (Cloern et al., 2002; Barr and Rees, 2003; Savage et al., 2010). These changes taken together indicate increases in autochthonous algal primary production rates and deposition, as N loading accelerated during recent agricultural intensification in the catchment (Cifuentes et al., 1988; Savage et al., 2010). Sediment indicators of recent eutrophication mirrored the evidence seen in regular ecological monitoring of the intertidal zone of NRE since 2001, showing increases in algal biomass and changes in faunal communities (Robertson and Stevens, 2013b; Stevens et al., 2022).

In addition to increases in N availability from the estuary catchment, increases in sediment δ¹⁵N in the mid-late 20^th century likely indicate contributions of N from treated sewage. Hydrodynamic modelling of the NRE suggests that N from the Invercargill WWTP main municipal wastewater discharge (in its modern location at Clifton) likely forms a sizable proportion of total water column N at the locations of both cores (Measures, 2016). Given this, it is possible that some of the increases in sediment δ¹⁵N resulted from the commencement of effluent discharge from the WWTP (in 1969), construction of a secondary effluent treatment processing facilities (in 1992), and addition of maturation ponds and wetlands (in 2004) that likely increased N removal from effluent. Typically, the nitrification-denitrification processes that cause gaseous losses of N during secondary and tertiary effluent processing result in δ¹⁵N-enrichment of the remaining effluent-N (Savage, 2005; Alldred et al., 2023). Also, a range of N-cycling processes within catchments exhibit concentration-dependent expression of discrimination against δ¹⁵N, leading to progressively greater δ¹⁵N-enrichment of inorganic N pools passing to estuaries and coastal waters, as catchment land use (and N application) intensifies (Bedard-Haughn et al., 2003; Needoba et al., 2004; Dudley et al., 2010; Kritee et al., 2012; Wells and Eyre, 2019). As a result, sediments and primary producers in eutrophic, N-enriched estuaries and coastal environments influenced by rivers tend to have higher δ¹⁵N values (McClelland et al., 1997; Voss et al., 2000).

4.2 Impact of land use on nutrient and sediment loads

Increased nutrient and sediment loading to the NRE since the early-to-mid-20^th century resulted largely from changing rural land use and land use practices. Human populations within the Invercargill urban area have increased only marginally between the early 1970’s (50,681 in 1971) and recent times (54,204 at the 2018 census). We estimate that N releases via urban wastewater discharge to the estuary also increased little over this time. In contrast, our results indicate that inputs of N from agricultural land within the NRE catchment approximately doubled over the same period. This increase in N loss from the catchment mostly resulted from agricultural intensification, particularly via increases in stock numbers per unit area of farmland, and conversion of land used for less intensive agriculture such as sheep and beef farming, to more intensive use, particularly dairy farming. Note that while the area occupied by dairy farming in the 2010’s was similar to or less than that in 1935 and 1978, there were far more dairy cows after 2010. These patterns of agricultural intensification seen in Southland since the mid-20^th century mirror those seen in other areas of New Zealand (Wynyard, 2016; Ministry for the Environment, 2023).

Land use intensification can also lead to an increase in sediment generation from the catchment, a trend that has been observed in various catchments across New Zealand (Haddadchi and Hicks, 2020). Changes in hydrology as a result of land reclamation in the Waihopai Arm of the NRE in the early 1900s have been suggested to be a cause of increased sedimentation in this region, however we find that sedimentation following the reclamations remains low compared to recent SAR and that changes in land use are the likely cause of increased sedimentation observed since the 1980s. SAR increases in more recent times likely indicate a positive feedback often seen between fine sediment and nutrient loading increases, with muddy sediments and high nutrient loading favouring nutrient retention and macroalgal outgrowth (Robertson and Savage, 2018) and with macroalgae, in turn, trapping more mud (Zeldis et al., 2019; Stevens et al., 2022).

4.3 Using sediment cores to determine historic N loads and potential N concentration

We found very high correlations between the time series of estuary-wide potential TN in the water column (determined from N load and dilution modelling) and sediment N content in the dated sediment core strata. This suggested that historic potential TN water column concentrations are predictable from sediment core data. However, the potential TN concentrations represent an average value over the whole estuary, while the cores were taken from the impacted depositional zone in the Lower Waihopai Arm of the estuary. Spatially-distributed nutrient sampling (Dudley et al., 2022) in NRE and hydrodynamic modelling using tracers as a proxy for nutrients (Measures, 2016) show that the inner parts of the estuary, such as the Waihopai Arm, have higher nutrient concentrations than ‘estuary-averaged’ values. Additionally, the best/high-temporal resolution sedimentary records come from high-SAR environments. Sediment mixing by infauna is less effective as SAR increases so that physical processes dominate. This suggests that a historic time series of ‘local’ potential TN predicted for the inner estuary and other locations with different SAR rates would have a different relation with sediment N content, likely with a higher slope, than we show here for ‘estuary averaged’ values. While recognising this caveat, the fact that estuary-wide historic potential TN appears predictable from profiles of sediment core N content in NRE indicates that similar coring work on other suitable (depositional) sites in NRE and other estuaries would be beneficial, to determine if the strong relationship found here between modelled TN loads and sediment core N content is a general phenomenon.

4.4 BBN hindcasts of NRE ecological health

The BBN resolved a complex interaction of effects, operating at different times and time scales, within the overall pattern of degradation. The ETI BBN showed that NRE ecological health has deteriorated across all the ecological indicators calculated, particularly from the mid-20^th century to the present day. Macroalgal EQR deteriorated mostly from the late 1980’s (through the C band to D band) at both core sites, associated with the strong increase in potential TN occurring at that time (exceeding 400 mg m^-3 TN). The EQR results from both sites showed identical trajectories of change, due to the estuary-wide value of potential TN driving the EQR node. Monitoring of the area of intertidal flats of the estuary showed that estuary area covered by > 50% macroalgae increased from 43 ha in 2001 to 364 ha in 2016 (Robertson et al., 2017), to 479.7 ha in 2022 (Roberts et al., 2022a), and a decrease in calculated EQR score for the Waihopai Arm from 0.418 to 0.308 from 2012 to 2022 (Roberts et al., 2022a). This was consistent with, although slightly lower than, our calculated decrease in macroalgal EQR from 0.30 to 0.21 over the same period, following the rapid decrease in EQR after 2008. Monitoring showed that the EQR decline was driven primarily by expansion in substrate coverage of the red alga species Agarophyton spp. (formerly Gracilaria spp.), and to a lesser extent by Ulva spp., which were limited in extent in the NRE before 2000 (Stevens et al., 2022). These macroalgae are typically associated with increases in sediment mud content, loss of indigenous seagrass, decreased sediment oxygenation, and impacted macrobenthos (Robertson and Stevens, 2007; Robertson and Stevens, 2013b; Robertson et al., 2017). Since 2007, Agarophyton spp. has become dominant in the depositional areas of the upper NRE (Forrest et al., 2022), underpinning the declines in ecological conditions seen in the sediment core and BBN indicators retrieved from that location.

The modelled BBN aRPD showed a steady decrease from a ‘low impact’ value of 4.4 cm (A band) to 0.5 cm (D band), below the ‘exhaustion threshold’ levels of Sutula et al. (2014). Measurements of redox potential in the western Waihopai Arm from 2013 show an aRPD of < 1 cm (Robertson et al., 2017), consistent with the present day values yielded by the BBN. Further measurements from 2016 showed the sediment was anoxic at 3 cm depth (redox potential < 150mV, the depth at which the primary indicator measurement for assessment is taken) in the Waihopai arm, with worsening of the general condition of the sediments observed, including surface anoxia rich in sulphide (Robertson et al., 2017). Thus, both the BBN predictions and the present observations of aRPD implied greatly reduced benthic habitat quality and volume for macrobenthic communities (Green et al., 2014).

The macrobenthic taxon abundance and richness index NZ AMBI steadily increased (worsened) over the time series after ~1960, reflecting its dependence on bulk properties (% TOC and % mud: Robertson et al., 2016). Its sudden increase after ~1998 resulted from steeply worsening macroalgal conditions, interacting with % TOC and % mud (Green et al., 2014). Macrobenthic state in NRE was not routinely measured prior to 2001, but after 2001 it was shown to have ‘polluted’ status at many sites (NZ AMBI > 4.3) (Robertson et al., 2016, 2017; Forrest et al., 2022) at sites with high mud content and macroalgal cover. The worsening in NZ AMBI was consistent with findings that New Zealand estuaries with less mud tend to be more favourable for healthy macrobenthic communities (Clark et al., 2020).

Sporadic aerial imagery shows intertidal seagrass beds present on the central flats of the NRE (south of the Waihopai) in the 1940’s (first images in 1946) through to the 60’s however these beds were not recorded on imagery in the 80’s (there is no whole estuary imagery for the 70’s and imagery from 1996 was taken at high tide so inconclusive, K. Roberts, unpubl. data March 2023). In 2001 seagrass was mapped as part of the National Estuary Monitoring Protocol development and was localised to the Waihopai Arm, sheltered areas on the eastern margin and at the entrance (Roberts et al., 2022b). Seagrass was then subsequently lost over the period 2001-2022, the largest losses occurred between 2001 and 2012, with a 85.2% decrease (-83.1 ha) in seagrass areas with >50% cover in the NRE between 2001 and 2022 with much of this decline occurring in the Waihopai Arm (Roberts et al., 2022a, b). It is likely that these historic seagrass beds were affected by reductions in water clarity associated with increased SAR, when SAR transitioned from 2.4 and 3.8 mm y^-1 to 17 and 34 mm y^-1 at the Upper and Central sites, respectively. The latter values greatly exceeded the 5 mm y^-1 SAR threshold described by Zabarte-Maeztu et al. (2021) and Benham et al. (2019) for Z. muelleri decline. These declines were also coincident with heavy macroalgal growth, particularly Agarophyton spp., and eutrophication (Stevens and Forrest, 2020; Stevens et al., 2022) driven by the steeply increasing potential TN.

The BBN ETI scores slightly overestimate the worsening condition of the NRE after 2012 compared to ETI scores calculated by Roberts et al. (2022a) in 2022 obtained via direct estuary monitoring in the Waihopai Arm (Table 5), however, we can still see the trend in estuarine trophic quality degradation in both the direct estuarine monitoring-based estimates and the BBN estimates.

Table 5

Table 5 Estuarine Trophic Index (ETI) scores calculated at the two sediment coring sites (Upper and Central) and Waihopai Arm scores calculated by Roberts et al. (2022a) for 2012, 2016, and 2022.

Another potential reason for the discrepancy between the monitoring-based ETI estimates and the BBN estimates is that the BBN ETI scores may approach, but not reach, the theoretical limits of the scores (0 and 1). As described by Zeldis and Plew (2022) the distribution of ETI scores predicted by the BBN become compressed at the upper and lower limits, with achieved maximum and minimum values of ~0.9 and ~0.2. This is an artifact of the probability distributions built into BBNs generally, which can spread predictions for indicators, and their resulting scores, over several bands (Uusitalo, 2007). Mean values for each indicator and the ETI score were calculated from those probability distributions (using the centre of each band weighted by the probability distribution across each band), tending to bring the ETI scores away from the limits of 0 or 1. Confidence intervals were calculated in a similar manner. The consequence of this is that confidence intervals were wide (Figure 8) due, in part, to the discretised structure of BBNs (Zeldis and Plew, 2022).

The EQR values calculated within the BBN can be assumed to represent the whole of the estuary (Zeldis and Plew, 2022), because of the use of estuary-wide potential TN estimates and spatially representative EQR (Water Framework Directive - United Kingdom Technical Advisory Group, 2014). However, scores for other BBN indicators were sensitive to the inputs acquired from the localised area where the cores were collected (Waihopai Arm) including % TOC and % mud. Waihopai Arm is recognised to be in relatively poor ecological condition in the present-day estuary (Stevens et al., 2022), so this needs to be considered when interpreting results. The indicators predicted from those inputs (i.e., aRPD, seagrass, macrobenthos) were likely ‘more impacted’ than would have been estimated if ‘estuary average’ core-derived values had been available. It would therefore be beneficial to expand this study using cores from other suitable areas in the NRE recognised to represent a wider range of potential impact to provide a better estuary-wide comparison.

4.5 Historic ecological health and prospects for remediation

The coring, hindcasted potential TN and BBN results show ‘from where we have come’ in terms of NRE ecological conditions dating from mid-19^th to early 21^st centuries. Severe ecological impact has arisen from both sediment and TN loading, with strong impacts from sediments appearing earlier (mid-20^th century) than those from TN (late 20^th/early 21^st century). These differing impacts are associated with differing land use profiles over time (land clearance, grazing from sheep and beef farming eliciting the earlier sedimentation and moderate TN increases, accelerated N runoff from intensified dairy farming accounting for recent rapid TN increases). The correlations of sediment properties with runoff and associations with environmental degradation resembled responses in Jamaica Bay Estuary (Wigand et al., 2014), although in that case the runoff variability was driven largely by human population size changes, rather than agricultural intensification.

Because the ecological responses in NRE show temporal associations with agricultural land use profiles, we believe that the results may give guidance in designating restoration aspirations and targets for limit-setting, by identifying times of pivotal changes in estuary state and associated catchment sediment and nutrient inputs. For example, identifying the sediment loading conditions that triggered turbidity-driven subtidal seagrass and macrobenthic declines prior to ~1990 could indicate the sediment runoff conditions required to re-establish those communities, or those required to prevent deterioration of similar, but currently less impacted, estuaries. The large decrease in seagrass cover (in the early 21^st century) and continued decline of macrobenthic health could indicate the potential TN levels required to prevent macroalgal outbreaks and their overgrowth of intertidal seagrass beds and impacts on macrobenthos in this and other estuaries.

Although our study indicates previous estuary ecological states and their prevailing drivers, the compounding effects of sedimentation and nutrient pressures also indicate that reversion to previous states may not follow a linear recovery trajectory, and that interactions between the pressures could be difficult to predict. There could be hysteresis in the recovery trajectory (Duarte et al., 2009, 2015; Hale et al., 2023), potentially following the ‘hysteresis’ or ‘hysteresis with threshold’ patterns of Kemp et al. (2009), wherein degradation and recovery do not follow similar but opposite trajectories (Duarte et al., 2015). The compounding stressors would likely have different time scales of remediation. If catchment sediment loading and SAR were significantly reduced, there could still be a long time-lag before the deep deposits of mud (with its high C and N contents) were washed naturally from NRE. Also, macroalgae trap mud in NRE (Stevens et al., 2022), so unless their excessive biomass was reduced, sediment winnowing could be slow. Reducing TN loading to reduce macroalgal eutrophication would likely be beneficial, but the benefit could be slowed by persistence of sediment-bound N in the muddy deposits (Engelsen et al., 2008; Borja et al., 2010; Zeldis et al., 2019) and a response lag in Agarophyton spp. reduction due to utilization of internal nutrient storage pools (Dudley et al., 2022). High porewater DIN concentrations (~160 µmol) occur in NRE sediments (Robertson and Savage, 2018), and these may be important (along with water column DIN) in supporting its eutrophic biomasses of macroalgae. The situation contrasts with the significant, rapid recovery of ecological condition in Avon-Heathcote Estuary (Christchurch, NZ) upon reduction of its N load (Barr et al., 2020; Zeldis et al., 2020), which was facilitated by the naturally coarse (sandy), well-flushed sediments of that estuary.

This discussion is not intended to discourage remediation efforts, but to indicate possible outcomes. There would undoubtedly be benefits accruing from better runoff and wastewater management in NRE, including efforts to simultaneously reduce sediment and nutrient loads, in terms of realising improved ecological states. These results are also informative by illustrating historic sediment and nutrient runoff conditions that could lead to NRE-like conditions in other similar, but currently less impacted, estuaries.

5 Summary and conclusion

In this work we have provided reconstructed NRE ecological health records based on historical biogeochemical properties derived from a novel combination of aged sediment core data, models of hindcasted nutrient loads and estuary potential N concentrations, and BBN predictions of current and historical ecological conditions. Our results suggest that catchment N load (including diffuse and point source loads) is predictable from sediment core N and, more broadly, that the combination of methods used is effective in hindcasting estuary health over centuries. However, work is needed to further validate the ecological state hindcasts, and to establish how transferable the results are, both to other parts of NRE and to other estuaries.

Using sediment cores to derive historic estuarine condition, in the absence of historic time series measurements or catchment land use information, could inform nutrient load limits for other estuaries, to protect them from impacts of excessive nutrient and sediment loads. Runoff conditions corresponding to when significant shifts in estuarine health occurred, and the corresponding states of land use, could be used to propose targets for nutrient and sediment input reduction and associated limit-setting in the catchment, or determination of upper limits to prevent further estuarine decline. The combination of sediment core sampling with modelling tools demonstrated here shows ‘from where we have come’ in terms of history of change.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Datasheet 2. Further inquiries can be directed to the corresponding author.

Author contributions

RH: Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. JZ: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. BD: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. AH: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. DP: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. US: Data curation, Formal analysis, Investigation, Methodology, Resources, Validation, Visualization, Writing – review & editing. AS: Formal analysis, Funding acquisition, Investigation, Validation, Visualization, Writing – review & editing. KR: Investigation, Resources, Writing – review & editing. SO: Investigation, Resources, Writing – review & editing. PV: Formal analysis, Investigation, Methodology, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This study was funded under the Coasts and Estuaries platform by the National Institute for Water & Atmospheric Research Ltd.’s Strategic Science Investment Fund (FWCE2106, FWCE2206, CECE2305, CECE2402) from New Zealand’s Ministry for Business, Innovation and Employment core funding.

Acknowledgments

Thanks to Environment Southland for in-kind support with fieldwork and provision of historic maps for digitization. Thanks to Invercargill City Council and Adrian Cocker for provision of water quality monitoring data. Stable Isotope Analysis was carried out by the NIWA Environmental and Ecological Stable Isotope Facility in Wellington, New Zealand: Sample acidification and sample weighing was carried out by Josette Delgado and Anna Kilimnik. Standards weighing was carried out by Anna Kilimnik and Rahul Peethambaran. Stable isotope analysis was completed by Josette Delgado and Rahul Peethambaran. Daily stable isotope data were compiled by Josette Delgado. Data processing and calculations were carried out by Julie Brown. Sample processing, data analyses and calculations were overseen, checked, and amended where necessary by Sarah Bury. We thank Amandine Sabadel for assistance with interpretation of stable isotopes. Radio-dating was carried out by ESR labs, and we thank Michael Lechermann and Levi Bourke for assistance. We thank Ron Ovendon for their assistance with sediment core slicing and NIWA Hamilton Water Quality Laboratory for sediment sample analysis.

Conflict of interest

Authors KR and SO were employed by the company Salt Ecology.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2024.1374869/full#supplementary-material

References

Alldred F. C., Gröcke D. R., Leung C. Y., Wright L. P., Banfield N. (2023). Diffuse and concentrated nitrogen sewage pollution in island environments with differing treatment systems. Sci. Rep. 13, 4838. doi: 10.1038/s41598-023-32105-6

PubMed Abstract | CrossRef Full Text | Google Scholar

AsureQuality (2008) Agribase land use boundaries. Available online at: https://www.asurequality.com/services/agribase/ (Accessed 2022).

Google Scholar

Ausseil A.-G. E., Dymond J. R., Weeks E. S. (2011). “Provision of Natural Habitat for Biodiversity: Quantifying Recent Trends in New Zealand,” in Biodiversity Loss in a Changing Planet. eds. Grillo O., Venora G. (Rijeka: IntechOpen), Ch. 9. doi: 10.5772/24969

CrossRef Full Text | Google Scholar

Barbier E. B., Hacker S. D., Kennedy C., Koch E. W., Stier A. C., Silliman B. R. (2011). The value of estuarine and coastal ecosystem services. Ecol. Monogr. 81, 169–193. doi: 10.1890/10-1510.1

CrossRef Full Text | Google Scholar

Barr N. G., Dudley B. D., Rogers K. M., Cornelisen C. D. (2013). Broad-scale patterns of tissue-δ15N and tissue-N indices in frondose Ulva spp.; developing a national baseline indicator of nitrogen-loading for coastal New Zealand. Mar. pollut. Bull. 67, 203–216. doi: 10.1016/j.marpolbul.2012.11.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Barr N. G., Rees T. (2003). Nitrogen status and metabolism in the green seaweed Enteromorpha intestinalis: an examination of three natural populations. Mar. Ecol. Prog. Ser. 249, 133–144. doi: 10.3354/meps249133

CrossRef Full Text | Google Scholar

Barr N., Zeldis J., Scheuer K., Schiel D. (2020). Macroalgal bioindicators of recovery from eutrophication in a tidal lagoon following wastewater diversion and earthquake disturbance. Estuaries Coasts 43, 240–255. doi: 10.1007/s12237-019-00654-7

CrossRef Full Text | Google Scholar

Bedard-Haughn A., Van Groenigen J., Van Kessel C. (2003). Tracing 15 N through landscapes: potential uses and precautions. J. Hydrology 272, 175–190. doi: 10.1016/S0022-1694(02)00263-9

CrossRef Full Text | Google Scholar

Benham C. F., Beavis S. G., Jackson E. L. (2019). Tolerance of tropical seagrasses Zostera muelleri and Halophila ovalis to burial: Toward an understanding of threshold effects. Estuarine Coast. Shelf Sci. 218, 131–138. doi: 10.1016/j.ecss.2018.11.005

CrossRef Full Text | Google Scholar

Bentley S. J., Swales A., Pyenson B., Dawe J. (2014). Sedimentation, bioturbation, and sedimentary fabric evolution on a modern mesotidal mudflat: A multi-tracer study of processes, rates, and scales. Estuarine Coast. Shelf Sci. 141, 58–68. doi: 10.1016/j.ecss.2014.02.004

CrossRef Full Text | Google Scholar

Borja Á., Dauer D. M., Elliott M., Simenstad C. A. (2010). Medium-and long-term recovery of estuarine and coastal ecosystems: patterns, rates and restoration effectiveness. Estuaries Coasts 33, 1249–1260. doi: 10.1007/s12237-010-9347-5

CrossRef Full Text | Google Scholar

Borja Á., Franco J., Pérez V. (2000). A marine biotic index to establish the ecological quality of soft-bottom benthos within European estuarine and coastal environments. Mar. pollut. Bull. 40, 1100–1114. doi: 10.1016/S0025-326X(00)00061-8

CrossRef Full Text | Google Scholar

Brown D. (2019). Geochemistry and isotopic composition of sediment cores to understand the lithological and anthropogenic controls on eutrophication in the New River Estuary, Southland, New Zealand (Ontario, Canada: Masters, Laurentian University Sudbury).

Google Scholar

Canfield D. E. (1994). Factors influencing organic carbon preservation in marine sediments. Chem. Geology 114, 315–329. doi: 10.1016/0009-2541(94)90061-2

CrossRef Full Text | Google Scholar

Cifuentes L., Sharp J., Fogel M. L. (1988). Stable carbon and nitrogen isotope biogeochemistry in the Delaware estuary. Limnology oceanography 33, 1102–1115. doi: 10.4319/lo.1988.33.5.1102

CrossRef Full Text | Google Scholar

Clark D. E., Hewitt J. E., Pilditch C. A., Ellis J. I. (2020). The development of a national approach to monitoring estuarine health based on multivariate analysis. Mar. pollut. Bull. 150, 110602. doi: 10.1016/j.marpolbul.2019.110602

PubMed Abstract | CrossRef Full Text | Google Scholar

Cloern J. E., Canuel E. A., Harris D. (2002). Stable carbon and nitrogen isotope composition of aquatic and terrestrial plants of the San Francisco Bay estuarine system. Limnology Oceanography 47, 713–729. doi: 10.4319/lo.2002.47.3.0713

CrossRef Full Text | Google Scholar

Conley D. J. (1999). Biogeochemical nutrient cycles and nutrient management strategies. In: Garnier J., Mouchel J. M. (eds) Man and River Systems. Developments in Hydrobiology. (Dordrecht: Springer) vol 146. doi: 10.1007/978-94-017-2163-9_10

CrossRef Full Text | Google Scholar

Cornelisen C., Wing S., Clark K., Bowman M., Frew R., Hurd C. (2007). Patterns in the d13C and d15N signature of Ulva pertusa: Interaction between physical gradients and nutrient source pools. Limnology Oceanography 52, 820–832. doi: 10.4319/lo.2007.52.2.0820

CrossRef Full Text | Google Scholar

Duarte C. M., Borja A., Carstensen J., Elliott M., Krause-Jensen D., Marbà N. (2015). Paradigms in the recovery of estuarine and coastal ecosystems. Estuaries Coasts 38, 1202–1212. doi: 10.1007/s12237-013-9750-9

CrossRef Full Text | Google Scholar

Duarte C. M., Conley D. J., Carstensen J., Sánchez-Camacho M. (2009). Return to Neverland: shifting baselines affect eutrophication restoration targets. Estuaries Coasts 32, 29–36. doi: 10.1007/s12237-008-9111-2

CrossRef Full Text | Google Scholar

Dudley B. D., Barr N. G., Plew D. R., Scheuer K. (2022). Experiments to parametrise a growth and nutrient storage model for Agarophyton spp. Estuarine Coast. Shelf Sci. 264, 107660. doi: 10.1016/j.ecss.2021.107660

CrossRef Full Text | Google Scholar

Dudley B. D., Barr N. G., Shima J. S. (2010). Influence of light intensity and nutrient source on δ13C and δ15N signatures in Ulva pertusa. Aquat. Biol. 9, 85–93. doi: 10.3354/ab00241

CrossRef Full Text | Google Scholar

Dudley B. D., Burge O. R., Plew D., Zeldis J. (2020). Effects of agricultural and urban land cover on New Zealand’s estuarine water quality. New Z. J. Mar. Freshw. Res. 54, 372–392. doi: 10.1080/00288330.2020.1729819

CrossRef Full Text | Google Scholar

Dudley B., Zeldis J., Burge O. (2017). New Zealand Coastal Water Quality Assessment. NIWA Report No 2016093CH. Prepared for Ministry for the Environment. (Christchurch, NZ: NIWA). 84.

Google Scholar

Durante L. M., Smith R. O., Kolodzey S., Mcmullin R. M., Salmond N. H., Schlieman C. D., et al. (2021). Oceanographic transport along frontal zones forms carbon, nitrogen, and oxygen isoscapes on the east coast of New Zealand: Implications for ecological studies. Continental Shelf Res. 216, 104368. doi: 10.1016/j.csr.2021.104368

CrossRef Full Text | Google Scholar

Elliott A. H., Semadeni-Davies A. F., Shankar U., Zeldis J. R., Wheeler D. M., Plew D. R., et al. (2016). A national-scale GIS-based system for modelling impacts of land use on water quality. Environ. Model. Software 86, 131–144. doi: 10.1016/j.envsoft.2016.09.011

CrossRef Full Text | Google Scholar

Engelsen A., Hulth S., Pihl L., Sundbäck K. (2008). Benthic trophic status and nutrient fluxes in shallow-water sediments. Estuarine Coast. Shelf Sci. 78, 783–795. doi: 10.1016/j.ecss.2008.02.018

CrossRef Full Text | Google Scholar

Forrest B. M., Roberts K. L., Stevens L. M., Rabel H. (2022). “Fine scale intertidal monitoring of new river estuary: synthesis of survey data 2001-2022,” in Salt Ecology Report 100, prepared for Environment Southland. (Nelson, NZ: Salt Ecology) 47 p.

Google Scholar

Fowler D., Coyle M., Skiba U., Sutton M. A., Cape J. N., Reis S., et al. (2013). The global nitrogen cycle in the twenty-first century. Philos. Trans. R. Soc. B: Biol. Sci. 368, 20130164. doi: 10.1098/rstb.2013.0164

CrossRef Full Text | Google Scholar

Grant D. (2023) Southland places - Invercargill', Te Ara - the Encyclopedia of New Zealand. Available online at: http://www.TeAra.govt.nz/en/southland-places/page-1 (Accessed 28 August 2023).

Google Scholar

Green L., Sutula M., Fong P. (2014). How much is too much? Identifying benchmarks of adverse effects of macroalgae on the macrofauna in intertidal flats. Ecol. Appl. 24, 300–314. doi: 10.1890/13-0524.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Haddadchi A., Hicks M. (2020). Understanding the effect of catchment characteristics on suspended sediment dynamics during flood events. Hydrological Processes 34, 1558–1574. doi: 10.1002/hyp.13682

CrossRef Full Text | Google Scholar

Hale R., Depree C., Broekhuizen N. (2023). Simulating fish farm enrichment and fallowing impacts reveals unequal biogeochemical recovery of benthic variables. Aquaculture Environ. Interact. 15, 115–131. doi: 10.3354/aei00456

CrossRef Full Text | Google Scholar

Handley S. J., Horrocks M., Willis T. J., Bradley A., Bury S., Brown J., et al. (2020). Historical analyses of coastal marine sediments reveal land-based impacts on ecosystem productivity. New Z. J. Ecol. 44, 3415. doi: 10.20417/nzjecol.44.18

CrossRef Full Text | Google Scholar

Hewitt J. E., Bell R., Costello M., Cummings V., Currie K., Ellis J., et al. (2014). “Development of a national marine environment monitoring programme (MEMP) for New Zealand,” in NNew Zealand Aquatic Environment and Biodiversity Report No. 141. (Wellington, NZ: Ministry for Primary Industries).

Google Scholar

Hilgendorf H. W. (1935). The grasslands of the South Island of New Zealand: An ecological survey. DSIR Bull. 47, 24.

Google Scholar

Howard-Williams C., Davies-Colley R., Rutherford K., Wilcock R. (2011). “Diffuse pollution and freshwater degradation: New Zealand perspectives”, in Bochove van E., Vanrolleghem P. A., Chambers P. A., Thériault G., Novotná B., Burkart M. R., et al (eds) Issues and Solutions to Diffuse Pollution, selected papers from the 14th international conference of the IWA Diffuse Pollution Specialist Group, DIPCON 2010.

Google Scholar

Howarth R. W., Marino R. (2006). Nitrogen as the limiting nutrient for eutrophication in coastal marine ecosystems: Evolving views over three decades. Limnology Oceanography 51, 364–376. doi: 10.4319/lo.2006.51.1_part_2.0364

CrossRef Full Text | Google Scholar

Hume T. M., Bell R. G., De Lange W. P., Healy T. R., Hicks D. M., Kirk R. M. (1992). Coastal oceanography and sedimentology in New Zealand 1967–91. New Z. J. Mar. Freshw. Res. 26, 1–36. doi: 10.1080/00288330.1992.9516497

CrossRef Full Text | Google Scholar

Jackson J. B. C., Kirby M. X., Berger W. H., Bjorndal K. A., Botsford L. W., Bourque B. J., et al. (2001). Historical overfishing and the recent collapse of coastal ecosystems. Science 293, 629–637. doi: 10.1126/science.1059199

PubMed Abstract | CrossRef Full Text | Google Scholar

Kemp W. M., Testa J. M., Conley D. J., Gilbert D., Hagy J. D. (2009). Temporal responses of coastal hypoxia to nutrient loading and physical controls. Biogeosciences 6, 2985–3008. doi: 10.5194/bg-6-2985-2009

CrossRef Full Text | Google Scholar

Kennish M. J. (2002). Environmental threats and environmental future of estuaries. Environ. Conserv. 29, 78–107. doi: 10.1017/S0376892902000061

CrossRef Full Text | Google Scholar

Kritee K., Sigman D. M., Granger J., Ward B. B., Jayakumar A., Deutsch C. (2012). Reduced isotope fractionation by denitrification under conditions relevant to the ocean. Geochimica Cosmochimica Acta 92, 243–259. doi: 10.1016/j.gca.2012.05.020

CrossRef Full Text | Google Scholar

Landcare Research Ltd. (2013). Land Cover Database. Available at: https://lris.scinfo.org.nz/layer/104400-lcdb-v50-land-cover-database-version-50-mainland-new-zealand/.

Google Scholar

Landcare Research Ltd. (2018). New Zealand Land Resource Inventory. Available at: https://lris.scinfo.org.nz/data/?q=NZLRI.

Google Scholar

Larned S. T., Snelder T., Unwin M. J., Mcbride G. B. (2016). Water quality in New Zealand rivers: current state and trends. New Z. J. Mar. Freshw. Res. 50, 389–417. doi: 10.1080/00288330.2016.1150309

CrossRef Full Text | Google Scholar

Leathwick J. (2001). New Zealand's potential forest pattern as predicted from current species-environment relationships. New Z. J. Bot. 39, 447–464. doi: 10.1080/0028825X.2001.9512748

CrossRef Full Text | Google Scholar

Leathwick J., Mcglone M., Walker S. (2012)New Zealand's potential vegetation pattern. In: Landcare Research. Available online at: https://lris.scinfo.org.nz/layer/48289-potential-vegetation-of-new-zealand/ (Accessed 25/05/2018).

Google Scholar

Leathwick J., Overton J., Mcleod M. (2003). An environmental domain classification of New Zealand and its use as a tool for biodiversity management. Conserv. Biol. 17, 1612–1623. doi: 10.1111/j.1523-1739.2003.00469.x

CrossRef Full Text | Google Scholar

Leathwick J., Wilson G., Rutledge D., Wardle P., Morgan F., Johnston K. (2002). Land environments of New Zealand (Auckland, New Zealand: Bateman).

Google Scholar

Le Guitton M., Soetaert K., Damsté J. S. S., Middelburg J. J. (2015). Biogeochemical consequences of vertical and lateral transport of particulate organic matter in the southern North Sea: A multiproxy approach. Estuarine Coast. Shelf Sci. 165, 117–127. doi: 10.1016/j.ecss.2015.09.010

CrossRef Full Text | Google Scholar

Le Moal M., Gascuel-Odoux C., Ménesguen A., Souchon Y., Étrillard C., Levain A., et al. (2019). Eutrophication: A new wine in an old bottle? Sci. Total Environ. 651, 1–11. doi: 10.1016/j.scitotenv.2018.09.139

PubMed Abstract | CrossRef Full Text | Google Scholar

Lohrer A. M., Hewitt J. E., Thrush S. F. (2006). Assessing far-field effects of terrigenous sediment loading in the coastal marine environment. Mar. Ecol. Prog. Ser. 315, 13–18. doi: 10.3354/meps315013

CrossRef Full Text | Google Scholar

Lotze H. K., Lenihan H. S., Bourque B. J., Bradbury R. H., Cooke R. G., Kay M. C., et al. (2006). Depletion, degradation, and recovery potential of estuaries and coastal seas. Science 312, 1806–1809. doi: 10.1126/science.1128035

PubMed Abstract | CrossRef Full Text | Google Scholar

Lotze H. K., Milewski I. (2004). Two centuries of multiple human impacts and successive changes in a North Atlantic food web. Ecol. Appl. 14, 1428–1447. doi: 10.1890/03-5027

CrossRef Full Text | Google Scholar

Lotze H. K., Reise K., Worm B., Van Beusekom J., Busch M., Ehlers A., et al. (2005). Human transformations of the Wadden Sea ecosystem through time: a synthesis. Helgoland Mar. Res. 59, 84–95. doi: 10.1007/s10152-004-0209-z

CrossRef Full Text | Google Scholar

Luo J., Ledgard S. (2021). New Zealand dairy farm systems and key environmental effects. Front. Agr. Sci. Eng. 8, 148–158. doi: 10.15302/J-FASE-2020372

CrossRef Full Text | Google Scholar

Matthews K. M. (1989). Radioactive fallout in the South Pacific: a history Part 1: deposition in New Zealand (Christchurch New Zealand: National Radiation Lab).

Google Scholar

McClelland J. W., Valiela I., Michener R. H. (1997). Nitrogen-stable isotope signatures in estuarine food webs: A record of increasing urbanization in coastal watersheds. Limnology Oceanography 42, 930–937. doi: 10.4319/lo.1997.42.5.0930

CrossRef Full Text | Google Scholar

Measures R. (2016). “New River Estuary Hydrodynamic Modelling. Model build, calibration and tracer simulations,” in NIWA Client Report No: 2016077CH. (Christchurch, NZ: NIWA).

Google Scholar

Ministry for the Environment (2007). Environment New Zealand 2007 (Wellington, New Zealand: Ministry for the Environment).

Google Scholar

Ministry for the Environment (2021). New Zealand’s Greenhouse Gas Inventory 1990–2019 (Wellington, New Zealand: Ministry for the Environment).

Google Scholar

Ministry for the Environment (2023). National Policy Statement for Freshwater Management 2020 (Wellington, New Zealand: Ministry for the Environment).

Google Scholar

Ministry for the Environment, and Stats Nz. (2019). New Zealand’s Environmental Reporting Series: Environment Aotearoa 2019. Wellington, NZ: Ministry for the Environment and Stats NZ.

Google Scholar

Morrison M., Lowe M., Parsons D., Usmar N., Mcleod I. (2009). “A review of land-based effects on coastal fisheries and supporting biodiversity in New Zealand,” in New Zealand Aquatic Environment and Biodiversity Report 37 (New Zealand Ministry of Fisheries, Wellington).

Google Scholar

Needoba J. A., Sigman D. M., Harrison P. J. (2004). The mechanism of isotope fractionation during algal nitrate assimilation as illuminated by the 15n/14n of intracellular nitrate. J. Phycology 40, 517–522. doi: 10.1111/j.1529-8817.2004.03172.x

CrossRef Full Text | Google Scholar

Plew D. R., Zeldis J. R., Dudley B. D., Whitehead A. L., Stevens L. M., Robertson B. M., et al. (2020). Assessing the eutrophic susceptibility of New Zealand estuaries. Estuaries Coasts 43, 2015–2033. doi: 10.1007/s12237-020-00729-w

CrossRef Full Text | Google Scholar

Plew D. R., Zeldis J. R., Shankar U., Elliott A. H. (2018). Using simple dilution models to predict New Zealand estuarine water quality. Estuaries Coasts 41, 1643–1659. doi: 10.1007/s12237-018-0387-6

CrossRef Full Text | Google Scholar

Pratt D. R., Lohrer A. M., Pilditch C. A., Thrush S. F. (2014). Changes in ecosystem function across sedimentary gradients in estuaries. Ecosystems 17, 182–194. doi: 10.1007/s10021-013-9716-6

CrossRef Full Text | Google Scholar

Rabalais N. N., Atilla N., Normandeau C., Eugene Turner R. (2004). Ecosystem history of Mississippi River-influenced continental shelf revealed through preserved phytoplankton pigments. Mar. pollut. Bull. 49, 537–547. doi: 10.1016/j.marpolbul.2004.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts K. L., Scott-Simmonds T., Stevens L. M., Forrest B. M. (2022a). Broad Scale Intertidal Habitat Mapping of New River Estuary. Salt Ecology Report 097. (Environment Southland: Salt Ecology, Nelson, Nz).

Google Scholar

Roberts K. L., Stevens L. M., Forrest B. M., Dudley B. D., Plew D. R., Shankar U., et al. (2022b). Use of a multi-metric macroalgal index to track changes in response to nutrient loads, New River Estuary. In Salt Ecology Report 085/ NIWA Client Report 2022153CH. (Salt Ecology, Nelson, NZ: Environment Southland by Salt Ecology and NIWA under an Envirolink medium advice grant). 50p.

Google Scholar

Roberts K. L., Ward N. (2021). Proposed reference conditions in Southland estuaries: review of historical data and literature (Invercargill: Environment southland).

Google Scholar

Robertson B., Savage C. (2018). Mud-entrained macroalgae utilise porewater and overlying water column nutrients to grow in a eutrophic intertidal estuary. Biogeochemistry 139, 53–68. doi: 10.1007/s10533-018-0454-x

CrossRef Full Text | Google Scholar

Robertson B. P., Savage C., Gardner J. P. A., Robertson B. M., Stevens L. M. (2016). Optimising a widely-used coastal health index through quantitative ecological group classifications and associated thresholds. Ecol. Indic. 69, 595–605. doi: 10.1016/j.ecolind.2016.04.003

CrossRef Full Text | Google Scholar

Robertson B., Stevens L. M. (2007). Broad Scale Habitat Mapping and Sedimentation Rate (Nelson, NZ: Wriggle Ltd.).

Google Scholar

Robertson B. M., Stevens L. M. (2013a). New River Estuary Preliminary Nutrient and Sediment Load Estimates 2012/13 (Nelson, NZ: Wriggle Ltd.).

Google Scholar

Robertson B. M., Stevens L. M. (2013b). New River Estuary: fine scale monitoring of highly eutrophic arms 2012/2013 (Nelson, New Zealand: Report prepared by Wriggle Coastal Management for Environment Southland).

Google Scholar

Robertson B. M., Stevens L. M., Ward N., Robertson B. P. (2017). “Condition of Southland’s shallow, intertidal dominated estuaries in relation to eutrophication and sedimentation: output 1: data analysis and technical assessment - habitat mapping, vulnerability assessment and monitoring recommendations related to issues of eutrophication and sedimentation,” in Report prepared for Environment Southland (Wriggle Coastal Management, Nelson).

Google Scholar

Sartory D. P. (1985). The determination of algal chlorophyllous pigments by high performance liquid chromatography and spectrophotometry. Water Res. 19, 605–610. doi: 10.1016/0043-1354(85)90066-1

CrossRef Full Text | Google Scholar

Savage C. (2005). Tracing the influence of sewage nitrogen in a coastal ecosystem using stable nitrogen isotopes. Ambio 34, 145–150. doi: 10.1579/0044-7447-34.2.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Savage C., Leavitt P. R., Elmgren R. (2010). Effects of land use, urbanization, and climate variability on coastal eutrophication in the Baltic Sea. Limnology Oceanography 55, 1033–1046. doi: 10.4319/lo.2010.55.3.1033

CrossRef Full Text | Google Scholar

Snelder T. H., Larned S. T., Mcdowell R. W. (2018). Anthropogenic increases of catchment nitrogen and phosphorus loads in New Zealand. New Z. J. Mar. Freshw. Res. 52, 336–361. doi: 10.1080/00288330.2017.1393758

CrossRef Full Text | Google Scholar

Stevens L. M., Forrest B. M. (2020). Macroalgal Monitoring of New River Estuary. Salt Ecology Report 054. (Salt Ecology, Nelson, NZ: Environment Southland).

Google Scholar

Stevens L. M., Forrest B. M., Dudley B. D., Plew D. R., Zeldis J. R., Shankar U., et al. (2022). Use of a multi-metric macroalgal index to document severe eutrophication in a New Zealand estuary. N. Z. J. Mar. Freshwater Res. 56 (3), 410–429. doi: 10.1080/00288330.2022.2093226

CrossRef Full Text | Google Scholar

Sun M., Aller R., Lee C. (1991). Early diagenesis of chlorophyll-a in Long Island Sound sediments: A measure of carbon flux and particle reworking. J. Mar. Res. 49, 379–401. doi: 10.1357/002224091784995927

CrossRef Full Text | Google Scholar

Sun M.-Y., Lee C., Aller R. C. (1993). Anoxic and oxic degradation of 14C-labeled chloropigments and a 14C-labeled diatom in Long Island Sound sediments. Limnology Oceanography 38, 1438–1451. doi: 10.4319/lo.1993.38.7.1438

CrossRef Full Text | Google Scholar

Sutula M., Green L., Cicchetti G., Detenbeck N., Fong P. (2014). Thresholds of adverse effects of Macroalgal abundance and sediment organic matter on benthic habitat quality in estuarine intertidal flats. Estuaries Coasts 37, 1532–1548. doi: 10.1007/s12237-014-9796-3

CrossRef Full Text | Google Scholar

Swales A., Gibbs M. M., Handley S., Olsen G., Ovenden R., Wadhwa S., et al. (2021). “Sources of fine sediment and contribution to sedimentation in the inner Pelorus Sound/Te Hoiere”. (NIWA Client report prepared for Marlborough District Council, 2021291HN, MDC17201.). (Christchurch, NZ: NIWA).

Google Scholar

Townsend M., Lohrer D. (2015). ANZECC Guidance for Estuary Sedimentation. Prepared for Ministry for the Environment, September 2015. NIWA Client Report No HAM2015-096. (Published by Hamilton, NZ: NIWA). 45.

Google Scholar

Uusitalo L. (2007). Advantages and challenges of Bayesian networks in environmental modelling. Ecol. Model. 203, 312–318. doi: 10.1016/j.ecolmodel.2006.11.033

CrossRef Full Text | Google Scholar

Verburg P. (2007). The need to correct for the Suess effect in the application of d13C as a productivity proxy in sediment of autotrophic Lake Tanganyika in the anthropocene. J. Paleolimnology 37, 591–602. doi: 10.1007/s10933-006-9056-z

CrossRef Full Text | Google Scholar

Vieillard A. M., Newell S. E., Thrush S. F. (2020). Recovering from bias: A call for further study of underrepresented tropical and low-nutrient estuaries. J. Geophysical Research: Biogeosciences 125, e2020JG005766. doi: 10.1029/2020JG005766

CrossRef Full Text | Google Scholar

Voss M., Larsen B., Leivuori M., Vallius H. (2000). Stable isotope signals of eutrophication in Baltic Sea sediments. J. Mar. Syst. 25, 287–298. doi: 10.1016/S0924-7963(00)00022-1

CrossRef Full Text | Google Scholar

Water Framework Directive - United Kingdom Technical Advisory Group (2014). UKTAG River Assessment Method Macrophytes and Phytobenthos. Macrophytes (River LEAFPACS2). (Stirling, Scotland: United Kingdom Advisory Group) (WFD-UKTAG).

Google Scholar

Wells N. S., Eyre B. D. (2019). δ15N patterns in three subtropical estuaries show switch from nitrogen “reactors” to “pipes” with increasing degradation. Limnology Oceanography 64, 860–876. doi: 10.1002/lno.11080

CrossRef Full Text | Google Scholar

Wigand C., Roman C. T., Davey E., Stolt M., Johnson R., Hanson A., et al. (2014). Below the disappearing marshes of an urban estuary: historic nitrogen trends and soil structure. Ecol. Appl. 24, 633–649. doi: 10.1890/13-0594.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wynyard M. A. (2016). The Price of Milk: Primitive accumulation and the New Zealand Dairy Industry 1814-2014. PhD thesis. (University of Auckland).

Google Scholar

Zabarte-Maeztu I., Matheson F. E., Manley-Harris M., Davies-Colley R. J., Hawes I. (2021). Fine sediment effects on seagrasses: A global review, quantitative synthesis and multi-stressor model. Mar. Environ. Res. 171, 105480. doi: 10.1016/j.marenvres.2021.105480

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeldis J. R., Depree C., Gongol C., South P. M., Marriner A., Schiel D. R. (2020). Trophic indicators of ecological resilience in a tidal lagoon estuary following wastewater diversion and earthquake disturbance. Estuaries Coasts 43, 223–239. doi: 10.1007/s12237-019-00637-8

CrossRef Full Text | Google Scholar

Zeldis J., Hicks M., Trustrum N., Orpin A., Nodder N., Probert K., et al. (2010). “New zealand continental margin fluxes,” in Carbon and nutrient fluxes in continental margins. Eds. Liu K., Atkinson L., Quiñones R., Talaue-Mcmanus L. (Springer-Verlag, Berlin Heidelberg), 273–287.

Google Scholar

Zeldis J., Measures R., Stevens L., Matheson F., Dudley B. (2019). Remediation Options for Southland Estuaries (NIWA, Christchurch, NZ: Environment Southland).

Google Scholar

Zeldis J. R., Plew D. R. (2022). Predicting and scoring estuary ecological health using a Bayesian belief network. Front. Mar. Sci. 9. doi: 10.3389/fmars.2022.898992

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeldis J., Storey R., Plew D., Whitehead A., Madarasz-Smith A., Oliver M., et al. (2017) The New Zealand Estuary Trophic Index (ETI) Tools: Tool 3 - Assessing Estuary Trophic State using a Bayesian Belief Network. Available online at: https://shiny.niwa.co.nz/Estuaries-Screening-Tool-3/.

Google Scholar