Assessing the Eutrophic Susceptibility of New Zealand Estuaries

Abstract

We developed a method to predict the susceptibility of New Zealand estuaries to eutrophication. This method predicts macroalgae and phytoplankton responses to potential nutrient concentrations and flushing times, obtained nationally from simple dilution models, a GIS land-use model and physical estuary properties. Macroalgal response was based on an empirically derived relationship between potential nitrogen concentrations and an established macroalgal index (EQR) and phytoplankton response using an analytical growth model. Intertidal area was used to determine which primary producer was likely to lead to eutrophic conditions within estuaries. We calculated the eutrophication susceptibility of 399 New Zealand estuaries and assigned them to susceptibility bands A (lowest expected impact) to D (highest expected impact). Twenty-seven percent of New Zealand estuaries have high or very high eutrophication susceptibilities (band C or D), mostly (63% of band C and D) due to macroalgae. The physical properties of estuaries strongly influence susceptibility to macroalgae or phytoplankton blooms, and estuaries with similar physical properties cluster spatially around New Zealand’s coasts. As a result, regional patterns in susceptibility are apparent due to a combination of estuary types and land use patterns. The few areas in New Zealand with consistently low estuary eutrophication susceptibilities are either undeveloped or have estuaries with short flushing times, low intertidal area and/or minimal tidal influx. Estuaries with conditions favourable for macroalgae are most at risk. Our approach provides estuary-integrated susceptibility scores likely to be of use as a regional or national screening tool to prioritise more in-depth estuary assessments, to evaluate likely responses to altered nutrient loading regimes and assist in developing management strategies for estuaries.

Introduction

Increased input of nutrients to land has caused worldwide increases in coastal eutrophication, the process whereby a water body becomes enriched with nutrients that stimulate excessive primary production (Fowler et al. 2013; Howarth 2008; Vitousek et al. 1997). In estuaries, common responses are prolific growth of phytoplankton and opportunistic macroalgae, changes in water chemistry and reduction in biodiversity (Howarth and Marino 2006; Morand and Briand 1996). Both macroalgae and phytoplankton blooms are considered primary symptoms of eutrophication (Bricker et al. 2003), which can cause secondary symptoms including changes in sediment chemistry, reductions in water clarity, reduced dissolved oxygen, reduced invertebrate diversity and reductions in sea grass. Consequently, assessments of the susceptibility of an estuary to developing eutrophic conditions may be based on determining if, and to what degree, macroalgae or phytoplankton growth may occur.

While macroalgae and phytoplankton blooms are both responses to high nutrient loads (Fox et al. 2008; Valiela et al. 1997; Woodland et al. 2015), phytoplankton biomass is also a function of estuary flushing time (Ferreira et al. 2005). Phytoplankton blooms are not likely to occur in estuaries that have flushing times less than the phytoplankton doubling or turn-over time (Cloern 1996; Ferreira et al. 2005). In tidally dominated shallow systems, where strong vertical mixing drives oxygen replenishment thereby mitigating the effects of phytoplankton respiration, the impacts of macroalgal blooms may be of greater concern than those of phytoplankton (Sutula et al. 2014). This is particularly the case in shallow systems with large intertidal areas, high tidal exchange and/or high freshwater flow driving strong mixing or short flushing times. However, there can be localised phytoplankton blooms in poorly mixed deeper areas, or if the estuary stratifies, which may need to be considered in eutrophication assessments. Deep sub-tidal estuaries tend to have limited area suitable for macroalgae to grow, but water column stratification is more likely. Under such conditions, eutrophication driven by phytoplankton is likely more of a concern than macroalgal eutrophication (Valiela et al. 1997). Thus, the expression of primary and secondary indicators is governed by physiographic factors which modify the influence of nutrient loading and eutrophic expression in different coastal systems (Cloern 2001; Hughes et al. 2011).

The physiographic factors that influence the relative susceptibility of estuaries should therefore be factored into predictions of eutrophication risk. The concentrations of nutrients in an estuary depend on the mixing between riverine and ocean waters within the estuary, which is determined by physical characteristics including the shape and size of the estuary, tidal range and freshwater inflow (Plew et al. 2018b). Consequently, predictions of how susceptible an estuary is to eutrophication can be based on their physical properties (volume, inflow, tidal range) and nutrient load (Bricker et al. 2003; Sanderson and Coade 2010). High susceptibility does not necessarily indicate that an estuary is eutrophic because there may be factors other than nutrient load and flushing that limit algal growth, such as light, suitable substrate or available intertidal area (Valiela et al. 1997). However, estuaries with a high susceptibility may be expected to be more likely to exhibit eutrophic conditions. Susceptibility may be considered synonymous with the pressure on an estuary and is often calculated as part of an integrated pressure, state and response assessment (Borja et al. 2008; Garmendia et al. 2012). Because susceptibility can be calculated without requiring measurements of ecological state, e.g. parameters such as macroalgae or phytoplankton biomass, oxygen concentration or sediment chemistry, susceptibility calculations can be done quickly and easily, and therefore be used to prioritise more detailed investigations.

A widely used approach to predict estuary susceptibility is the Assessment of Estuarine Eutrophic Status (ASSETS) approach (Bricker et al. 1999; Bricker et al. 2003) where the physical susceptibility based on volume (dilution potential), tide and freshwater inflow (flushing potential) is combined with nutrient loads, giving a human influence factor. ASSETS was originally developed for large estuaries in the United States and overestimates the sensitivity of small estuaries to nutrient loads (Garmendia et al. 2012). A variant of this method developed for the Basque Country (WFD-BC) demonstrates better sensitivity for small estuaries due to differences in how the dilution potential is calculated (Garmendia et al. 2012). Recently, Plew et al. (2018b) used a combination of simple dilution models of mixing between fresh water and ocean waters and land-use models to predict nutrient concentrations and flushing times in New Zealand estuaries. This approach was developed to provide a screening tool for assessing the effects of catchment land use change, applicable across a range of estuary types and sizes. Because of the relationships between nutrient concentration, flushing time and algal biomass, we expected that this approach should have capability to predict ecological response with regard to macroalgae and phytoplankton.

In this paper, we present a method of predicting the susceptibility of an estuary to macroalgal and phytoplankton blooms. This susceptibility assessment approach was developed to assist New Zealand regional government manage and protect estuary health and is the basis of a freely available web-based tool (https://shiny.niwa.co.nz/Estuaries-Screening-Tool-1/, Zeldis et al. 2017a). Our method uses field observations and simple growth models to relate predicted estuary nutrient concentrations and flushing times to growth of macroalgae and/or phytoplankton. These predictions are an indicator of the susceptibility of an estuary to eutrophication. The dilution modelling approach of Plew et al. (2018b) is used to make predictions of nutrient concentrations and flushing times of New Zealand estuaries, which are in turn used to make predictions of susceptibility to eutrophication. The predicted eutrophication susceptibilities are calculated for 399 New Zealand estuaries.

Methods

Estuary Typology

For this study, we make use of two estuary typologies. The classification of an estuary does not affect our analysis but provides a useful context for interpreting patterns in susceptibility. The first is the typology used in the New Zealand Estuary Trophic Index (ETI) (Plew et al. 2018b; Robertson et al. 2016a, b). The ETI uses four primary types: coastal lakes (normally closed to the sea), shallow intertidally dominated estuaries (SIDE), shallow short residence-time tidal river estuaries (SSRTRE) and deep sub-tidally dominated estuaries (DSDE). Subtypes of SIDE and SSRTRE that intermittently close to the sea are referred to as intermittently closed and open estuaries (ICOE). The normal state of ICOEs is open, in contrast to coastal lakes which are usually or always closed. The ETI typology is less detailed than other typologies (Hume et al. 2007, 2016) but provides sufficient distinction between the major characteristics and behaviours of New Zealand estuaries (Plew et al. 2018b). In general terms, salinity is lowest for coastal lakes, followed by SSRTRE, SIDE, then DSDE. Flushing times are shortest for SSRTRE (a few hours to 1–2 days), followed by SIDE (2 days to a week) and DSDE (typically 1 week to 1 month). Coastal lakes generally have long flushing times (weeks to months).

The other typology used in this study is the New Zealand Coastal Hydrosystem classification (NZCHS) (Hume et al. 2016). This typology incorporates the ETI types but with more granularity, with 11 main classes, some of which contain subclasses. The 11 classifications span from lacustrine through to riverine, estuarine and marine systems. The names that describe the systems in the NZCHS are illustrative and are described in detail by Hume et al. (2016). The terms Waituna-type lagoon and Hāpua-type lagoon incorporate the indigenous Māori words for those systems. Waituna are lacustrine systems that are typically large (several km²), shallow (typically 2–3 m) coastal lagoons that are normally closed to the sea. Hāpua are narrow, elongated, shallow river-mouth lagoons that typically run parallel to the sea with a coarse (mixed sand and gravel) barrier beach formed by alongshore sediment transport. Hāpua are open to the sea, but while they may experience tidal backwater effects, tidal inflows seldom occur and saltwater entry mostly occurs via spray and wave overtopping during storm events (Hume et al. 2016). The relationships between the ETI and NZCHS estuary types are described in Hume (2018) and summarised in Table 1. Maps of New Zealand showing the locations of estuaries classified by each typology are given in Fig. 1.

Table 1 Relationships between New Zealand Coastal Hydrosystem (NZCHS) Classes and Estuary Trophic Index (ETI) estuary types. ETI type is that which most commonly corresponds to each NZCHS class

Fig. 1

Distribution of estuary types classified according to the a New Zealand Estuary Trophic Index (ETI) and b the New Zealand Coastal Hydrosystem (NZCHS)

The classification of estuaries is based on data available in the New Zealand Coastal Explorer Database (Hume et al. 2007) and the NZCHS (Hume et al. 2016), which draws heavily on the former. Neither of these contain information as to which systems are ICOEs. Therefore, we assume that all systems, other than coastal lakes, are open to the sea. We note that some systems, such as coastal lakes, are freshwater and not estuarine, but for the purpose of this study, all coastal hydrosystems are called “estuaries”.

Prediction of Estuary Nutrient Concentrations and Flushing Time

Plew et al. (2018b) used simple dilution models to estimate potential total nitrogen (TN), nitrate (NO₃), total phosphorus (TP), and dissolved reactive phosphorus (DRP) concentrations in, and flushing times of, estuaries and coastal water bodies across New Zealand. Potential concentrations are defined as the concentration that would occur in the absence of uptake by algae, or losses or gains due to non-conservative processes such as denitrification (Plew et al. 2018b). A dilution factor D was derived for each estuary, which allowed the potential concentration in the estuary C to be calculated from the concentration in the inflow C_R and ocean C_O

$$ C=\frac{C_R+{C}_O\left(D-1\right)}{D}. $$

(1)

The dilution factor is the reciprocal of the fraction of the estuary volume occupied by freshwater, i.e. D = 1/f = S_o /(S_o-S) where S is the volume averaged salinity of the estuary at high tide and S_o salinity of ocean water. Because measurements of S are not available for most New Zealand, we used various dilution models to calculate dilution factor, depending on the characteristics of the estuary (Plew et al. 2018b). The dilution model used for the majority of estuaries includes a tuning parameter which accounts for return flow and incomplete mixing. Plew et al. (2018b) used a regression-based approach to predict this tuning factor as a function of freshwater inflow and tidal prism. This regression is intended to be used where the tuning factor or dilution cannot be measured directly (for example, derived from measurements of salinity) and introduces a source of error, particularly for estuaries with high ratios of freshwater inflow to tidal prism. In this study, we calculate dilution for all estuaries following Plew et al. (2018b) and estimate volume-averaged estuary salinity from the dilution value. A summary of estuary properties and estimated dilution factors calculated for 399 New Zealand estuaries is available from the ETI website https://shiny.niwa.co.nz/Estuaries-Screening-Tool-1 or by request from the authors.

Concentrations of nitrogen and phosphorus in the river inflows were obtained using a GIS land use model (Elliott et al. 2016) and estuary properties from a New Zealand database (Hume et al. 2007). Land use is based on a 2008 baseline year using the New Zealand national land cover database version 3.0 (LCDB3, Landcare Research, http://www.lcdb.scinfo.org.nz), AgriBase Rural Properties database (AsureQuality, New Zealand, 2008 baseline year), and the Land Environments of New Zealand (Leathwick 2002). Ocean salinity (S_O) and nutrient concentrations (C_O) have been interpolated from the CARS 2009 climatology (CSIRO 2011). Comparisons between predicted and observed catchment nutrient loads and river concentrations are reported in previous studies (Elliott et al. 2016; Oehler and Elliott 2011; Semadeni-Davies et al. 2020a, b), while Plew et al. (2018b) compare the resulting predicted estuarine potential NO₃ and DRP concentrations with observations.

Flushing time (T_F) is related to dilution, estuary volume high tide (V) and freshwater inflow (Q_F)

$$ {T}_F=\frac{V}{DQ_F}. $$

(2)

Estuary volume at high tide was obtained from the estuary database (Hume et al. 2007), freshwater inflow from the GIS land use model and dilution calculated for each estuary as described above.

Macroalgal Susceptibility

Susceptibility to nuisance macroalgal blooms is determined using bandings developed from a comparison of potential TN concentrations with observations of macroalgae from 21 estuaries (Robertson et al. 2016b; Zeldis et al. 2017b). Macroalgae are unlikely to show phosphorus limitation for N/P molar ratios less than 30 (Atkinson and Smith 1983); consequently, we assume that nitrogen is likely to be the limiting nutrient for macroalgae growth in most estuaries. Our macroalgae bandings are based on thresholds from the Opportunistic Macroalgal Blooming Tool (OBMT) (Water Framework Directive-United Kingdom Advisory Group 2014). Macroalgal levels are assessed using Ecological Quality Rating (EQR), which is a combined metric based on both biomass and spatial measures. EQR is calculated from observations of % cover of available intertidal habitat, affected area with > 5% macroalgae cover, average biomass (wet-weight), and % cover with algae > 3 cm deep (Robertson et al. 2016b). Biomass was measured by collecting algae growing on the surface of the sediment within a defined area (e.g. 0.25 × 0.25 m quadrat) and placing it into a sieve bag. The algae were rinsed to remove sediment. Any non-algal material including stones, shells and large invertebrate fauna (e.g. crabs, shellfish) were also removed. Remaining algae were then hand squeezed until water stopped running and the wet weight of algae was recorded to the nearest 10 g using a 1-kg Pesola light-line spring scale (Robertson et al. 2016b; Water Framework Directive-United Kingdom Advisory Group 2014).

Biomass thresholds included in the OMBT were lowered for use in NZ based on unpublished data from > 25 shallow well-flushed intertidal NZ estuaries (Robertson et al. 2016b) and the results from similar estuaries in California. Sutula et al. (2014) reported that in eight Californian estuaries, macroalgal biomass of 1450 g m⁻² wet weight, total organic carbon of 1.1% and sediment total nitrogen of 0.1% were thresholds associated with anoxic conditions near the surface (aRPD < 10 mm). Green et al. (2014) reported significant and rapid negative effects on benthic invertebrate abundance and species richness at macroalgal abundances as low as 840–930 g m⁻² wet weight in two Californian estuaries. McLaughlin et al. (2014) reviewed Californian biomass thresholds and found the elimination of surface deposit feeders in the range of 700–800 g m⁻². As the Californian results were consistent with NZ findings, the latter thresholds were used to lower the OMBT good/moderate threshold from ≤ 500 to ≤ 200 g m⁻², the moderate/poor threshold from ≤ 1000 to ≤ 500 g m⁻² and the poor/bad threshold from > 3000 to > 1450 g m⁻². These thresholds are considered to provide an early warning of nutrient related impacts in NZ prior to the establishment of adverse enrichment conditions that are likely difficult to reverse. The OMBT bandings are 0–100, 100–500, 500–1000, 1000–3000, and > 3000 g wet weight m⁻² for bands of high, good, moderate, poor and bad, respectively. The modified bandings used in New Zealand are 0–100, 100–200, 200–500, 500–1450, and > 1450 g wet weight m⁻².

EQR scores range from 0 (severely impacted) to 1 (no impact). While the OMBT has 5 bandings for EQR, we combine the lowest two categories and use 4 bandings (A–D). A description of the expected ecological condition corresponding to each banding is given in Table 2.

Table 2 Macroalgal bands with corresponding EQR ratings, potential total nitrogen (TN) and potential nitrate (NO₃) ranges, and a description of expected ecological state for each band. Potential TN and NO₃ concentrations are based on annual loads and annual mean flow. Descriptions of expected ecological states are adapted from Robertson et al. (2016b)

We have developed bandings for potential TN and potential NO₃ corresponding to EQR bands by fitting linear regressions between predicted potential concentrations and observed EQR (Fig. 2, underlying data are provided in supplementary data). These regressions were used to calculate potential TN and NO₃ concentrations corresponding to EQR thresholds of 0.8, 0.6 and 0.4, which are the thresholds between A-B, B-C and C-D bands. We used observed annual nitrogen loads and annual mean flows to calculate potential concentrations for the estuaries with EQR observations (but we used modelled TN loads in subsequent susceptibility predictions). EQR observations are from peak growth (summer) periods. Our bandings therefore relate annual loads and flows to summer macroalgae response. These thresholds are reported in Table 2. The potential NO₃ thresholds calculated from the regressions are 18% lower than for TN, but the R² for the TN vs EQR and NO₃ vs EQR relationships are nearly identical (R² = 0.71).

Fig. 2

Observations of macroalgae Ecological Quality Rating plotted against a potential total nitrogen (TN) and b potential NO₃ concentrations for 21 New Zealand Estuaries. Data from Robertson et al. (2016b) and Plew et al. (2018b)

Because estuarine macroalgal growth is inhibited by low salinity conditions (Martins et al. 1999), we apply a macroalgae susceptibility band of A if the estuary salinity (calculated from the dilution modelling) is less than 5 ppt, irrespective of potential nutrient concentrations.

Phytoplankton Susceptibility

Phytoplankton biomass in an estuary is predicted using an analytical growth model which is described in the Appendix. In this model, we assume that phytoplankton growth is potentially limited by nitrogen and/or phosphorus concentration. Nitrogen is normally the limiting nutrient for phytoplankton growth in estuaries (Boynton et al. 1982; Howarth and Marino 2006). However, phosphorus limitation can occur, particularly when nitrogen loads to estuarine or coastal waters are high relative to phosphorus (Downing 1997; Harrison et al. 1990). The phytoplankton model is based on the dilution modelling approach (Plew et al. 2018b) and is similar to that of Ferreira et al. (2005) in that it applies to the estuary as a whole, and that other factors such as mortality, grazing, or sinking of phytoplankton cells are not explicitly accounted for. Nor does the model consider variability of flushing time or residence time within the estuary. Our model differs from Ferreira et al. (2005) in the way nutrient limitation is defined, and in parameterising exchange with the ocean. The model uses potential total nitrogen, total phosphorus and flushing times calculated for each estuary using the dilution modelling approach described above to predict the maximum likely phytoplankton biomass (potential phytoplankton) in the estuary. Phytoplankton biomass is calculated as chlorophyll-a (Chla) concentration as this is how phytoplankton is most commonly measured.

Because the half saturation coefficient for phosphorus is much lower than nitrogen (see Table 3), it is seldom that phosphorus will limit the rate of growth, although it may limit the maximum biomass. Figure 3 shows modelled potential chlorophyll-a concentrations as a function of flushing time and potential nitrate concentration with the default model parameters given in Table 3, for the case when phosphorus is not limiting. The model shows that phytoplankton does not accumulate for flushing times of less than 3.3 days. Above 3.3 days, the potential phytoplankton concentration increases with both flushing time and TN concentration but is largely determined by potential TN concentration for flushing times beyond 20 days.

Table 3 Parameters used in the phytoplankton model. The phytoplankton model is described in the Appendix

Fig. 3

Contours of predicted chlorophyll-a concentrations (μg/l) as a function of the potential total nitrogen concentration and estuary flushing time. The model shows chlorophyll concentrations when phosphorus is not limiting and assumes a specific growth rate µ = 0.3 day⁻¹ and a half saturation coefficient for growth response to nitrogen of 35 mg m⁻³. Phytoplankton does not accumulate in the estuary if the flushing time is shorter than the doubling time (1/µ), and steady-state concentrations are proportional to potential TN concentrations as flushing time increases

Phytoplankton blooms are generally a summer-time occurrence due to warmer water temperatures and longer daylength. While this is also true for macroalgal blooms, the macroalgal susceptibility is based on observations, and seasonality is implicitly accounted for by calibrating the prediction of macroalgal blooms using annual mean flows and nutrient loads. However, the phytoplankton susceptibility is predicted from an uncalibrated model, and therefore, we explicitly allow for seasonality in inflow as this affects both dilution (and therefore nutrient concentration) and flushing time (which is important for phytoplankton). Flows in most New Zealand rivers are below mean values during summer months (Statistics New Zealand 2017), with the majority having February (late summer) flows in the range 50–75% of mean annual flows. We assume that the inflow to estuaries is 60% of the mean annual flow when calculating potential nutrient concentrations and flushing times as inputs to the phytoplankton model. We also assume that nutrient concentrations in the inflows are the same as annual mean values (calculated from annual load and mean annual inflow). The seasonal flow adjustment has a small effect on dilution (generally increasing for most systems) and increases flushing time, while reducing potential nutrient concentrations in the estuary. A comparison between modelled Chla and observed 90th percentiles of observations is made in section “Predicted Versus Observed Phytoplankton Response” below.

The susceptibility bands for phytoplankton in estuarine systems are based on those proposed for the NZ ETI (Robertson et al. 2016b). There is limited observational data on chlorophyll-a across New Zealand estuaries, so the NZ ETI bandings for scoring estuary eutrophic state for phytoplankton are largely based on response thresholds for Basque estuaries (Revilla et al. 2010). Basque estuaries are generally shallow and well drained like the majority of New Zealand estuaries; hence, their use here rather than bandings based on US data which are more representative of deeper, less well flushed systems. The NZ ETI bandings are derived for 90th percentiles of monthly observations. The same threshold bandings are used for susceptibility using the modelled potential chlorophyll concentration (Table 4).

Many coastal hydrosystems are freshwater systems (e.g. coastal lakes or lagoons, so not strictly “estuarine”) or have low salinities that would suppress estuarine phytoplankton growth. For freshwater or brackish (salinity < 5 ppt) systems, we apply bandings from the New Zealand National Policy Statement for Freshwater Management for the maximum chlorophyll-a concentrations in lakes (Ministry for the Environment 2018).

Overall Susceptibility Banding

Primary symptoms of estuary eutrophication are high biomass of macroalgae or phytoplankton. However, these may not necessarily result in secondary symptoms of eutrophication. For example, high phytoplankton concentrations are unlikely to result in low oxygen levels or significant light attenuation if an estuary is shallow and well mixed. A high susceptibility to macroalgae blooms might not result in eutrophic conditions if there are little suitable shallow or intertidal areas available for the macroalgae to grow. Based on our experience with New Zealand estuaries, we use % intertidal area to determine whether macroalgae or phytoplankton (or both) are of concern (Table 5). The overall susceptibility of estuaries with > 40% intertidal areas is determined by the susceptibility to macroalgae, because these estuaries (mostly SIDEs or tidal lagoons) tend to have large areas favourable for macroalgal growth. The susceptibility for estuaries with < 5% intertidal area is determined from the phytoplankton susceptibility. Such estuaries have little area suitable for nuisance macroalgae and are typically large DSDE systems such as fjords, sounds or embayments where phytoplankton blooms are of most concern. For estuaries with intermediate intertidal area (5–40%), we consider both macroalgal and phytoplankton susceptibility and use the worst of these as the overall susceptibility. Our justification for this precautionary approach of taking the maximum rather than averaging (e.g. Garmendia et al. 2012) is that we consider it possible for an estuary to be considered highly eutrophic even if only one of these primary symptoms is present.

Table 4 Phytoplankton susceptibility bands based on ranges of chlorophyll-a in high (> 30 ppt), mid (5–30 ppt) and low (< 5 ppt) salinity estuaries, and a description of expected ecological state for each band adapted from Robertson et al. (2016b) and Ministry for the Environment (2018). The bandings were developed for the 90th percentile of monthly observations but applied to the model using annual N loads and 60% of mean annual flows

Validation Data

All available EQR data were used to develop the predictive model for macroalgae response. Consequently, no independent data are available for assessing the macroalgae predictor. Chlorophyll-a concentrations (Chla) for 25 estuaries were extracted from a compilation of monitoring data collected by Regional Councils (Dudley and Jones-Todd 2019). 90th percentiles of Chla concentration were calculated for each estuary. The 90th percentiles were averaged across all sites within each estuary if observations were collected at multiple locations.

Results

Predicted Versus Observed Phytoplankton Response

Observed estuary-averaged 90th percentile chlorophyll-a concentrations are positively (r = 0.47) and significantly correlated (p = 0.018) with the modelled potential Chla concentrations (Fig. 4, underlying data are provided in supplementary data). The model generally overpredicts Chla relative to the observed 90th percentiles. With regard to bandings (Table 3), the predicted Chla susceptibility band matches the observed band for 46% of estuaries, is higher than observed for 50% of estuaries and is lower for only 4% of estuaries (1 estuary). The model is expected to return higher values than observed because it predicts the maximum likely Chla concentrations under optimum conditions.

Fig. 4

Comparison of modelled potential chlorophyll maximum and observed 90th percentile chlorophyll a concentration in 25 New Zealand estuaries. The solid line shows a linear regression between predicted and observed chlorophyll-a, and the dashed line is a 1:1 slope

Macroalgae and Phytoplankton Susceptibilities

Macroalgae susceptibilities of the different NZETI estuary types are shown in the maps in Fig. 5. These maps show both the spatial distributions of estuary types, but also indicate how patterns of susceptibility to macroalgae vary around the country. These maps do not consider whether macroalgae determines overall eutrophication susceptibility for each estuary, but rather indicate the likelihood of macroalgal blooms occurring if there is suitable intertidal area available. Figure 5a shows that many SIDEs have a high (band C) or very high (band D) macroalgae susceptibility. SSRTREs around the North Island and the south-east coast of the South Island have high macroalgae susceptibility, while those along the western coast of the South Island have low or moderate susceptibilities (bands A and B), largely because salinity in, and nutrient loads to, those systems tends to be low compared to elsewhere (Fig. 5b). With few exceptions, DSDEs have low or moderate susceptibility to macroalgae (Fig. 5c), and coastal lakes have low susceptibility because they are freshwater-dominated systems (Fig. 5d). Epiphytic growths of macroalgae on freshwater macrophytes are known to occur in coastal lakes but are not considered within the current assessment.

Fig. 5

Susceptibility of New Zealand estuaries to macroalgae blooms (A = minimal, B = moderate, C = high, D = very high) for different estuary types classified according to the New Zealand Estuary Trophic Index. a SIDE = shallow intertidally dominated estuaries, b SSRTRE = shallow short residence time river estuaries, c DSDE = deep sub-tidally dominated estuaries, d coastal lakes = freshwater systems normally closed to the sea

Phytoplankton susceptibilities are mapped in Fig. 6. A surprisingly high proportion of SIDEs have high phytoplankton susceptibilities (band D, Fig. 6a), although as previously described, they are generally well mixed, so phytoplankton blooms are generally of less concern than macroalgae. SSRTEs mostly have low susceptibilities to phytoplankton (band A, Fig. 6b) due to their short flushing times. DSDEs show regional patterns in susceptibility bands. Most of the fjords along the south-west South Island have low/moderate susceptibilities, as do the sounds (deep drowned valleys) and coastal embayments at the north of the South Island.

Fig. 6

Susceptibility of New Zealand estuaries to phytoplankton blooms (A = minimal, B = moderate, C = high, D = very high) for different estuary types classified according to the New Zealand Estuary Trophic Index. a SIDE = shallow intertidally dominated estuaries, b SSRTRE = shallow short residence time river estuaries, c DSDE = deep sub-tidally dominated estuaries, d Coastal lakes = freshwater systems normally closed to the sea

As previously shown, most SSRTREs have very low phytoplankton susceptibility. Compared to nearby estuaries, a greater proportion of SIDEs on the south-east coast of New Zealand have high or very high phytoplankton susceptibility (C or D) than other parts of New Zealand. Coastal lakes with low phytoplankton susceptibility are mostly located on the west coast of the South Island and south-east coast of the North Island. While it appears that many coastal embayments on Banks Peninsula (the prominent point mid-way up the eastern coast of the South Island) have high susceptibility (band C), this is an artefact of the plotting, and the inset in Fig. 6c shows that most of the coastal embayments in this area have only moderate phytoplankton susceptibility (band B).

Predicted Eutrophication Susceptibility

The overall eutrophication susceptibility scores, which are determined from either macroalgae or phytoplankton susceptibility after considering intertidal area, are plotted for 399 New Zealand estuaries in Fig. 7. Estuaries with high or very high susceptibility (C or D band) can be found throughout the country but are particularly prevalent along the south-east coast of the South Island, and much of the North Island. Systems along the western coast of the South Island tend to have the lowest susceptibilities. Overall, 40% of systems have low susceptibility (band A), 32% have moderate susceptibility (band B), 15% have high susceptibility (C band), and 12% have very high susceptibility (D band).

Fig. 7

Overall eutrophication susceptibility of 399 New Zealand estuaries. A = minimal, B = moderate, C = high, D = very high

A breakdown of susceptibility bands by ETI estuary type (Fig. 8) shows that the susceptibility of New Zealand estuaries differs between estuary type. Eutrophication susceptibility is high or very high (C or D band) for 42% of SIDEs, for 27% of coastal lakes and 23% of SSRTEs, while only 12% of DSDE have high or very high susceptibility. SSTREs have the highest proportion of estuaries with low susceptibility (69% are A band), which is due to a combination of short flushing times preventing phytoplankton accumulation and low salinity inhibiting macroalgae growth.

Fig. 8

Susceptibility bands by a ETI type and b NZCHS type. A = minimal, B = moderate, C = high, D = very high

NZCHS types with the greatest proportions of high or very high susceptibilities are tidal lagoons (48%), Waituna-type lagoons (38%), deep drowned valleys, beach streams and tidal river mouths (26–28%). Eutrophication expression in tidal lagoons is mostly via macroalgae, while phytoplankton drives the eutrophication susceptibility in the other NZCHS types that have large proportions of high susceptibilities. No freshwater river mouths (0%) and few Hāpua-type lagoons (5%) have high susceptibility due to their low salinity and short flushing times. Only a small proportion of coastal embayments (11%) have high susceptibility, largely due to their high dilution. No fjords have high susceptibility due to a combination of high dilution and being mostly located along the south-west coast of the South Island where land cover is mostly indigenous forest.

The susceptibility factor (macroalgae or phytoplankton) that determines the overall susceptibility score for an estuary is determined by intertidal area, or by the greater of the two indicators if intertidal area is between 5 and 40%. Figure 9 shows which of the two susceptibility factors determines overall susceptibility by estuary type. If the intertidal area is between 5 and 40% and both indicators are equal, then “both” is reported. The susceptibility of ETI coastal lakes is determined by phytoplankton (Fig. 9a). SSRTREs are mostly phytoplankton-dominated, although as Fig. 8 shows, phytoplankton susceptibilities are low in most of these systems because of their short flushing times. Susceptibilities of SIDEs are nearly completely determined by macroalgae, other than in a few larger systems with intermediate intertidal areas and long flushing times where phytoplankton bandings can be high. DSDEs are mostly susceptible to phytoplankton, although some systems classified as DSDEs have high intertidal area where macroalgae can become established. More of New Zealand estuaries are susceptible to phytoplankton (207) than macroalgae (167) based on predicted nutrient concentrations and flushing times, with a further 25 systems equally susceptible to both. Of the 54 estuaries with very high eutrophication susceptibilities (band D), most (41) of these are due to very high macroalgal susceptibility scores. Macroalgae determined overall susceptibility for 63% of estuaries in bands C and D combined and co-determined for another 3%. This suggests that among New Zealand estuaries, nuisance macroalgal blooms are likely more common than phytoplankton blooms.

Fig. 9

Breakdown of which susceptibility factor (macroalgae, phytoplankton or both) determines the overall susceptibility of each estuary by a ETI type and b NZCHS type. “Both” indicates that the two susceptibility factors are equal and determine the overall susceptibility for an estuary

Figure 9b shows a similar breakdown by NZCHS type. In general, systems with low salinities (damp sand plain lakes, Waituna and Hāpua type lagoons, beach streams and freshwater river mouths) or large systems (deep drowned valleys, fjords and coastal embayments) are mostly susceptible to phytoplankton, while shallow tidally dominated systems such as tidal lagoons and shallow drowned valleys are mostly susceptible to macroalgae.

Nitrogen or Phosphorus Limitation

Predicted total N/P molar ratios are < 30:1 for 176 of the 185 estuaries for which macroalgal is the primary or co-determinant of susceptibility (Fig. 10a). The TN/TP molar ratio is < 16:1 for 185 of the 234 estuaries for which phytoplankton is the primary or co-determinant of susceptibility (Fig. 10b).

Fig. 10

Total nitrogen to phosphorus molar ratios in estuaries where susceptibility is determined by a macroalgae or b phytoplankton. Dashed vertical lines show N/P ratios of a 30 and b 16

Discussion

Seasonality

The nitrogen bandings for macroalgae susceptibility were derived from a regression fit between observed summer macroalgae EQR observations and potential TN or NO₃ concentrations calculated from annual loads and mean flows. Thus, our prediction for summer macroalgae impact is calibrated to annual loads. However, there is known to be seasonality in both flows and loads. Differences in this seasonality among estuaries may cause some of the spread in the data seen in Fig. 2, with different estuaries receiving a greater or lesser portion of nitrogen loads during summer months when growth occurs. Recent work shows that summer nutrient concentrations in terminal reaches of New Zealand rivers are lower than in winter (Whitehead et al. 2019), and lower nutrient threshold concentrations may be derived if summer loads were used. However, the catchment land-use model from which nutrient loads were obtained (CLUES) provides annual loads and does not yet provide seasonal resolution. The catchment model allows different land-use scenarios to be considered, providing a valuable management tool for assessing likely impacts on estuaries. Therefore, we developed our susceptibility bandings using potential TN or NO₃ concentrations calculated from annual loads and mean flows. A significant portion of the annual load will be delivered to estuaries during flood events and under such conditions would be rapidly discharged to the sea. Nitrogen in particulate organic matter that settles in the estuary may become mineralised during the growth season and become available to macroalgae via pore water (Robertson and Savage 2018), and the proportion of particulate nutrient load may also influence the macroalgae response relative to annual load within estuaries.

A constant adjustment to inflow was used to account for seasonality when calculating phytoplankton susceptibility (this flow adjustment was applied to the phytoplankton model but not to the macroalgal prediction because, as described above, this was calibrated to annual loads). Reducing the inflow increases the dilution factor, leading to reduced potential nutrient concentrations within estuaries, but increases flushing time. These two effects (reduced nutrient concentration and increased flushing time) have opposing influences on Chla (Fig. 3), but on net reduced predicted phytoplankton growth. Compared to using mean flows, assuming mean summer flows were 60% of annual mean flows reduced predicted Chla concentrations by 15% on average. Yet as Fig. 4 shows, the predicted Chla concentrations are higher than the 90th percentile of observations in most estuaries. We did not allow for the reduced summer nutrient concentrations in rivers (Whitehead et al. 2019), which would further reduce in-estuary potential nutrient concentrations and may bring predicted and observed Chla values into closer agreement.

Typology and Regional Effects

Estuary typology has an effect on eutrophication susceptibility, with estuaries that have high intertidal areas (ETI SIDEs or NZCHS tidal lagoons) having the most ‘high’ or ‘very high’ (C or D band, respectively) susceptibility scores (Fig. 8). Macroalgae is the primary expression of eutrophication in such systems. SIDE and tidal lagoons tend to have moderate dilution, with flushing times that can make them susceptible to phytoplankton blooms (Plew et al. 2018b), but as noted previously, they are generally shallow and well mixed so that the secondary effects of excess phytoplankton such as water column oxygen depletion are of less concern than the effects of macroalgae. The estuary types with lowest susceptibilities are systems such as Hāpua-type lagoons and freshwater river mouths (generally subclasses of SSRTRE) which, despite having poor dilution and thus high nutrient concentrations, have little intertidal area for macroalgae to grow on, low salinity which limits growth of macroalgae, and often insufficient residence times for phytoplankton to accumulate. There are, however, many observed instances in New Zealand Hāpua where stratified bottom waters become trapped in deeper holes or where mouths constrict, and localised eutrophication is evident. Overall, high or very high susceptibility scores were mostly due to macroalgae, either solely (63% of C and D bands) or in combination with phytoplankton (3%), and this proportion is higher still when considering only band D (76% due to macroalgae, 6% combined macroalgae and phytoplankton). Although more estuaries have their eutrophication susceptibility determined by phytoplankton than macroalgae, our results indicate that macroalgae blooms are likely to be a more common issue in New Zealand. Estuaries with characteristics suitable for macroalgal growth (i.e. shallow or with large intertidal areas) are most at risk of eutrophication in New Zealand.

Overlaying the differences in susceptibility due to estuary types are regional patterns in nutrient loads. Snelder et al. (2017) showed that anthropogenic increases in nutrient loads from New Zealand catchments were associated with intensive agriculture. The regions of New Zealand with the highest catchment nutrient yields are the south and east coast of the South Island, and most of the North Island. Areas with the lowest increases in nutrient loads (particularly NO₃) are the west coast of the South Island, particularly the Fiordland region, where land cover is close to natural condition. This interaction between estuary type and increases in nutrient loads from catchments is an important consideration for whether land use changes will have significant impacts on estuaries. Catchments where estuaries have short flushing times and low intertidal areas are likely to tolerate greater increases in loads than those with long flushing times and high intertidal areas.

Nutrient Limitation

Our macroalgae model assumed that nitrogen was the limiting nutrient. Macroalgae are unlikely to show phosphorus limitation for N:P molar ratios less than 30:1 (Atkinson and Smith 1983), and as Fig. 10a shows, this condition was met for 95% of estuaries for which macroalgae determines overall susceptibility. Of the 9 macroalgae susceptible estuaries for which N/P > 30:1, 7 had a D band for susceptibility and the other two are C band. For the estuaries with N/P > 30:1, we assess whether phosphorus limitation would be likely to result in a lower susceptibility band by reducing the potential N concentration to give an N/P molar ratio of 30:1, then obtaining the macroalgal susceptibility band from this modified potential N concentration. When this is done, the susceptibility bands remain the same for all the estuaries, showing that neglecting phosphorus had little effect on our macroalgal susceptibility predictions. All estuaries with an N/P ratio > 35:1 have D bands, and nutrient levels of both N and P are likely so high that neither limit growth. It is valid, therefore, to assume that nitrogen is nearly always the limiting nutrient for macroalgae growth in New Zealand estuaries. This conclusion is supported by experimental assays in New Zealand estuaries (Barr 2007; Robertson and Savage 2018).

The phytoplankton model assumed internal N/P ratios of 16:1, the typical Redfield ratio (Redfield 1958). However, nitrogen or phosphorus limitation was determined by comparing the maximum potential biomass based on either nutrient independently (eqs. 9 and 10) and this effectively determined whether phytoplankton growth was N or P limited. Because the half saturation coefficients also influenced maximum biomass, we found P limitation occurring at potential TN/TP > 21, N limitation occurring at TN/TP ratios < 18.5, with colimitation between these values. These ratios are slightly higher than the Redfield ratio. Other studies report strong nitrogen limitation at TN/TP < 20:1 (molar ratio), while P limitation consistently occurs when TN/TP > 50 (Guildford and Hecky 2000; Smith 2006). Nitrogen was the limiting nutrient for phytoplankton in the majority (81%) of New Zealand estuaries susceptible to phytoplankton.

Susceptibility Vs State

The goal of this study was to develop a method to predict the susceptibility to eutrophication based on nutrient loads. Our method consists of an empirically driven predictor of macroalgae impact, and an analytical model to predict phytoplankton biomass (expressed as chlorophyll-a). The method attempts to predict the eutrophic state of estuaries; however, the predicted state or susceptibility may not match observed state for several reasons. These include the accuracy of the catchment load model, the simplified dilution modelling approach which ignores important processes such as wind-driven circulations and spatial variability, and the accuracy of the estuary data contained in the underlying databases (volume, tidal prism, inflows, intertidal area). Sediments also influence whether or not eutrophication will occur, inhibiting growth through toxicity and/or by reducing light penetration (Lawson et al. 2007), but they can also be an important nutrient source (Kamer et al. 2004; Robertson and Savage 2018).

Inspection of Fig. 2 indicates that our regression-based predictor for macroalgae is likely to separate high and very highly impacted estuaries (C-D) from low impact (A-B), but there is no clear separation between EQR for estuaries that are classified as A or B band based on potential TN or nitrate concentration. The threshold between B and C bands indicates the most ecologically significant shift in estuary eutrophic state, from systems that are showing some signs of nutrient enrichment to those where ecosystem function is compromised. Our NO₃ threshold between B and C bands of 175 mg m⁻³ (12.5 μM) is based on annual loads and mean flows. As noted in the development of the phytoplankton model, flows in most New Zealand rivers are below mean annual values during summer months (Statistics New Zealand 2017), with the majority having February (late summer) flows in the range 0.5–0.75 of mean annual flows. Assuming loads scale in a similar ratio to flows, our summer B-C NO₃ threshold would reduce to ~ 6–9 μM, which is typical of half saturation coefficients (Runcie et al. 2003; Valiela et al. 1997; Wang et al. 2014) for macroalgae found in New Zealand estuaries, such as Ulva sp. and Gracilaria sp. (Barr et al. 2013; Robertson et al. 2016a). Thus, the B-C threshold appears consistent with nutrient concentrations that would limit growth rates of macroalgae.

The phytoplankton model appears to perform poorly against observations (Fig. 4). However, the observations are also only a representation of in-estuary conditions as samples are generally taken from the surface, at few locations and close to the shore. Thus, they seldom provide an unbiased estimate of the volume-averaged Chla concentrations in the estuary and are of limited use for calibrating or validating the phytoplankton model. There are also estuaries where the observations show low levels of Chla while the model predicts 0 values due to short flushing times. In many cases, this can be due to Chla in the river inflows, which are not accounted for in the model which assumes growth only occurs within the estuary. The model also assumes a constant maximum specific growth rate, fixed ratios of Chla to tissue nitrogen and tissue nitrogen to tissue phosphorus, and half saturation coefficients for nitrogen and phosphorus. All of these vary both between species and within species depending on conditions (Cloern et al. 1995; Eppley et al. 1969; Gibbs and Vant 1997; Laws et al. 2011a; Laws et al. 2011b; Riegman et al. 2000; Vant and Budd 1993). The model also calculates the maximum Chla that might be obtained based on the predicted nutrients (i.e. a potential Chla), and this biomass might not be obtained when the observations were collected, or if other factors not accounted for in the model limit phytoplankton growth. Despite this, the significant correlation between observed and predicted Chla indicates that the model has some predictive power, although further refinement is desirable. It is preferable for a susceptibility assessment to be precautionary and overpredict impact rather than to underpredict. Our model returned higher Chla than observed for 15 of the 25 estuaries where data were available (Fig. 4). As noted above (“Predicted Versus Observed Phytoplankton Response”), with regard to bandings, the phytoplankton model only predicted a lower banding than observed for 1 of the 25 estuaries. We note that eutrophic conditions may only express in localised areas (e.g. stratified bottom waters), thus might not be reflected in observations. While only a small portion of an estuary may be affected, it can nevertheless cause significant ecological degradation, and thus, the susceptibility assessment highlights where further investigations of eutrophic expression may be needed.

Limitations and Recommendations

Our susceptibility assessments rely on two main data sources. The first is the Coastal Explorer database (Hume et al. 2007). The data in Coastal Explorer were obtained from a variety of sources including bathymetry charts, aerial photographs, tidal models and various estuary studies (Hume et al. 2007). We have found these data to be inaccurate or outdated for many estuaries. These inaccuracies may lead to erroneous predictions of susceptibility. Also, Coastal Explorer classifies large systems as a single estuary; usually as a DSDE. Consequently, smaller SIDEs and SSRTREs that are present within these systems are missing from our analysis. This under-represents estuaries in parts of New Zealand, but also underestimates trophic state which is more likely to be elevated in localised sub-estuaries.

The Coastal Explorer database (Hume et al. 2007) and the NZCHS (Hume et al. 2016) do not clearly identify whether estuaries intermittently open and close to the sea (ICOE). We have assumed that, unless the tidal prism was 0, estuaries were open to the sea. Plew et al. (2018b) provide a means of predicting potential nutrient concentrations in ICOEs based on closure length, and that approach can be used to assess the likely eutrophic condition of ICOEs in their closed state.

The other main data source is the CLUES catchment model (Elliott et al. 2016 and references therein) which provided mean annual nutrient loads and flows. CLUES does not yet simulate groundwater or effects of irrigation on flows, or subsurface nutrient decay. Although the nutrient load models are calibrated to measured river concentrations, lag effects from land-use intensification in some catchments may result in underestimation of ultimate loadings. Modelling was based on a land cover layer that is over 10 years old (LCDB3), so may not provide load estimates that incorporate land use changes or management practices implemented since that time.

Table 5 Overall susceptibility to eutrophication is determined from macroalgal or phytoplankton based on % intertidal area

With these limitations in mind, the susceptibility assessment method developed here is most useful for (1) identifying which estuaries are likely to be at risk of eutrophication (i.e. a regional or national screening tool such as Plew et al. (2018a)), (2) investigating how estuaries may respond to changes in nutrient loads and (3) identifying whether macroalgae, phytoplankton or both are the likely primary response to nutrient loads. Within New Zealand, Regional Councils are generally responsible for making decisions that affect the health of estuaries. Our susceptibility assessment method can be used by Regional Councils to prioritise monitoring, to assist in developing catchment management plans that provide desired outcomes for estuary health (i.e. nutrient load limits), and as an indicator of the likely magnitude of impacts of activities that affect nutrient loading to estuaries. Here, we used a catchment model to estimate nutrient loads to estuaries on a national scale. This catchment model can be used to calculate nutrient loads under different land use scenarios. Alternatively, measured loads could be used where available. Also, the likely impacts on estuary health of the addition or removal of point sources, such as waste water discharges, can be considered via the effect on susceptibility. With regard to impacts on estuary health, we consider our approach best suited as first-order screening assessment that might trigger more detailed investigations if changes in susceptibility are significant. More accurate predictions of susceptibility or state for individual estuaries can be obtained by improving the performance of the underlying dilution models by calibrating these with estuary specific data such as bathymetry, tidal prism and salinity (Plew et al. 2018b), using observed nutrient loads and inflows, and making comparison with observed conditions within the estuaries. The dilution modelling approach used here provides steady-state, estuary-averaged assessments (i.e. no temporal or spatial resolution). More complex approaches can include coupled hydrodynamic-ecological models (e.g. Cerco and Noel 2013; del Barrio et al. 2014; Testa et al. 2013), although the high resolution and detail of such an approach comes at the cost of requiring considerable input data, complex process parameterisation, computational time, calibration, validation and modeller skill (Ganju et al. 2016); making them difficult to apply on a regional or national scale.

Appendix: Phytoplankton model

The rate of change of phytoplankton in the estuary is a balance between growth rate and advection. Averaged over a tidal period T,

$$ V\frac{dM}{dt}=Q{M}_{in}-\frac{\left(P+ QT/2\right)}{T}M+\frac{\left(P- QT/2\right)}{T} BM+\frac{\left(P- QT/2\right)}{T}\left(1-B\right){M}_o+\mu MV\frac{N}{N+{N}_s}\frac{K}{K+{K}_s}. $$

(3)

M is the phytoplankton concentration in the estuary, measured in terms of nitrogen (g N m⁻³), M_in is the phytoplankton concentration in the freshwater inflow, which we assume to be 0 for estuarine species. M_O is the phytoplankton concentration in the coastal region where the estuary is located, Q (m³ day⁻¹) the freshwater inflow, T (d) the tidal period, P the tidal prism (m³), V the estuary volume at high tide (m³), B the tuning factor that corrects for return flow and incomplete mixing (Plew et al. 2018b), N and K the nitrogen and phosphorus concentration in the estuary, N_S and K_s are half saturation coefficients for nitrogen and phosphorus limited growth respectively, and μ a net specific growth rate (day⁻¹) that include both growth and mortality. The term on the left of eq. (3) is the rate of change of phytoplankton biomass in the estuary. The terms on the right hand side account for, in order: the input of phytoplankton from the freshwater sources, discharge of phytoplankton to the ocean on the outgoing tide, return flow of phytoplankton from the ocean back to estuary (some of the phytoplankton discharged on the outgoing tide returns to the estuary on the subsequent ebb tide), input of oceanic phytoplankton on the incoming tide and the growth of phytoplankton within the estuary. The final term in eq. (3) has the effect of reducing the specific growth rate at low concentrations of N and P. Equation (3) (and 4–5) assume that the length of the ebb and flood tides are approximately equal, i.e. T/2.

A similar mass balance applies to nitrogen

$$ V\frac{dN}{dt}=Q{N}_{in}-\frac{\left(P+ QT/2\right)}{T}N+\frac{\left(P- QT/2\right)}{T} BN+\frac{\left(P- QT/2\right)}{T}\left(1-B\right){N}_o-\mu MV\frac{N}{N+{N}_s}\frac{K}{K+{K}_s} $$

(4)

where N_in is the nitrogen concentration in the freshwater inflow and N_o the nitrogen concentration in the ocean.

And to phosphorus

$$ V\frac{dK}{dt}=Q{K}_{in}-\frac{\left(P+ QT/2\right)}{T}K+\frac{\left(P- QT/2\right)}{T} BK+\frac{\left(P- QT/2\right)}{T}\left(1-B\right){K}_o-\mu \beta MV\frac{N}{N+{N}_s}\frac{K}{K+{K}_s}. $$

(5)

K, K_o and K_in are the phosphorus concentrations in the estuary, ocean and freshwater inflow, respectively, and β is the ratio of tissue phosphorus to nitrogen in the algae, which is assumed to be constant. The widely used Redfield ratio of N/P = 16:1 (molar) is assumed here, which equates to β = 0.138 g P/g N.

Equations 3–5 are solved simultaneously for steady-state conditions (i.e. dN/dt = dM/dt = dK/dt = 0), making use of eq. (1) to define the potential nitrogen and phosphorus concentrations N_p and K_p and substituting the tuning factor B with the dilution factor D as follows (Plew et al. 2018b).

$$ B=\frac{P+\frac{QT}{2}-D}{P-\frac{QT}{2}}. $$

(6)

Because we are interested in estuarine phytoplankton, we assume that the concentration of phytoplankton in the ocean can be neglected, i.e. M_o = 0. The solution is obtained by solving the following quadratic

$$ M=\frac{-b\pm \sqrt{b^2-4 ac}}{2a}. $$

(7)

where the coefficients are

$$ a=\beta \left(u{T}_F-1\right) $$

$$ b={K}_s+{N}_s\beta -\left({N}_P\beta +{K}_p\right)\left(u{T}_F-1\right) $$

$$ c={N}_P{K}_p\left(u{T}_F-1\right)-{N}_P{K}_s-{N}_s{K}_p-{N}_s{K}_s $$

A valid solution is obtained when the flushing time is sufficient that phytoplankton growth exceeds the net flux of phytoplankton out of the estuary. This minimum flushing time is calculated as

$$ {T}_F>\frac{1}{\mu}\max \left(1+\frac{N_s}{N_p},1+\frac{K_s}{K_p}\right). $$

(8)

Else, M = 0.

A further restriction is placed on the model that the nutrient concentrations in the estuary cannot be negative. This restriction is imposed by calculating the maximum biomass that could be obtained if only one nutrient (either nitrogen or phosphorus) limits growth. If only nitrogen or phosphorus limits growth, then a solution for biomass can be obtained as

$$ {M}_N={N}_p-\frac{N_s}{u{T}_F-1}. $$

(9)

Or

$$ {M}_P=\frac{1}{\beta}\left({K}_P-\frac{K_s}{u{T}_F-1}\right). $$

(10)

The phytoplankton model then reduces to

$$ Chla=\frac{1}{\alpha}\min \left(M,{M}_N,{M}_P\right)\ \mathrm{for}\ {T}_F>\frac{1}{\mu}\max \left(1+\frac{N_s}{N_p},1+\frac{K_s}{K_p}\right) $$

(11)

Else

$$ Chla=0 $$

The coefficient α is the ratio between chlorophyll-a and tissue nitrogen concentration, which is assumed constant.

The values of the various parameters and their sources are given in Table 3. A typical nitrogen half-saturation coefficient for phytoplankton growth is N_s = 35 mg m⁻³ (Eppley et al. 1969). The half saturation for phosphorus is much lower, typically 0.1 mg m⁻³ (Laws et al. 2011a, b; Riegman et al. 2000). Reported phytoplankton growth rates vary widely, depending on phytoplankton species, and factors including (but not limited to) temperature, light and salinity (Alpine and Cloern 1988; Eppley 1972; Grzebyk and Berland 1996). Net growth rates of < 0.2–0.4 day⁻¹ have been reported in New Zealand estuaries (Gibbs and Vant 1997; Vant and Budd 1993). We adopt a mid-range value of k = 0.3 day⁻¹.

Phytoplankton concentrations are converted to Chla concentrations (μg/l) using the standard Redfield ratio (Redfield 1958), which for phytoplankton is approximately 106:16 for carbon/nitrogen (molecular ratio). Carbon to Chla ratios vary but a typical value is 50 μg C l⁻¹:1 μg Chla l⁻¹ (Cloern et al. 1995). Therefore, the conversion between phytoplankton nitrogen and Chla is α = 8.8 μg N μg⁻¹ Chla.