Effect of climate variability on water quality trends in New Zealand rivers

Keywords: climate, El Niño–Southern Oscillation, New Zealand, rivers, trends, water quality.

Detection and characterisation of temporal trends are basic objectives of the water quality monitoring programs in operation around the world (Dawe 2006; Rosen and Lapham 2008; Kauffman et al. 2011; Ryberg et al. 2018). The analytical methods used to assess trends are well established, and quality-assured water quality datasets are widely available (e.g. Larned et al. 2016; Civan et al. 2018; Shoda et al. 2019; Helsel et al. 2020). Trend analysis serves multiple purposes, including evaluating changes in environmental state, assessing the effectiveness of management actions and policies, evaluating relationships between environmental conditions and the factors that influence them (i.e. driver or explanatory variables) and providing early warning of environmental problems (McMellor and Underwood 2014; Tomperi et al. 2016; Ministry for the Environment and Statistics New Zealand 2020; Murphy 2020).

In contrast to the widespread assessment of water quality trends, robust attribution of those trends to environmental drivers is less common. Attributing changes in water quality to environmental drivers is challenging due to uncertainty about causation and the interactive effects of multiple drivers, both anthropogenic and natural (Morse and Wollheim 2014; Dupas et al. 2016; Erb et al. 2016; Jackson et al. 2016; Ryberg et al. 2018). In the absence of empirical evidence about drivers, degrading trends in water quality are often casually attributed to anthropogenic drivers, such as agricultural intensification and urban development, based on anecdotal information (e.g. Johnson et al. 2009; Luo et al. 2011). Despite the challenges, arresting and reversing water quality degradation is predicated on the identification of contributing drivers, and predictive relationships between those drivers and water quality conditions.

It is reasonable to expect that water quality trends in most catchments are influenced by both natural and anthropogenic drivers, because both are ubiquitous at regional scales (e.g. natural climate variability, anthropogenic atmospheric nutrient deposition) and at catchment scales (e.g. natural plant succession, anthropogenic land use change). In several recent studies, temporal changes in coastal and fresh water quality were attributed to natural climate variability generated by climatic processes such as the El Niño–Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO; Mellander et al. 2018; Harding et al. 2019; Murphy and Sprague 2019). These processes produce regional-scale temporal variation in temperature and rainfall that influence catchment processes (e.g. contaminant mobilisation, transport, storage, transformation and dilution), which, in turn, affect water quality conditions (e.g. Gascuel-Odoux et al. 2010; Mellander et al. 2018).

There is continuing concern in New Zealand about water quality in freshwater and coastal environments (Suthanthangjai et al. 2013; Ministry for the Environment and Statistics New Zealand 2019). Water quality trend assessments of regional and national extent are conducted and reported as required by the New Zealand Resource Management Act and the Environmental Reporting Act (Ministry for the Environment and Statistics New Zealand 2017). Assessments of trends in water quality at 600–800 long-term river water quality sites across New Zealand have been conducted at 2- to 5-year intervals since the late 1990s (Larned et al. 2004, 2016; Ministry for the Environment and Statistics New Zealand 2015, 2017, 2019). The analysed trends covered a wide range of time windows and durations, and the reported patterns between successive analyses have been inconsistent. For example, two recent national State of Environment reports reported that 37 and 15% of monitored river sites in New Zealand had increasing nitrate trends for the 10-year windows ending in 2013 and 2017, and that 55 and 46% of the sites had increasing nitrate trends for the 20-year windows with the same end dates (Ministry for the Environment and Statistics New Zealand 2017, 2019). These inconsistencies have not been explained, but the results of an earlier study of water quality trends in New Zealand rivers suggests that climate variability is involved (Scarsbrook et al. 2003). Scarsbrook et al. (2003) found that trends in water quality variables for time windows of 5 and 10 years were reasonably consistent with trends in SOI over the same time windows. That study concluded that attributing changes in water quality to human management is difficult because of the confounding effect of climate variability. However, there has been no investigation of the degree to which climate variability can explain variation in water quality trends pertaining to different time windows and durations.

ENSO is a major source of interannual climatic variation in New Zealand. ENSO reflects a continuum of variation in sea level atmospheric pressure in the tropical Pacific Ocean, the extremes of which are designated El Niño and La Niña (Allan et al. 1996). Variation in ENSO strength is quasi-cyclic, with El Niño and La Niña events occurring on average every 2–7 years. ENSO strength is associated with regionally variable effects on climate conditions in New Zealand (Mullan 1996; Salinger and Mullan 1999). During El Niño events, New Zealand experiences stronger and more frequent westerly winds in summer, leading to dry conditions in eastern regions and increased rainfall in western regions. In winter, El Niño events are associated with more frequent southerly winds and colder than average temperatures. In spring and autumn, south-westerlies tend to be stronger and more frequent during El Niño events, producing higher than average rainfall in the south of the country and decreased rainfall in the east. La Niña events have weaker effects on New Zealand’s climate and are associated with more frequent north-easterly winds, higher than average rainfall in the north-east and decreased rainfall in the south and south-west. Variation in ENSO strength has been linked with regionally variable effects on hydrology in New Zealand (McKerchar et al. 1998; Mosley 2000; Chiew and McMahon 2002). Scarsbrook et al. (2003) showed that ENSO strength was associated with variation in river water quality observations at individual sites. Scarsbrook et al. (2003) also observed that the directions of national-scale trends (based on aggregating individual site trends) for some water quality variables corresponded to trends in ENSO strength.

In this study, we investigated relationships between trends in 10 river water quality variables at 56 New Zealand river monitoring sites and variation in ENSO strength over a 31-year period. The monitoring sites represented large main-stem rivers draining catchments with a wide range of environmental conditions, including near natural and highly modified land cover (Julian et al. 2017). We predicted that the strength and direction of water quality trends measured for specific time windows at both individual river monitoring sites and in aggregate would be associated with ENSO. We used statistical analyses to test these predictions and evaluate the proportions of variation in site and aggregated trends that could be explained by variation in ENSO. Our aim was to produce models that quantified the component of water quality trend strength and direction that can be attributed to ENSO variation for both site and aggregate trends.

The methods involved multiple steps aimed at providing four primary results, namely to quantify the proportion of variation in site and aggregate water quality trends explained by SOI trends, to explain variation in water quality trend responses to SOI between sites and to explain variation in aggregate water quality trend responses to SOI between variables, as shown schematically in Fig. 1.

Fig. 1.  Schematic overview of the methods. The input data are shown as trapezoids, processing steps are shown as rectangles and results are shown as ovals.

This study was based on data from the New Zealand National River Water Quality Network (NRWQN; Smith and Maasdam 1994; Davies-Colley et al. 2011). For the period from 1989 to 2019 (31 years), there was continuous monitoring of 56 sites located on 31 mostly large rivers (mean daily flow ranging from 0.8 to 511 m³ s^–1; median flow 96 m³ s^–1; Fig. 2). Monitoring records for the NRWQN sites comprise consistent monthly observations of 10 physicochemical variables (including discharge), which are referred to hereafter as ‘water quality variables’ (Table 1), with few censored or missing values. The exceptions were total nitrogen (TN) and ammoniacal nitrogen (NH₄-N), for which there are no observations during 1994 (Davies-Colley et al. 2011).

Fig. 2.  Locations of NRWQN monitoring sites and their catchments. The boundaries and names of New Zealand’s jurisdictional regions are indicated.

Table 1.  Water quality variables observed at the New Zealand National River Water Quality Network sites
Medians and ranges are for all 56 sites. NTU, nephelometric turbidity units

The values of water quality observations in rivers often vary in response to flow fluctuations due to two main mechanisms (Helsel et al. 2020). First, observations may decrease with flow due to dilution if the contaminant is delivered to the river at a constant rate. Second, observations may increase with flow if the contaminant is mobilised in the catchment in association with rainfall or in-stream in association with increasing turbulence. The responses of water quality observations to flow may be a combination of the two mechanisms. When trend analyses are conducted on river water quality data, the observations are often adjusted to account for the effect of flow (flow adjustment) for two reasons. First, flow adjustment can decrease extraneous variation (noise) in the variable of interest, so that trends (i.e. signals) are more confidently detected. Second, flow adjustment can remove the component of a water quality trend that can be attributed to a trend in flow (Helsel et al. 2020).

Flow adjustment is based on regressing the observations against their corresponding flows and then using the residuals of the regression model in the subsequent trend analysis. A variety of regression models may be used for flow adjustment (Helsel et al. 2020). We chose to use ordinary least-squares regression with log (base 10)-transformed values for both the water quality variable observations and flow observations. This form of regression is constrained to be monotonic (increasing or decreasing), which, in our experience, is generally a parsimonious and physically realistic model of the observation–flow relationship. In our experience, more flexible regression methods, such as locally weighted regression scatterplot smoothing, can produce unrealistic models.

Flow adjustment of the observations of each variable at each site, apart from instantaneous discharge at the time of sampling (FLOW), was performed by fitting a regression model to the entire time series. If this model was statistically significant (α = 0.05), the flow-adjusted observations were derived as the residuals of this model. Because the subsequent trend analysis was non-parametric, it was not necessary to back-transform the residuals. If the model was not significant, the variable was not flow adjusted (for further details, see Table S1 of the Supplementary material). At the end of this process, the observations, whether flow adjusted or not, were used in all subsequent analyses and are referred to as the water quality variables. To ensure consistency with studies in which flow adjustments are not used, we repeated all analyses using non-flow adjusted (‘raw’) data, and these results are provided in the Supplementary material.

We evaluated trends in each of the water quality variables at each monitoring site for rolling windows of 5-, 10- and 15-year duration starting in 1989 and incrementing by 1 year to a final window ending in 2019. These durations were chosen because they are typical of the timescales that are widely used in water quality trend assessments in New Zealand and elsewhere (Larned et al. 2016; Oelsner et al. 2017). This resulted in 27, 22 and 17 time windows of 5-, 10- and 15-year duration respectively. Trends were assessed for all site, variable and time window combinations for which data were available for 90% of months and 90% of years.

Trend strength and direction were evaluated using the Mann–Kendall trend test or, in the case that there was a statistically significant difference in the observations grouped by their sampling month (Kruskal–Wallis test, α ≤ 0.05), the seasonal Kendall test (Hirsch et al. 1982). The two procedures are robust to non-normal data distributions, outliers and non-linear trends and are widely used to evaluate the magnitude and significance of environmental trends. Both tests are based on the Mann–Kendall S statistic, which is the difference between the number of positive and negative changes in the water quality variable between all pairs of observations. We converted S to Kendall’s tau (τ) by dividing it by the total number of pairs of observations to provide a standardised measure of trend strength and direction. Kendall’s τ takes values between –1 and +1; a positive value indicates that the observations increased through time, whereas a negative value indicates that the observations decreased through time. Kendall’s τ for each site, water quality variable, duration and time window is denoted τ_w.

We quantified the aggregate trend for each water quality variable, duration and time window as the proportion of individual site trends with positive τ_w values; we refer to this statistic as P_w.

We used two independent variables based on the SOI to assess the degree to which fluctuations in trend strength and direction between time windows could be explained by variation in ENSO. The SOI is calculated as the normalised anomalies of the monthly mean sea level pressure difference between Tahiti and Darwin, Australia (Salinger and Mullan 1999). Monthly values of the SOI for the period 1989–2017 were obtained from the Australian Bureau of Meteorology (http://www.bom.gov.au/climate/current/soi2.shtml) and we used the Troup convention, whereby index values are multiplied by 10. The SOI typically ranges from –30 to 30 and is quasi-periodic with a typical period of 3–7 years. El Niño conditions are defined as an SOI less than –8 over ≥3 months and La Niña conditions are defined as an SOI >8.

The first independent variable described the trend in the SOI for a given time window. Although the SOI is quasi-cyclic, SOI values exhibit either an increasing or decreasing trend when a time window is defined by arbitrary starting and ending dates. We refer to the trend in the SOI for a specific time window (w) as the SOI trend, denoted by δSOI_w. We evaluated δSOI_w by linear regression of monthly SOI values against their respective dates for the same rolling windows of 5-, 10- and 15-year duration as the water quality trends.

The second independent variable described the extent to which there was an association between the cyclic fluctuations in the SOI and the corresponding monthly water quality observations at each site. This variable was the Pearson correlation between the monthly water quality observations and corresponding SOI values, which we calculated in four steps (Fig. 3). First, we selected all the observations of the water quality variable in the entire 31-year period of record at each site. Second, because regular seasonal variation in both the water quality measurements and the SOI values obscured their longer-term interannual variation, we deseasonalised both sets of values using classical seasonal decomposition (Hyndman and Athanasopoulos 2021). This analysis used centred annual moving averages to remove the variability associated with regular monthly changes for both the SOI and all water quality variables at each site. Third, because we wanted to assess only the extent to which cyclical fluctuations in the water quality observations and SOI were associated, we removed any secular temporal trend (i.e. any monotonic trend that was exhibited through the entire time series) for both the SOI and the water quality variables. Detrending was performed by regressing the deseasonalised values of the water quality variable and SOI against time using linear models and using the residuals of each model to represent the deseasonalised and detrended time series. Fourth, we calculated the Pearson correlation coefficient between the deseasonalised and detrended time series of the water quality variable and SOI. The correlation coefficients are referred to hereafter as r_o.

Fig. 3.  Example of the calculation of r_o. The plot shows TN concentrations at Hurunui River at SH1 in the Canterbury region (blue points and lines) and SOI values (red points and lines). Data are shown for the entire period of record (1989–2019). The points are the individual monthly observations of TN and the SOI, with lines indicating deseasonalised data (bold) and detrended data (fine). The dashed straight lines indicate the trends in the deseasonalised SOI and TN for the period of record. The resulting r_o value for TN at this site was 0.46.

We quantified the extent to which fluctuations in water quality trend strength and direction could be explained by the SOI for each site, water quality variable and duration in two steps (Fig. 1). First, because we wanted to assess only the extent to which cyclic fluctuations in water quality trends at each site were associated with corresponding fluctuations in the SOI trends, we removed any secular temporal trend that was exhibited in both sets of trends. For each water quality variable and duration, we used linear regression to remove the secular trend that was exhibited through the entire time series of τ_w values. We also used linear regression to remove any secular trend that was exhibited through the entire time series of the δSOI_w values. We regressed the τ_w and δSOI_w values against their corresponding end years and retained both sets of regression residuals (denoted τ_w′ and δSOI_w′ respectively). Second, we fitted the resulting data to models of the form:

where β₁ is a fitted coefficient that represents the response of the site water quality trend to the SOI trend. The intercept of Eqn 1 is 0, which represents our expectation that there will be no water quality trend response when the SOI trend for a time window is 0. The individual data points in the above model were not independent because they pertained to a rolling window through time. The models were therefore fitted using generalised least-squares regression, which provides robust estimates of the regression coefficients when there is temporal autocorrelation of the model residuals (Venables and Ripley 2013). In each model, we accounted for temporal autocorrelation by including a first-order autocorrelation structure. We used a pseudo-R² value (hereafter R²) to quantify the variation in τ_w′ explained by the model.

For each water quality variable and duration, we inspected the association between the variation in τ_w′ explained by the models (i.e. R² values) and the corresponding absolute values of r_o. We compared R² values for sites belonging to four categories defined by the absolute value of r_o, namely ≥0 (i.e. all sites), ≥0.1, ≥0.2 and ≥0.3. We expected that sites with larger absolute r_o values would have higher R² values because water quality observations at these sites have higher fidelity to fluctuations in the SOI.

We expected that the between-site variation in the response of water quality trends to the SOI trend for each variable, indicated by the fitted coefficients β₁, would be explained by the site r_o values. We modelled the between-site variation in responses of the water quality trends to the SOI trend by fitting 30 models (one for each water quality variable and duration) of the form:

where β₁ represents the site regression coefficients from Eqn 1 and β₂ is a fitted coefficient. The intercept of Eqn 2 is 0, which represents our expectation that there will be no water quality trend response to SOI trends for site and variable combinations with r_o of 0. We examined the R² values and overall significance (F-test P-values) of these models. Because we were examining multiple models (i.e. 10 variables × 3 time window durations), we controlled for false discovery by adjusting the P-values (Benjamini and Hochberg 1995).

We quantified the extent to which fluctuations in the aggregate water quality trends could be explained by the corresponding SOI trends (δSOI_w) in two steps. These steps are analogous to the steps described above to quantify the variation in the water quality trends at individual monitoring sites explained by SOI trends (Fig. 1). First, for each water quality variable and duration, we used linear regression to remove any secular temporal trend that was exhibited through the entire time series of both the aggregate water quality trends (i.e. the P_w values) and the corresponding SOI trends (i.e. the δSOI_w values). We regressed the P_w and δSOI_w values against their corresponding end years and retained both sets of regression residuals (denoted P_w′ and δSOI_w′ respectively). Second, for each water quality variable and duration, we fitted a linear regression model of the form:

where β₃ is a fitted coefficient that represents the response of the aggregate water quality trend to the SOI trend. The intercept of Eqn 3 is 0, which represents our expectation that there will be no P_w response when the SOI trend for a time window is 0. Because the individual data points in the above model were not independent, the models were fitted using generalised least-squares regression. Temporal autocorrelation was accounted for by including a first-order autocorrelation structure. We used a pseudo-R² value (hereafter R²) to quantify the variation in P_w′ explained by the models.

We expected that the between-variable variation in the responses of the aggregate trends to the SOI trend, indicated by the fitted coefficients β₃, would be explained by the median of the corresponding r_o values (denoted ). This is because, when is positive, the majority of sites’ responses to the SOI trend in each time window (δSOI_w) will be consistent with the sign of δSOI_w. Therefore, the aggregate trend (i.e. P_w) will increase and decrease when δSOI_w is positive and negative respectively. Conversely, when is 0, approximately equal numbers of sites will respond positively and negatively to the SOI trend in every time window. Therefore, the opposing influences of the SOI trend on variables with positive and negative values will tend to cancel out and will have little relationship with δSOI_w. We modelled the between-variable variation in responses of the aggregate quality trends to the SOI trend by fitting three models (one for each duration) of the form:

where β₃ represents the regression coefficients for each water quality variable from Eqn 3, and β₄ is a fitted coefficient. The intercept of Eqn 4 is 0, which represents our expectation that there will be no aggregate trend response to SOI trends for variables with of 0. We examined the R² values and overall significance (F-test P-values, controlled for false discovery) of these models.

Trend strength and direction at individual sites (i.e. τ_w) varied across all water quality variables and time windows (Fig. 4). The aggregate trend (i.e. P_w) oscillated with advancing time window for all variables and durations. For example, for the 5-year duration and turbidity (TURB), the median of values of τ_w were negative and, consequently, P_w values were <0.5 for four sets of time windows (windows ending 1993 and 1994, 1999–2002, 2005–07 and after 2016). The frequency of the oscillations in the aggregate trend directions decreased with increasing duration. For example, TURB had negative median values of τ_w and P_w values <0.5, for only two sets of time windows for the 10-year duration and one set for the 15-year duration (Fig. 4).

Fig. 4.  Distributions of Kendall tau values (τ_w) for all sites, for each water quality variable and assessment window for the 5-, 10- and 15-year durations and the proportion of positive trends (P_w). The central bar in each box indicates the median τ_w value, the box indicates the interquartile range and the whiskers extend to the 5th and 95th percentile values. The horizontal dotted line indicates a τ value of 0. Each panel corresponds to a water quality variable (rows) and a time window duration (columns).

The amplitude of the oscillation of the P_w values differed between variables. For example, P_w for dissolved reactive phosphorus (DRP) and temperature (TEMP) for the 10-year time window varied from ~0.2 to 0.9, but electrical conductivity (COND) and nitrate and nitrite nitrogen (NNN) varied from ~0.5 to 0.8 (Fig. 4). For some variables, there were obvious secular trends in the P_w values through the time series (i.e. a monotonic trend in the trends) that were a more dominant temporal pattern than the oscillations. This was particularly evident for the 15-year duration for NH₄-N (increasing P_w values with advancing time window) and total phosphorus (TP; decreasing P_w values with advancing time window).

The SOI trend (i.e. linear trend in the SOI; δSOI_w) differed between time windows, and alternated between increasing and decreasing, for all three durations (Fig. 5). The amplitude of the variation in the linear trend in the SOI decreased with increasing time window duration (Fig. 5).

Fig. 5.  SOI trends, defined by linear trends in the SOI, for time windows of 5, 10 and 15-year duration. Each point represents δSOI_w for the time window indicated by the end year (x-axis).

Correlations between the monthly deseasonalised and detrended water quality variable observations and monthly deseasonalised and detrended SOI values (r_o) varied in direction across sites and variables, ranging from –0.45 to 0.69 (Table 2). The median values of r_o were close to 0 for all variables, except TEMP and DRP, for which median values of r_o were >0.2. This indicates that TEMP and DRP were dominated by positive site r_o values, whereas there was a relatively even split between sites with positive and negative r_o values for the other variables. The proportion of sites with absolute values of r_o > 0.1 were lowest for TN and TP and were highest for TEMP, DRP and COND (Table 2).

Table 2.  Distribution characteristics of site r_o values for each water quality variable

Examples of the relationship between τ_w′ and δSOI_w′ for TN are shown in Fig. 6 for two sites for time windows of 10-year duration. The r_o values were negative for the Waipa River site (Fig. 6a) and positive for the Hurunui River site (Fig. 6b). The relationships in Fig. 6 were consistent with our expectation that fluctuations in water quality trends at individual sites are explained by the SOI trend (i.e. Eqn 1) and that variation in the water quality trend response between sites is explained by the respective r_o values (i.e. Eqn 2).

Fig. 6.  Examples of relationship between water quality trends (τ_w) and SOI trends (δSOI_w) for TN at two river monitoring sites for the 22 10-year duration time windows. (a) Waipa River at Otewa in the Waikato region, which has an r_o value of –0.36 and R² and β₁ values for Eqn 1 of 0.67 and –0.14 respectively. (b) Hurunui River at State Highway in the Canterbury region, which has an r_o value of 0.49 and R² and β₁ values for Eqn 1 of 0.72 and 0.10 respectively. Points are the individual water quality trends for each window with end years shown on the x-axis (blue) and the corresponding individual SOI trends (red).

The amount of variation in τ_w′ that was explained by δSOI_w′ (i.e. model represented by Eqn 1) varied for each combination of site, water quality variable and duration (Fig. 7). R² values tended to be greatest for the 10-year duration and least for the 15-year duration. When averaged across all variables, 13, 32 and 30% of sites had R² values >0.25 for the 5-, 10- and 15-year time window durations respectively. R² values tended to increase with increasing absolute values of r_o (Fig. 7). When averaged across all variables, 14, 36 and 30% of sites with absolute r_o values ≥ 0.1 had R² values that exceeded 0.25 for the 5-, 10- and 15-year time window durations respectively; this increased to 19, 38 and 30% respectively for absolute r_o values ≥0.2 and to 31, 44 and 37 respectively for absolute r_o values ≥0.3 (Fig. 7). These results are consistent with our expectation that the explanation of between-time window variation in trends (i.e. R² values) would be higher for sites with higher r_o values.

Fig. 7.  Cumulative frequency distribution of R² values measuring the strength of the association between τ_w and δSOI_w (Eqn 1). The distributions are shown for four categories defined by the absolute value of r_o (≥0; i.e. all sites, ≥0.1, ≥0.2, ≥0.3). Each panel corresponds to a water quality variable (rows) and a time window (columns).

After detrending (i.e. removal of secular trends), the r_o values explained between 1 and 79% of the variation in the between-site response of water quality trends to the SOI trends (i.e. β₁ coefficients fitted to Eqn 1; Fig. 8). This result was consistent with our expectation that the β₁ coefficients, indicating the direction of the response of the water quality trends to the SOI trend at each site, would be consistent with the site’s r_o values. The variation in the β₁ coefficients explained by r_o (R²) was consistently lower for the 15-year duration than for the 5- and 10-year durations (Fig. 8). The TURB, NH₄-N and TP models for all durations were not statistically significant (adjusted P > 0.05).

Fig. 8.  Fitted β₁ coefficients for Eqn 1 indicating the relationship between τ_w and δSOI_w for each site compared to r_o values (Eqn 2). Each panel corresponds to a water quality variable (rows) and a time window duration (columns). The red line indicates a regression fitted to these data and the annotations on each panel indicate the regression equation (i.e. Eqn 2), R² and adjusted P-values.

The magnitude of the variation in the aggregate water quality trends (P_w′) associated with δSOI_w′ (i.e. β₃ in Eqn 3) varied for each combination of water quality variable and duration (Fig. 9). The explained variation in P_w′ was greatest for DRP for the 5- and 10-year durations (R² values of 0.63 and 0.60 respectively) and TEMP and NNN for the 5-year duration (R² = 0.41 and 0.40 respectively).

Fig. 9.  Fitted β₃ coefficients for Eqn 3 indicating the relationship between P_w and δSOI_w for each variable compared to values (Eqn 4). Each panel corresponds to a duration. The red line indicates a regression fitted to these data and the annotations on each panel indicate the regression equation (i.e. Eqn 4), R² and adjusted P-values. Points are scaled according to each model’s R² values for the models represented by Eqn 3.

The values explained between 52 and 80% of the variation in the between-variable response of the detrended aggregate water quality trends to SOI trends (Eqn 4; Fig. 9). This result was consistent with our expectation that the response of P_w′ to the SOI trend would be consistent with the values. The variation in the β₃ coefficients explained by (R²) was lower for the 15-year duration than for the 5- and 10-year durations, but all models were statistically significant (adjusted P < 0.05; Fig. 9).

Five results of the analyses support our prediction that the strength and direction of water quality trends in New Zealand are associated with ENSO variation. First, for each water quality variable, fluctuations in individual observations of the variables were associated with fluctuations in ENSO indicated by the SOI (r_o; Table 2). Second, fluctuations in water quality trends between time windows (i.e. τ_w′) at individual sites were partly explained by the corresponding SOI trends (Fig. 6, 7). Third, between-site variation in the response of water quality trends to SOI trends was partly explained by r_o values (Fig. 8). Fourth, fluctuations in aggregate trends between time windows (i.e. P_w′) were partly explained by the corresponding SOI trends (Fig. 9). Fifth, variation in the response of aggregate trends for individual water quality variables was partly explained by the median of the site r_o values (; Fig. 9). These results were very similar when the analyses were performed using flow-adjusted and non-flow-adjusted water quality observations (see Fig. S1–S3).

Because the ENSO process is quasi-cyclic with irregular phases occurring on average every 2–7 years, there tend to be monotonic trends in the SOI at interannual time scales that we refer to here as ‘SOI trends’ (δSOI_w; Fig. 5). Partial explanation of the fluctuations in water quality trends between time windows at individual sites (i.e. τ_w) by the associated SOI trends (δSOI_w) indicates that, during periods when the SOI trend is positive, water quality trends tend to increase at sites with positive r_o values. The opposite applies to sites with negative r_o values. This pattern is reversed during periods when the SOI trend is negative.

Oscillations in aggregate water quality trends (P_w; Fig. 4) are likely to be caused by interannual fluctuations in rainfall and temperature over large spatial scales, which are driven by climatic processes such as ENSO and NAO over large spatial scales (Mullan 1996; Salinger and Mullan 1999; Mellander et al. 2018). In turn, the fluctuations in rainfall and temperature drive synchronous water quality responses (i.e. changes at the same time across sites). However, the strength and direction of those water quality responses depend on the specific response to the SOI at each site (i.e. r_o values). Our results also show that the strength of oscillations in aggregate trends, represented by the magnitude of the β₃ coefficient in Eqn 3, is influenced by the distribution of r_o values across the individual sites. Oscillations of greater magnitude tend to be associated with water quality variables with larger absolute values of (Fig. 9). These results indicate that some of the inconsistencies in the results of successive water quality trend analyses noted above (Ministry for the Environment and Statistics New Zealand 2017, 2019) are due to the effects of ENSO variation.

The methods set out in this study can be used to account for the effects of climate variability on individual and aggregate water quality trends. For an individual water quality variable and site, the combination of the r_o value and the SOI trend over a given time window provides an indication of whether climate variability has compensated for or amplified the effects of anthropogenic drivers. For example, for a 10-year duration time window and an SOI trend (δSOI_w) of 2.3 year^–1, the contribution of the SOI trend to a DRP trend at a site with an r_o value of 0.25 can be calculated in two steps. In Step 1, the coefficient β₁ for the site can be calculated from Eqn 2, and the β₂ coefficient shown in Fig. 8 is 0.045 (i.e. 0.18 × 0.25). In Step 2, the contribution of the SOI trend to the τ_w value for the site can be calculated from Eqn 1 as 0.1 (i.e. 0.045 × 2.3). This indicates that ENSO variation made a positive contribution to the assessed DRP trend that increased τ_w by 0.1 compared with a time window with a δSOI_w of 0 or compared with a site with an r_o value of 0. Similarly, the contribution of the same SOI trend to aggregate trends can be estimated using Eqn 4 and 3. For example, DRP in this study had an of 0.21 (Table 2) and the β₄ coefficient for a time window of 10-year duration was 0.41 (Fig. 9). Therefore, the coefficient β₃ is 0.086 (0.41 × 0.21) and the contribution of the SOI trend to the P_w is 0.2 (0.086 × 2.3). That is, for this 10-year window, the influence of ENSO increases the proportion of sites with positive trends in DRP by 20%.

Another way our results can be used to account for climate variability on water quality trends is to use r_o values to represent the response of trends for individual sites and variables to the SOI trend for a given time window. If measures of anthropogenic change in the catchments of sites are also available, statistical models can be used to elucidate the influence of those changes on τ_w by partialling out the variation attributable to the SOI trend (Snelder et al. 2021).

Attribution statements made based on our models could be made more robust by including uncertainty estimates based on the uncertainties of the underlying models. However, attribution of water quality trends to SOI trends should not be construed as comprehensively accounting for the effect of climate variability on water quality trends. Although ENSO strength has an important effect on New Zealand’s climate, it generally accounts for less than 25% of the year-to-year variation in seasonal rainfall and temperature (Salinger and Mullan 1999). Therefore, the SOI is an imprecise representation of the proximate climatic drivers of trends (i.e. rainfall and temperature) at any site. More research is needed to identify rainfall and temperature metrics that best explain variation in water quality trend strength and direction between time windows.

The variable r_o is a measure of the specific responses of a water quality variable at a site to fluctuations in ENSO. Understanding the mechanisms that produce the observed correlations may lead to methods that can more comprehensively account for the effect of climate variability on water quality trends. Previous studies have shown that the responses of hydrological metrics (e.g. mean, high and low flows) at river monitoring sites to fluctuations in the SOI are broadly consistent with the known regional effects of ENSO on rainfall (e.g. McKerchar et al. 1998; Mosley 2000). However, although the hydrological responses to variation in ENSO exhibit a degree of regional coherence, these can differ strongly between rivers that are geographically close (Mosley 2000). Similarly, Scarsbrook et al. (2003) showed there were regional patterns in water quality responses to the SOI, but that responses can differ strongly between sites that are geographically close. These findings indicate that flow responses to climate variability are mediated by hydrological processes such as water storage by catchments, which are also likely to mediate water quality responses. Scarsbrook et al. (2003) reported that the responses of some water quality variables to SOI variation were significantly correlated with the flow response to the SOI, but that this was not true for NNN, DRP, NH₄-N and coloured dissolved organic matter. This suggests that water quality responses to climate variability, for at least some variables, are mediated by catchment processes such as storage, mobilisation and biogeochemical transformation (Gascuel-Odoux et al. 2010; Green et al. 2014; Mellander et al. 2018). More research is needed to determine the causal links between climate variability and water quality responses. The benefits of this research would be more robust attribution of those trends to environmental drivers and an increased ability to assess the effectiveness of management actions and policies.

Our study shows that the effects of climate variability may amplify or counteract the effects of other drivers of water quality trends, even when those trends are assessed over time windows that are longer than ENSO cycles. Assessing and reporting water quality trends without robust attempts to identify the causes can lead to attribution of the trends to anthropogenic drivers alone, without supporting evidence (Ryberg et al. 2018). Speculative attribution could lead to the application of ineffective management actions to mitigate anthropogenic drivers and, equally, a failure to apply actions because a degrading trend is masked by the effects of climate variability.

In the present study, we used water quality trend directions and magnitudes assessed using Kendall’s τ, which is based on an underlying model of a monotonic trend (Helsel et al. 2020). This approach is a practical way to assess short-term changes in water quality or for analyses associated with hypothesised water quality changes. However, the underlying model of a monotonic trend obscures cyclic temporal patterns, and trend assessments of short duration may be dominated by the effects of climate. As time series of water quality data lengthen, studies that aim to describe patterns in water quality changes should incorporate the entire times series and use more flexible statistical models that detect and fit cyclic components, such as weighted regressions on time, discharge, and season (WRTDS; Hirsch et al. 2010, 2015). This type of more flexible modelling approach to trend analysis will provide a more complete picture of water quality changes as cyclic variation superimposed on longer-term secular trends.

The authors declare that they have no conflicts of interest.

This work was funded by the Ministry for the Environment and the Suitability program of the Our Land and Water National Science Challenge (Ministry of Business, Innovation and Employment contract C10X1507).