Reference and current Trophic Level Index of New Zealand lakes: benchmarks to inform lake management and assessment

ABSTRACT

Knowledge of trophic status is fundamental to understanding the condition and function of lake ecosystems. We developed regression models to predict chlorophyll a concentrations (chl a) in New Zealand lakes for reference and current states, based on an existing dataset of total nitrogen (TN) and total phosphorus (TP) concentrations for 1031 lakes. Models were then developed to predict Secchi depth based on chl a and a sediment resuspension term applicable to shallow lakes. Estimates of all four Trophic Level Index (TLI) variables (chl a, TN, TP and Secchi depth) were analysed to estimate reference and current state TLI for the nationally representative sample of 1031 lakes. There was a trend of eutrophication between reference and current states, with systematic differences among lake geomorphic types. Mean chl a increased 3.5-fold (2.42 mg m⁻³ vs. 8.32 mg m⁻³) and mean Secchi depth decreased (indicating lower clarity) by approximately one-third (9.62 m vs. 6.48 m) between reference and current states. On average, TLI increased by 0.67, with the TLI increase >1 in approximately one-third (31%) of lakes. This study informs the status of lake ecosystems in NZ and provides benchmarks to guide management and assessment.

KEYWORDS:

• Baseline
• chlorophyll
• eutrophication
• nutrients
• Secchi depth
• trophic status
• water quality

Introduction

Knowledge of trophic status is fundamental to understanding the condition and function of lake ecosystems. Lake trophic status can be considered a ‘master variable’ that corresponds to important lake attributes such as the prevalence of harmful cyanobacteria blooms, depletion rate of deep-water oxygen, biomass of aquatic vascular plants, invertebrate and vertebrate species composition, and perceived aesthetic values (Smith 2009; Angradi et al. 2018). In New Zealand (NZ), multiple lake characteristics have been shown to correlate with lake trophic status, including denitrification potential (Bruesewitz et al. 2011), bottom-sediment accumulation rate (Trolle et al. 2008; Santoso et al. 2017), CO₂ efflux (Santoso et al. 2017), prokaryotic picoplankton density (Burns and Stockner 1991), crustacean zooplankton density (Chapman et al. 1985), rotifer assemblage composition (Duggan et al. 2001), food chain length (Kelly and Schallenberg 2019), the mass and quality of consumer trophic resources (Collier et al. 2018), and the distribution of kōura (Paranephrops planifrons) (Kusabs et al. 2015).

In NZ, lake trophic status is commonly measured using the Trophic Level Index (TLI), which comprises a value of 0 to ∼7 that relates to seven classes ranging from ultra-microtrophic to hypertrophic (Burns et al. 1999). The TLI is calculated using annual average measurements of chlorophyll a concentration (chl a), Secchi depth (Secchi), and concentrations of total nitrogen (TN) and total phosphorus (TP). The TLI is analogous to the Trophic State Index (TSI; Carlson 1977) applied in the USA although, unlike the TSI, the TLI includes nitrogen concentrations in the calculation. This inclusion reflects that nitrogen, in addition to phosphorus, frequently limits phytoplankton biomass accumulation in NZ lakes, with nitrogen limitation seemingly more prevalent in NZ than in the northern hemisphere (White 1983; Abell et al. 2010). The TLI is routinely measured by NZ regional councils, with measurements also used to characterise state and trends of lake water quality at a national scale (Gluckman 2017; Ministry for the Environment & Stats NZ 2019). National legislation requires objectives to be set for three of the four constituent TLI variables in lakes (New Zealand Government 2017), and some regional agencies have enacted statutory rules that require establishment of TLI targets for important lakes, with changes to catchment management required if targets are not met (BOPRC 2008).

The four constituent TLI variables are highly inter-correlated, although the relationships between variables varies among lakes depending on lake characteristics (Burns et al. 1999). Chl a is strongly controlled by nutrient concentrations in temperate lakes and is commonly predicted based on TN and/or TP concentrations, although the strength and form of the relationship varies depending on factors that include nutrient limitation status and top-down controls such as zooplankton grazing (Abell et al. 2012). Secchi is inversely correlated with chl a, although the relationship between these variables is affected by the presence of other light attenuating substances (Davies-Colley and Smith 2001). In particular, Secchi is often lower in shallow lakes due to reduced clarity caused by wind driven resuspension of fine sediments (Hamilton and Mitchell 1997).

Many lakes in NZ have been affected by anthropogenic (cultural) eutrophication, with nutrient and sediment pollution from pastoral agriculture identified as a major causal factor (Gluckman 2017). Quantifying and managing human influences on lake trophic status requires knowledge of natural variability in trophic status among lakes, thereby allowing the magnitude of anthropogenic eutrophication to be quantified and providing context for setting management objectives and targets. Lake trophic status varies naturally due to a range of catchment and within-lake factors. At the catchment scale, geology can affect lake trophic status (Dillon and Kirchner 1975). In NZ, there is a geological influence associated with elevated contributions of phosphorus from areas with volcanic geology such as the Taupō volcanic zone (Timperley 1983), whereas streams that drain areas with alluvium or sedimentary geology have naturally lower phosphorus concentrations (McDowell et al. 2013).

Catchment climate and hydrology are also important factors that influence lake trophic status by regulating the erosion and transport of nutrients and light attenuating substances (Lohse et al. 2009). At the lake scale, basin morphometry affects lake trophic status, and morphometric variables can be inter-correlated (Håkanson 2005). Hydraulic residence time (τ_w; influenced by morphometry and hydrology) affects particle settling and flushing losses of plankton, with τ_w inversely correlated with lake trophic status (Vollenweider 1971). Shallower lakes tend to have lower τ_w (Håkanson 2005) and are naturally susceptible to sediment resuspension by wind waves, which can increase trophic status by mobilising nutrients and light attenuating sediments (Hamilton and Mitchell 1997). Lake fetch (correlated with area) affects the mixing depth (Davies-Colley 1988; Hanna 1990) and the potential for wind driven resuspension of sediments (Hamilton and Mitchell 1997); thus, sediment resuspension is typicaly greatest in large shallow lakes, especially those that are elongated with the longest axis aligned parallel to the dominant wind direction.

Reference states relating to pristine conditions in the absence of human disturbance provide benchmarks against which to assess environmental impacts and inform management (Hawkins et al. 2010). Due to the influences of catchment and within-lake factors, reference conditions for lake trophic status variables should ideally be developed for lakes on an individual basis, although reference conditions for similar lakes can be grouped to develop ranges of values for lake categories (e.g. Carvalho et al. 2008; Poikāne et al. 2010). Recently, Abell et al. (2019) used a modelling approach to estimate annual mean TN and TP concentrations for the surface waters of ∼1030 NZ lakes, with concentrations estimated separately for reference and current states. Their approach built on work by McDowell et al. (2013), who derived median reference concentrations for nutrients in categories of NZ streams and rivers that accounted for climate, topography and geology. This previous work provided an opportunity to extend existing research to predict other trophic status variables in a nationally representative sample of lakes across NZ.

The objective of this study was to use the TN and TP concentrations estimated by Abell et al. (2019), alongside other variables, to predict reference and current values of chl a and Secchi, which could then be used to estimate reference and current TLI values for a sample of >1000 NZ lakes. We hypothesised that predicting reference and current lake trophic status for individual lake types could provide benchmarks to guide lake management and enhance understanding of the magnitude of anthropogenic changes to lake water quality in NZ.

Methods

Overview of analysis sequence and model fitting

Analysis steps are summarised in Figure 1 and described further in the sections below. All analyses and data manipulation were conducted using R (R Core Team 2017). Linear models (lm function) were fitted to predict chl a and Secchi, with the following metrics used to compare performance when trialling models: standard error of the estimate, sum of squared errors (SSE), R², Akaike’s Information Criterion (AIC; Burnham and Anderson 2002) and the root-mean-square error (RMSE) of a 10-fold cross validation. Overall, AIC (a measure of model fit) and RMSE (a measure of model error) were the primary performance metrics used to select optimum models. Bias was examined by visually inspecting observed and predicted values plotted alongside a 1:1 line.

Figure 1. Summary of analysis steps to predict current and reference Trophic Level Index (TLI) values for 1031 lakes in New Zealand. Acronyms are defined in the text.

A detailed residual analysis was undertaken to check that model assumptions were met and to confirm the necessity of transformations (required for most variables). Residual analysis involved reviewing a range of diagnostic plots (e.g. Q-Q plots) to examine whether residuals were normally distributed, and violin plots to compare residuals among depth and geomorphic classes (hypothesised covariates). Residual plots and detailed results of model evaluation are provided in Supplementary Information.

National Water Quality Monitoring Network dataset

Lake water quality data (Table 1) were obtained from the National Water Quality Monitoring Network (NWQMN) database compiled by NZ regional government agencies. Data were provided via the Land, Air, Water Aotearoa partnership, with additional information about lake geomorphic type and morphometric variables, including maximum depth (z_max) and fetch, obtained from the Freshwater Ecosystems of New Zealand (FENZ) database (Leathwick et al. 2010). Measurements from the NWQMN were collected from the surface or within the epilimnion in stratified lakes, with sampling typically undertaken monthly (Burns et al. 2000; Davies-Colley et al. 2012). The NWQMN database (1976–2014) was first screened, primarily to remove measurements collected prior to 2000 and to identify lakes with ≥10 measurements of TN, TP and chlorophyll a (see Abell et al. (2019) for details of additional screening criteria). This sample was used to develop models to predict chl a (n = 73; Table 1), and then screened further to identify lakes with ≥10 measurements of Secchi (not available for 30 of the lakes) to create a sample to develop models to predict Secchi (n = 43; Table 1).

Table 1. Summary of NWQMN lake water quality data used to develop models to predict annual mean chlorophyll a (chl a) concentrations and annual mean Secchi depth (Secchi).

Dependent variable modelled	Variable	Units	Number of lakes	Number of measurements per lake (mean in parentheses)	Minimum of annual mean values	Maximum of annual mean values	Mean of annual mean values	Mean coefficient of variation (%)
Chl a	TN	mg m^−3	73	10–261 (66)	39	3266	567	45
	TP	mg m^−3	73	10–261 (68)	4.7	349.0	44.8	93
	Chl a	mg m^−3	73	10–262 (68)	0.41	184.59	15.37	100
Ref_ Chl a	TN	mg m^−3	66	10–261 (65)	39	1440	367	38
	TP	mg m^−3	66	10–261 (66)	4.7	145.9	22.9	38
	Chl a	mg m^−3	66	10–262 (66)	0.41	44.54	6.67	100
Secchi	Chl a	mg m^−3	43	10–45 (34)	0.41	184.59	12.22	95
	Secchi	m	43	10–261 (79)	0.11	15.28	4.83	36
Ref_Secchi	Chl a	mg m^−3	16	51–262 (129)	1.03	17.56	3.93	85
	Secchi	m	16	48–261 (149)	2.64	15.28	8.26	25

Note: ‘Ref_’ denotes reference values.

Reference and current nutrient concentrations

We used the existing dataset of modelled TN and TP concentrations in NZ lake surface waters developed during a previous study by Abell et al. (2019), which we briefly summarise in this section. The authors used nutrient mass loading models (‘box models’) to estimate annual average lake nutrient concentrations that corresponded to reference (REF_TN_lake and REF_TP_lake) and current states (TN_lake and TP_lake). Estimates of both sets of TN and TP concentrations were derived for a sample of 1031 lakes included in the Waters of National Importance (WONI) database (Figure 2), which comprises spatial, land use and physical data for all 3820 mapped lakes in NZ with area >1 ha (Snelder 2006; Leathwick et al. 2010).

Figure 2. Locations and geomorphic type of the 1031 lakes for which current and reference Trophic Level Index values were estimated in this study. Adapted from Abell et al. (2019).

Models were developed using a training dataset comprising measurements for lakes included in the NWQMN database. Following a model screening process, models were selected to predict TN and TP concentrations based on lake inflow concentrations and morphometric variables (τ_w for TP and maximum depth (z_max) for TN). To predict current concentrations, lake inflow concentrations were estimated using the Catchment Land Use for Environmental Sustainability (CLUES) model (Elliott et al. 2016). To predict reference concentrations, lake inflow concentrations were derived from the stream reference concentrations estimated by McDowell et al. (2013) and streamflow modelled with CLUES.

Chlorophyll a concentrations

Current chlorophyll a concentrations

A model to predict chl a was fitted using data in the NWQMN dataset (Table 1). Chl a was strongly correlated with TN (Spearman’s r = 0.91) and TP (Spearman’s r = 0.90). Based on these relationships, the following combinations of predictor variables were trialled: (1) TN; (2) TP; (3) TN + TP, and; (4) TN + TP + TN·TP, with the latter model including an interaction term between TN and TP. It was necessary to log₁₀-transform all variables to meet linear model assumptions. Based on performance, the best chl a model included the terms TN and TP (the third model listed above), with RMSE of 0.21 mg m⁻³ (in transformed units) and r² of 0.89 (Table 2). Model performance (AIC) could not be improved by adding depth (z_max) or geomorphic type and a residual analysis confirmed that the residuals were not correlated with these hypothesised covariates.

Table 2. Summary of models used to predict current and reference state values of chlorophyll a and Secchi depth in 1031 New Zealand lakes.

Dependent variable	Model	SEE	RMSE	R^2	Notes
log10chl a (mg m^−3)	$-1.80+0.70\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}T{N}_{lake}+0.55\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}T{P}_{lake}$	0.21	0.21	0.89	chl a, mean lake chlorophyll a concentration in mg m^−3 TNlake (or TPlake), mean in-lake surface TN (or TP) concentration in mg m^−3
log10Ref_Chl a (mg m^−3)	$-1.73+0.65\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}Ref\mathrm{\_}T{N}_{lake}+0.59\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}Ref\mathrm{\_}T{P}_{lake}$	0.20	0.21	0.84	Ref_, denotes reference state values
Secchi^0.5 (m)	$3.46-0.74\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}Chla-0.79\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}Chla\cdot d-0.35\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}\left({\displaystyle \frac{Fetch\cdot U^2}{z_{max}}}\right)\cdot \left(1-d\right)$	0.32	1.68	0.89	Secchi, mean lake Secchi depth in m, Fetch, maximum lake fetch in m, U, mean windspeed at lake (m s^−1), zmax, maximum lake depth in m, d, equal to 0 if zmax<20 m or equal to 1 if zmax ≥ 20 m
Ref_Secchi^0.5 (m)	$3.42-1.46\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}Ref\mathrm{\_}Chla$	0.35	2.48	0.70

Note: SEE: standard error of the estimate; RMSE: root mean square error of a 10-fold cross validation.

Confidence intervals were calculated for the chl a model estimates with a parametric bootstrap (n = 1000) that accounted for the propagation of errors in multiple models (note that, for unmonitored lakes, chl a was modelled based on nutrient concentrations that were in turn modelled based on other variables). On each iteration i of the bootstrap, new values of chl a ($chl{a}_i^{\ast }$) were drawn randomly from a normal distribution with mean equal to chl a estimated with the fitted model and standard deviation equal to the standard error of the estimate of the fitted model. We then refit the original model to $chl{a}_i^{\ast }$, adding model uncertainty to the coefficient estimates (${\text{β}}_i^{\ast }$). Finally, we predicted chl a* values for each iteration of the bootstrap as:(1) $chl{a}_i^{\ast }={\text{β}}_{0,i}^{\ast }+{\text{β}}_{1,i}^{\ast }T{N}_{lake,i}^{\ast }+{\text{β}}_{2,i}^{\ast }T{P}_{lake,i}^{\ast }$(1) where $T{N}_{lake,i}^{\ast }$ and $T{P}_{lake,i}^{\ast }$ are the bootstrapped predicted TN_lake and TP_lake values for iteration i from the previous bootstrap exercise. Specifically, values of $T{N}_{lake,i}^{\ast }$ and $T{P}_{lake,i}^{\ast }$ were recalculated to account for estimated error in the lake inflow nutrient concentrations used by Abell et al. (2019) by recalculating TN_lake and TP_lake using lake inflow nutrient concentrations that were adjusted by adding a random error term based on the standard error of lake inflow nutrient concentrations (see Abell et al. 2019). Each iteration of the bootstrap therefore accounted for uncertainty in the model (e.g. structural error) through the refitting of the model coefficients, as well as the uncertainty in the predicted lake inflow nutrient concentration (errors in variables) that were used to predict $T{N}_{lake}^{\ast }$ and $T{P}_{lake}^{\ast }$ during the previous step. The ranges of the 95% confidence intervals were then defined based on the 2.5% and 97.5% quantiles of the 1000 bootstrapped iterations.

Reference chlorophyll a concentrations

A model to predict annual mean reference chl a concentrations (Ref_chl a) was fitted using the same approach as described above except a smaller (n = 66 vs. 73) NWQMN dataset was used. This smaller sample was created by omitting lakes with estimates of TN_lake and TP_lake that were higher than the upper bounds of the 95% confidence intervals of Ref_TN_lake and Ref_TP_lake estimates for the sample of 1031 lakes; i.e. it was considered logical to omit lakes with current nutrient concentrations greater than those expected in NZ lakes under a reference state when predicting Ref_chl a. Model performance and parameter coefficients for this model were highly similar to the model used to predict current chl a concentrations (Table 2). Confidence intervals were calculated for Ref_chl a estimates with a parametric bootstrap, as described above for chl a.

Secchi depth

Current Secchi depth

Models were trialled to estimate mean Secchi (m). Chl a was included as a predictor variable in candidate models as there is generally a strong correlation between Secchi and chl a in temperate lakes (Vant and Davies-Colley 1984; Abell et al. 2012). Concentrations of non-algal suspended sediment and coloured dissolved organic matter influence Secchi (Vant and Davies-Colley 1984; Swift et al. 2006). Resuspension of non-algal sediment is expected to increase in significance with declining depth in shallow lakes. For deep lakes, where the depth of most of the bed exceeds a threshold related to potential wave height, wind-driven resuspension is expected to be negligible, although we recognise that high levels of light-attenuating fine sediment (inorganic glacial flour) occur in several glacial lakes (Rose et al. 2014), which may be an exception to this. The influence of non-algal suspended sediment on Secchi in shallow lakes is generally expected to be greater in lakes under current conditions than under a reference state because many shallow lakes in NZ have undergone a unidirectional regime shift involving increases in turbidity (Schallenberg and Sorell 2009).

Examination of the relationship between Secchi and chl a showed no evidence of a systematic difference in the relationship among lake geomorphic types, but there was a systematic difference between shallow and deep lakes; in general, the Secchi value at a given chl a value was higher (i.e. greater clarity) in deep lakes compared with shallow lakes. To represent the potential influence of sediment resuspension in shallow lakes, we trialled including $z_{max}$ as a predictor variable, in addition to the following variables that we hypothesised would be positively correlated with wind wave height: surface area in ha (Area), maximum lake fetch in m (Fetch) and the dynamic ratio calculated as ${\displaystyle \frac{Are{a}^{0.5}}{z_{max}}}$ (Håkanson 1982). Further, we also trialled inclusion of a ‘resuspension’ term, defined as${\displaystyle \frac{Fetch\cdot U^2}{z_{max}}}$ where U is mean windspeed (m s⁻¹) at each lake based on analysis of regional climate data by Leathwick et al. (2010). The resuspension term, with units of m² s⁻² (equivalent to J kg⁻¹; the energy provided relative to weight), was hypothesised to be proportional to the mean wind wave height in a lake.

Unlike wind-driven resuspension, it was not possible to include continuous variables in the candidate models to represent the effect of coloured dissolved organic matter on Secchi. To consider the potential influence of this factor, we instead analysed the residual values from a Secchi = f(chl a) linear regression model to confirm whether peat lakes and/or riverine lakes were outliers based on a hypothesis that organic-rich soils in catchments of these lakes could lead to high coloured dissolved organic matter concentrations (Vant and Davies-Colley 1984). We also examined whether glacial lakes were outliers, which we hypothesised would be the case due to the influence of fine glacial sediments on light scattering. However, this data exploration did not indicate that lakes with these geomorphic types were outliers, although the analysis was constrained by small sample sizes (e.g. only four peat lakes were included). Based on this, the potential effect of coloured dissolved organic matter or glacial sediments on Secchi was not modelled and therefore this factor is implicitly included in the model error (quantified using bootstrapping, as described below).

Residual analysis and a Box Cox power test (Box and Cox 1964) showed that linear model assumptions (Quinn and Keough 2002) were best met when Secchi was square-root transformed; i.e. Secchi^0.5 was the dependent variable in the models that were trialled. Comparison of model performance among candidate models that included the predictor variables described above showed that the best performing model included chl a and the resuspension term $\left({\displaystyle \frac{Fetch\cdot U^2}{z_{max}}}\right)$ as predictor variables. This model was therefore applied to predict Secchi. This model was piecewise as the fitted parameter values varied depending on whether a lake is shallow or deep, with the resuspension term only used for shallow lakes (Table 2). Shallow (n = 466) and deep (n = 565) lakes were distinguished based on a threshold of z_max = 20 m, which was determined iteratively by trailing a range of thresholds and selecting the threshold that optimised model performance. Predictions of Secchi^0.5 had RMSE of 1.68 m and r² of 0.89 (Table 2).

For both Secchi and Ref_Secchi (see below), model uncertainty was quantified by refitting the data to response variables drawn from a normal distribution with variance equal to the variance of the estimated values. Confidence intervals for model predictions were quantified using bootstrapping methods consistent with those described for chl a. For each lake, developing these predictions entailed recalculating Secchi model coefficients based on bootstrapped values (Secchi*) randomly drawn from a normal distribution centred on Secchi, with a standard deviation equal to the standard error of the estimate for the model. Bootstrapped chl a* values calculated in the previous step using Equation 1 were then used in the model shown in Table 2. 95% confidence intervals were calculated as the 2.5% and 97.5% quantiles of 1000 bootstrapped iterations.

Reference Secchi depth

To predict reference condition Secchi depth (Ref_Secchi), we assumed that chl a was the sole determinant of Ref_Secchi. We recognise that sediment resuspension would also presumably influence Ref_Secchi in shallow lakes; however, we expect that the magnitude of this effect would be lower as the presence of healthy macrophyte communities would stabilise sediments. Also, sediment loads would generally be lower, and benthivorous non-indigenous fish such as common carp (Cyprinis carpio) would be absent, under a reference state. For all lakes, Ref_Secchi was therefore estimated using a linear regression with chl a as the sole predictor variable. The model was fitted to a sub-sample comprising only deep lakes (defined based on z_max > 20 m, as described above) and with chl a ≤ 33 mg m⁻³, which corresponds to the upper bound of the confidence intervals fitted to the Ref_Chl a estimates. Ref_Secchi was square root transformed during model fitting. The resulting regression model has similar predictive performance and fitted parameter values as the model described above to predict Secchi in deep lakes (Table 2).

TLI predictions

The optimum models to predict chl a and Ref_Chl a were applied with the current and reference TN and TP concentrations for 1031 lakes estimated by Abell et al. (2019). The transformed model predictions of chl a and Ref_Chl a were back-transformed, with a correction factor used where relevant to convert log₁₀-transformed variables back to a linear scale to account for log-transformation bias (Sprugel 1983). Secchi and Ref_Secchi were then derived using the optimum models based on chl a and Ref_Chl a estimates. Finally, the current and reference TN and TP concentrations estimated by Abell et al. (2019) were combined with the chl a and Secchi estimates derived in this study to estimate TLI for 1031 lakes throughout NZ, using the following equations in Burns et al. (1999), where TLI is the average of the four components $TLn$, $TLp$, $TLc$ and $TLs$:(2) $TLn=-3.61+3.01\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}TN$(2) (3) $TLp=0.218+2.92\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}TP$(3) (4) $TLc=2.22+2.54\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}Chla$(4) (5) $TLs=5.10+2.27\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}\left(\frac{1}{Secchi}-\frac{1}{40}\right)$(5)

Current TLI was estimated by substituting estimated current values of water quality variables into Equations 2–5 whereas reference TLI was estimated using estimated reference values of water quality variables. Confidence intervals for TLI predictions were estimated by recalculating TLI 1000 times for each lake using the bootstrapped distributions of all individual TLI variables. 95% confidence intervals for lake-specific TLI estimates were then estimated based on the 2.5% and 97.5% quantiles of the bootstrapped values.

Results

There were marked differences between reference state and current values of chl a; e.g. mean chl a for current state (8.32 mg m⁻³) was ∼3.5 times higher than for reference state (2.42 mg m⁻³; Table 3). Predicted current chl a concentrations exceeded reference concentrations in 88% of lakes (see Supplementary Information for results for all lakes, including 95% confidence intervals). For 22% of lakes, the lower bound of the 95% confidence interval for the chl a prediction exceeded the upper bound of the 95% confidence interval for the Ref_Chl a prediction; i.e. chl a was higher than Ref_Chl a and the confidence intervals did not overlap in just over one-fifth of the lakes.

Table 3. Summary of estimated values of Trophic Level Index (TLI) variables for 1031 lakes.

Variable	Units	State	Percentiles					Mean
			10th	25th	50th	75th	90th
Chlorophyll a	mg m^−3	Reference	1.10	1.17	1.53	2.81	5.11	2.42
		Current	1.39	1.88	2.98	7.52	20.56	8.32
Secchi depth	m	Reference	6.01	8.01	10.35	11.44	11.71	9.62
		Current	2.12	3.59	5.88	9.37	10.90	6.48
TLI	TLI units	Reference	2.31	2.36	2.60	3.16	3.73	2.82
		Current	2.53	2.84	3.28	4.18	5.11	3.58

Mean Secchi under a reference state (9.62 m) was ∼1.5 times higher (indicating greater clarity) than the current state (6.48 m; Table 3). Predicted current Secchi was lower than reference Secchi in 89% of lakes. For 41% of lakes, the upper bound of the 95% confidence interval for the Secchi prediction was lower than the lower bound of the 95% confidence interval for the Ref_Secchi prediction; thus, current Secchi was higher and the confidence intervals did not overlap in approximately two-fifths of the lakes.

Predictions of current TLI exceeded reference TLI in 88% of lakes (Figure 3). For 27% of lakes, the lower bound of the 95% confidence interval for the current TLI prediction exceeded the upper bound of the 95% confidence interval for the reference TLI prediction. That is, in about one-quarter of the lakes, current TLI was greater than reference TLI with no overlap in the confidence intervals. Overall, the median and mean TLI values for a reference state both corresponded to an oligotrophic lake state (TLI = 2.0–3.0), whereas the median and mean TLI values for a current state both corresponded to a mesotrophic state (TLI = 3.0–4.0; Table 3).

Figure 3. Predicted TLI for 1031 New Zealand lakes corresponding to current and reference states. Vertical lines denote 95% confidence intervals. Lakes are ordered by ascending current TLI values.

The predominant lake trophic status under a reference state was oligotrophic (68% of lakes), while the most common lake trophic status under current conditions was mesotrophic (38% of lakes) followed by oligotrophic (32% of lakes). Approximately 5% of lakes corresponded to a eutrophic status under a reference state, while the current status of 16% of lakes was eutrophic. No lakes corresponded to supertrophic or hypertrophic conditions under a reference state, while the current status of approximately 15% of lakes corresponded to these trophic status classifications. Reference TLI values varied among lakes with different geomorphic types (Figure 4). Reference TLI values were generally lowest for glacial, landslide, volcanic and dam-formed lakes; median reference TLI values for each of these four geomorphic types corresponded to an oligotrophic status. Reference TLI values were generally highest for aeolian (dune), peat and shoreline lakes, with median reference TLI values corresponding to a mesotrophic status.

Figure 4. Estimated TLI for New Zealand lakes corresponding to reference and current conditions, differentiated by lake geomorphic type. Sample sizes are shown in parentheses. Only types with sample size ≥25 are presented and lakes with an undefined geomorphic type (n = 42) are not shown.

Change in TLI (current minus reference TLI) ranged from −0.84 to 3.31 (Figure 5). The mean change was an increase of 0.67 and the median change was 0.77, with a standard deviation of 0.67. TLI was estimated to have decreased in 11% of lakes, with glacial lakes comprising a large proportion of these lakes (Figure 5). TLI was estimated to have increased by 0–1 TLI units in 58% of lakes, by 1–2 TLI units in 26% of lakes and by >2 TLI units in the remaining 5% of lakes (Figure 5). The lakes with the greatest increases in TLI (>1 TLI unit) encompassed multiple geomorphic types, although riverine lakes were frequently among the lakes with the greatest departure from reference state (Figure 5). Further, the greatest increases in TLI tended to correspond to shallower lakes and there was a negative correlation between z_max and change in TLI (Spearman’s R = −0.52).

Figure 5. Frequency (counts) distribution of estimated departure from reference Trophic Level Index (TLI) for 1031 New Zealand lakes categorised based on geomorphic type.

Discussion

By building on work of other researchers (Snelder 2006; Leathwick et al. 2010; McDowell et al. 2013; Elliott et al. 2016; Abell et al. 2019), this study has derived estimates of chl a, Secchi and TLI for a national-scale sample of 1031 NZ lakes. These estimates are lake-specific and account for variability in factors that include geology, lake morphometry, hydrology, land cover and climate. This specificity means that we expect the reference state predictions to be more accurate than estimates derived from percentiles of measured data; e.g. as evaluated (alongside other methods) by Dodds et al. (2006) in relation to lakes and reservoirs in Kansas (USA) and by Huo et al. (2012) for lakes in southwest China. However, our results are likely less accurate than analysis of water samples collected from minimally disturbed lakes that are representative of the lake types considered; e.g. as undertaken to inform the National Lake Assessment to develop reference nutrient and chl a concentrations for different ecoregions in the USA (USEPA 2016), or as used to develop reference chl a concentrations for lake types in Europe based on analysing a dataset of 540 lakes deemed to be undisturbed (Carvalho et al. 2008). However, the pervasive extent of land use intensification and associated increases in catchment nutrient loads in NZ (Gluckman 2017; Snelder et al. 2018) means it is not possible to identify and sample adequate numbers of lakes that currently approximate a reference state in order to generate statistically reliable estimates that are representative of different typologies (see Schallenberg et al. (2018) for further review of the applicability of techniques to estimate lake reference conditions). Our methods are an extension of landscape-context statistical modelling approaches, which are recommended to derive reference condition estimates for lakes that lack reference state or paleolimnological measurements (Soranno et al. 2011). By incorporating databases and results from spatial modelling tools that are specific to NZ (Snelder 2006; Leathwick et al. 2010; McDowell et al. 2013; Elliott et al. 2016), and by striving to account for variability in lake processes, the TLI estimates derived in this study apply to a large and nationally representative sample of lakes, yet also reflect environmental variability and specific characteristics of lakes that are not routinely sampled.

The accuracy of the chl a, Secchi and TLI estimates will vary among the sample, reflecting that limited information was available for several factors. In particular, although a resuspension term was included in the Secchi model applied to shallow lakes, the models did not account for variability in coloured dissolved organic matter (e.g. expected to be high in peat lakes) or inorganic glacial flour that can influence water clarity (Vant and Davies-Colley 1984; Rose et al. 2014). This was justified as our analysis showed that developing separate Secchi – chl a relationships for these lake types did not improve model performance (based on AIC). However, as acknowledged above, evaluation of these groups was constrained by small sample sizes. Further, we acknowledge that geomorphic type is based on lake formation and limnological characteristics can vary substantially within individual groups; e.g. several ‘glacial’ lakes included in the sample are subject to limited influence by glacial flour and had average Secchi > 10 m (e.g. lakes Te Anau and Manapouri). With further analysis, catchment predictors such as wetland coverage might provide a suitable proxy for coloured dissolved organic matter in lakes or, alternatively, satellite remote sensing could be used to estimate optical characteristics of lake water (Lehmann et al. 2019). However, such characteristics are likely to be substantially different (but unknown) for several lake types under pristine conditions, limiting the value of these approaches to estimate reference state conditions. Alternatively, with greater statistical power (i.e. a larger sample size of monitored lakes), it may be possible to develop improved models that account for systematic variability among lake geomorphic types in light attenuation characteristics; e.g. to develop different models to predict Secchi in peat or glacial lakes. Overall, the TLI estimates are highly dependent on the estimated nutrient concentrations; e.g. the chl a component of TLI ($TLc$; Equation 4) is solely a function of the estimated nutrient concentrations, with chl a then used as a key predictor of Secchi. A simpler approach to estimating TLI is to take the average of $TLn$ (Equation 2) and $TLp$ (Equation 3); however, this would not allow lake specific chl a and Secchi values to be estimated and would also yield less accurate (albeit similar) TLI estimates (note that lake-specific values of three variables, in addition to chl a, are used to estimate Secchi in shallow lakes). A quantitative comparison of the two methods is presented in the Supplementary Information.

There is also uncertainty regarding the estimate that TLI has decreased in 11% of lakes, based on analysis of the median estimated TLI values for each lake. In some rare instances, improvements in trophic status (oligotrophication) may have occurred in NZ lakes as a consequence of hydrological manipulations; e.g. in reservoirs where historical dam construction has increased depth and water residence time of a pre-existing lake, thereby increasing nutrient attenuation due to settling and reducing the influence of sediment resuspension. However, our methods assumed that there was no change in hydrology or morphometry between reference and current states, hence the results are insensitive to such changes. Instead, the estimate that TLI decreased in 11% of lakes (Figure 5) may be biased high and potentially reflects a small number of lakes with estimated current nutrient concentrations that are anomalously low due to under-estimation of lake inflow nutrient concentrations using CLUES (especially TN; see Abell et al. (2019) for details of how nutrient concentrations were estimated). Uncertainty in TLI estimates is accounted for by the confidence intervals fitted using bootstrapping and it is notable that all confidence intervals on the left-hand side (low end of TLI range) of Figure 3 overlap, whereas this is not the case for the confidence intervals on the right-hand side (high end of TLI range) of Figure 3. This separation indicates that trophic status has changed little in lakes that currently have the lowest TLI values, whereas there has been a clear shift to higher TLI in the most eutrophic lakes in NZ. Uncertainty is higher (confidence intervals are wider) for the reference TLI estimates than the current TLI estimates (Figure 3) because the fit was poorer for the models used to predict reference state variables than models used to predict current state variables (Table 2).

Differences among geomorphic types in reference TLI estimates (Figure 4) are intuitive; generally, reference state TLI values were lowest in glacial lakes, which are typically deep and located in high elevation South Island catchments with low fertility soils. By contrast, reference state TLI values were highest in shoreline, dune, and peat lakes, which are typically shallower and present in lowland areas that receive higher subsidies of organic material from upstream catchments. Such trends are consistent with those described by Abell et al. (2019) in relation to nutrients (see that paper for further discussion).

The finding that change in TLI was generally greater for shallow lakes is consistent with recognition that many shallow lakes in NZ, especially in lowland areas, have been particularly influenced by human pressures (Schallenberg and Sorell 2009). This finding is also consistent with Özkundakci et al. (2014) who showed that indicators of ecological integrity were only weakly associated with human pressure indicators in deep NZ lakes. Our demonstration that the Secchi – chl a relationship varies between shallow and deep lakes (generally lower Secchi at a specific chl a value in shallow lakes) is intuitive given that resuspension of non-algal suspended sediment can make a substantial contribution to light attenuation in shallow lakes (e.g. Phlips et al. 1995). Variability among lakes in the Secchi – chl a relationship has long been recognised (Lorenzen 1980) and multiple studies have derived and contrasted Secchi – chl a relationships for different lake types (e.g. Canfield and Bachmann 1981; Mazumder and Havens 1998). However, previous studies that have used regression methods to predict Secchi in groups of lakes have generally not considered the contribution of non-algal particles; therefore, our approach of including a resuspension term in the model to predict Secchi in shallow lakes represents an improvement on the use of standard Secchi = f(chl a) regression models as it seeks to also account for the influence of sediment resuspension in a semi-mechanistic way. There is uncertainty about the threshold used to define ‘shallow’ lakes (z_max < 20 m); although this threshold was shown to minimise error in the estimates, it intuitively seems high as we are aware of lakes with z_max that approaches the threshold yet are not considered functionally shallow. It is possible therefore that this threshold is biased high to specifically account for sediment resuspension and its determination may have been confounded by the presence of multiple lakes in our sample with high TLI (low Secchi and high chl a) and moderate z_max (e.g. Lake Okaro, Bay of Plenty; z_max = 18 m) that are distinct from lakes with much lower TLI (high Secchi and low chl a) and much higher z_max, such as glacial lakes in the South Island with z_max > 100 m (e.g. Lake Wanaka, Otago). We recognise that mean depth ($\overline{\mathrm{z}}$) is superior to z_max for characterising lake typology based on depth; however, $\overline{\mathrm{z}}$ had not been estimated for most of the lakes in the WONI sample. For reference, the average ratio of z_max to $\overline{\mathrm{z}}$ for lakes with a defined $\overline{\mathrm{z}}$ in the NWQMN sample was 2.8 (n = 37, standard deviation = 1.7), indicating that, on average, a lake with z_max of 20 m has $\overline{\mathrm{z}}$ of ∼7 m. The definition of the fetch variable used in our Secchi model should be noted because we chose to use maximum lake fetch (i.e. maximum length) rather than proxies based on the square root of lake surface area (A^0.5) or half of the sum of the effective length (L) and effective width (W) (Davies-Colley 1988; L and W defined by Welch 1948). Compared with values of fetch derived from other proxies (A, L, W), our approach may differentially alter (reduce) fetch values of lakes with circular morphology compared with those that are long and narrow. Detailed studies to examine directionally dependent fetch would need to be integrated with comprehensive bathymetric mapping at an individual lake scale to improve definition of the fetch variable, and were considered to be beyond the scope of the current study. We expect that any model error associated with our definition of fetch is similar or smaller than error associated with estimates of mean wind speed for each lake.

Reservoirs (i.e. dam-formed lakes in Figures 4 and 5) were retained in the lake sample, although we recognise that a reference state is not directly applicable to these systems. We chose to retain reservoirs in the analysis as they include highly valued lakes that require management (Schallenberg et al. 2018). Reference condition estimates for a specific reservoir provide an indication of reference state that would apply to a natural lake that otherwise has similar morphometric and catchment characteristics. Thus, the finding that there are multiple reservoirs with current TLI greater than a reference value (by up to ∼2 TLI units; Figure 5) indicates that external nutrient load reductions in these systems could be undertaken to reduce TLI.

The estimates derived in this study provide benchmarks to inform lake and catchment management and restoration. These estimates may be complemented by broader metrics of freshwater ecological integrity (Özkundakci et al. 2014; Schallenberg et al. 2018; Kelly and Schallenberg 2019) to assess the state of lake ecosystems and develop restoration targets. Given the uncertainties described above, we recommend caution if using lake-specific reference values (see Supplementary Information) to inform restoration targets for individual lakes. Instead, it is preferable to consider ranges of values (including confidence intervals) for relevant lake classes (e.g. geomorphic types), as well as other lines of evidence that could include paleolimnological studies of similar lakes or traditional ecological knowledge. Additionally, other approaches, including deterministic lake ecosystem modelling (Schallenberg et al. 2017; Lehmann and Hamilton 2018) or extensions of other statistical modelling approaches could be used to predict reference water quality, particularly in unmonitored NZ lakes (Snelder et al. 2016), while remote sensing can be used to derive alternate estimates of current water quality in unmonitored lakes (Lehmann et al. 2019). Further, major changes in some catchments severely hinder the attainment of restoration targets that correspond to reference conditions prior to human disturbance; e.g. land management changes in agricultural areas such as the Waikato region have greatly enriched benthic sediment nutrient stores and considerably reduced depth and area of some lowland lakes (Hamilton et al. 2010). In such highly disturbed systems, conditions that correspond to a slightly altered state are likely to be the most ecologically ambitious targets that could be achieved (i.e. ‘Anthropocene baselines’; Kopf et al. 2015), with climate change a ‘shifting driver’ (Gillon et al. 2016) expected to modify the baseline ecological conditions that environmental managers can aim for in many lake ecosystems (Bennion et al. 2011; Me et al. 2018). A challenge for further research is therefore to develop restoration targets that are sufficiently protective, yet also technically and politically feasible in the face of global change.

Conclusion

This study quantified reference and current trophic state in a nationally representative sample of 1031 NZ lakes. It demonstrated a general trend of eutrophication in a broad sample of NZ lakes, reflecting pressures that include nutrient and sediment pollution associated with land use intensification. This study informs the status of lake ecosystems in NZ and provides valuable benchmarks to inform future management and assessment. These benchmarks will help to better quantify baseline ecological conditions and targets for restoration that will be achievable with global change.

Acknowledgements

We thank the New Zealand Ministry of the Environment, Regional Councils and their field staff who provided lake monitoring data via the Land, Air, Water Aotearoa partnership. Kevin Collier provided direction and insightful reviews that improved the manuscript. Britta Deichmann helped to make the map.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was funded by the New Zealand Ministry of Business, Innovation and Employment via theEnhancing the Health and Resilience of New Zealand Lakes programme [grant number UOWX1503].