Attribution of river water-quality trends to agricultural land use and climate variability in New Zealand

Keywords: climate, land use, New Zealand, rivers, trends, water quality.

Deteriorating water quality and freshwater ecosystem health have been attributed to agricultural land use for over a century (Brooks 1916; Turner and Rabalais 2003). Efforts to arrest and reverse this deterioration often take the form of agricultural land use regulations, such as input and output controls, and mitigation and intervention requirements (Duncan 2014; Wiering et al. 2020). These regulatory actions are only effective if the observed deterioration is, in fact, due to agricultural land use and not to other anthropogenic or natural drivers (e.g. urbanisation, climate variability). Therefore, effective regulation of agricultural land use needs to be based on accurate attribution of water quality impacts to specific agricultural land uses and practices.

One of the stated aims of many freshwater quality networks is to identify relationships between water quality and land use (e.g. Davies-Colley et al. 2011; Behmel et al. 2016). However, the data from these networks are often used more narrowly to describe water quality state and detect temporal trends (e.g. Larned et al. 2016; Oelsner et al. 2017). Investigations of the environmental drivers that cause those trends are rare, as are investigations that rigorously attribute water quality trends to trends in drivers. By ‘rigorous attribution’, we mean quantitative analyses of relationships between water quality trends and putative drivers, and consideration of multiple alternative drivers (Ryberg 2017; Ryberg et al. 2018; Murphy 2020). Weaker alternatives to rigorous attribution include qualitative reasoning, references to previous studies and simple speculation (e.g. Luo et al. 2011; Joy et al. 2019). When attribution is weak, there is a risk that efforts to improve water quality by managing the presumed driver will be ineffective. Similar risks have been identified for cases of weak attribution of trends in river flows and floods (Merz et al. 2012; Harrigan et al. 2014).

Despite the need, examples of rigorous attribution of water quality trends to drivers are rare for two reasons. First, suitable data to characterise spatiotemporal variation in environmental drivers are scarce and fragmented (Merz et al. 2012; Diamantini et al. 2018; Mellander et al. 2018; Ryberg et al. 2018). Suitable data consist of time series of measurements of driver magnitude, with durations and frequencies that correspond to water quality time series and are spatially congruent with water quality monitoring sites (e.g. Howden and Burt 2009; Bouraoui and Grizzetti 2011). Second, water quality is generally influenced by multiple environmental drivers, including anthropogenic drivers, such as land use, and natural drivers, such as climate variability. There may be additive, compensatory or synergistic interactions among these drivers, making it difficult to reliably attribute water quality responses to specific land use or management actions (Morse and Wollheim 2014; Dupas et al. 2016; Ryberg et al. 2018; Choquette et al. 2019; Murphy 2020).

The situation has progressed recently, with rapid growth in the generation and provision of data for characterising drivers at large spatial scales and improved analytical methods. In the past decade, attribution studies have been performed using statistical analyses of relationships between water quality trends and temporal changes in drivers (Diamantini et al. 2018; Mellander et al. 2018), structural equation models (Ryberg 2017; Ryberg et al. 2018), joint hydrology–water quality empirical models (Choquette et al. 2019; Murphy and Sprague 2019) and deterministic models (Gascuel-Odoux et al. 2010).

Here we report the results of a large-scale statistical analysis of relationships between agricultural land use drivers and water quality trends in New Zealand rivers. Agriculture is the dominant land use in New Zealand, accounting for 49% of the land area, of which over 80% is used for pastoral farming; the remaining agricultural land is used for plantation forestry, cropping and horticulture. The negative effects of pastoral land use on water quality in freshwater and coastal environments are longstanding public issues (Suthanthangjai et al. 2013; Ministry for the Environment and Statistics New Zealand 2019). Repeated studies have shown that water quality state (i.e. the characteristic current condition) in New Zealand rivers is strongly explained by the proportion of pastoral agricultural land cover in the upstream catchment (Larned et al. 2004, 2016; Davies-Colley 2013), but these studies did not include attribution of water quality trends (i.e. temporal changes in condition) to land use or other drivers.

After over a century of predominately low-intensity pastoral agriculture in New Zealand, a period of agricultural intensification and diversification commenced in the early 1980s, partly in response to economic deregulation and the removal of farm subsidies (Smith and Montgomery 2004; MacLeod and Moller 2006). The changes included increased fertiliser and supplementary feed input, expansion of irrigation schemes, growth of the dairy, deer farming and plantation forestry sectors and contraction of the sheep and beef sectors (MacLeod and Moller 2006). In response to public concerns about water quality impacts, regulatory reforms have been enacted by the New Zealand government over the past decade that aim to prevent and reverse river water quality degradation associated with agriculture (New Zealand Government 2014, 2017, 2020a). These reforms increased requirements for the agriculture sector to reduce impacts on aquatic receiving environments through measures such as improved management of fertiliser, livestock effluent and irrigation water, reducing stock access to streams, riparian protection and tree planting on erodible hill country (Monaghan et al. 2021). In addition, there have been industry-led initiatives such as the Dairying and Clean Streams accord, which requires the exclusion of dairy cattle from channels and riparian zones of perennial rivers greater than 1 m wide (Bewsell et al. 2007).

River water quality monitoring programs are well established in New Zealand and time series of observations have increased over time; there are now more than 1000 sites at which a range of water quality variables have been observed on a monthly or quarterly basis for up to 30 years (Larned et al. 2016). These river water quality data have been used to produce numerous reports on water quality state and trends, and a smaller number of quantitative comparisons of water quality trends and environmental drivers. In the latter cases, a single environmental driver was considered in some studies (e.g. Hamill and McBride 2003; Scarsbrook et al. 2003), and multiple drivers were considered in others, but only as static descriptions of land uses (e.g. Ballantine and Davies-Colley 2014; Julian et al. 2017). We know of no previous investigations that met the requirements set out above for rigorous attribution (i.e. quantitative analyses that consider multiple alternative drivers and contemporaneous variation in those drivers).

Climate variability is a complication in attributing changes in water quality to anthropogenic drivers (Chanat and Yang 2018; Ryberg et al. 2018; Choquette et al. 2019; Murphy and Sprague 2019). At catchment scales, water quality responses to all environmental drivers are mediated by biophysical processes, including contaminant mobilisation, transport, attenuation and dilution (Mellander et al. 2018). These processes are strongly influenced by interannual climate variability, producing variation in water quality that may mask responses to anthropogenic drivers (Choquette et al. 2019). The primary sources of interannual climate variability are global-scale cyclic processes such as the El Niño–Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO); these processes produce temporal variations in temperature and rainfall, which, in turn, influence the biophysical processes identified above (e.g. Gascuel-Odoux et al. 2010; Harding et al. 2019; Mellander et al. 2018).

In New Zealand, ENSO is the dominant mode of interannual variation in rainfall and air temperature (Salinger and Mullan 1999; Ummenhofer et al. 2009). Consequently, ENSO also influences interannual variation in river flow regimes and water quality, and these ENSO effects may obscure water quality trends caused by anthropogenic drivers such as land use change (Mosley 2000; Chiew and McMahon 2002; Scarsbrook et al. 2003; Zeldis et al. 2008).

In this study, an extensive spatiotemporal dataset describing agricultural land use and land cover was combined with data describing variation in ENSO strength and observations of six water quality variables made between 1990 and 2017 at 937 river monitoring sites across New Zealand. We used these data to investigate potential drivers of trends over discrete time windows of 10 and 20 years within the 28-year study period.

The methods involved multiple steps aimed at providing two primary results, which are 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.

River water quality trends were assessed for six variables, namely dissolved reactive phosphorus (DRP; mg m^–3), ammoniacal nitrogen (NH₄-N; mg m^–3), nitrate nitrogen (NO₃-N; mg m^–3), total nitrogen (TN; mg m^–3), total phosphorus (TP; mg m^–3) and visual clarity (CLAR; m). Water quality data for sites distributed across New Zealand were acquired from the 16 regional councils and the National Institute of Water and Atmosphere (NIWA). Each council and NIWA operate a network of river monitoring sites at which water quality sampling is performed monthly or quarterly. We acquired all available data from January 1990 to December 2017. These dates cover a period from when regular water quality monitoring in most regions commenced to the most recent acquisition and federation of all data from the 16 regional councils. As part of federating the datasets, we undertook data checks and corrections described by Larned et al. (2016). For most variables, there was variation in analytical methods across the data collecting agencies. For each variable, only data corresponding to the most widely used and comparable procedures were retained and the remaining data were omitted, as described by Larned et al. (2016). The individual datasets contained censored values, indicating that reported values were below the analytical detection limit or above the reporting limit. Censored values were identified in all datasets, and these entries were consistently indicated by the combination of the reported values and a flag indicating the type of censoring. The individual datasets were then combined into a single consistent dataset comprising 918 389 observed values of the variables at 1051 sites (Fig. 2).

Fig. 2.  Map showing the location of all sites used in this study and the boundaries of jurisdictional regions.

We quantified variation in the ENSO cycle with the SOI, which is calculated as the normalised anomalies of the monthly mean sea level pressure differences between Tahiti, French Polynesia and Darwin, Australia (Salinger and Mullan 1999). The SOI ranges from –3 to 3 and is quasi-periodic with a typical period of 3–7 years. An El Niño phase is defined as SOI <0 and a La Niña phase is defined as SOI >0. Monthly values of the SOI for the period from 1989 to 2017 were obtained from the Australian Bureau of Meteorology (http://www.bom.gov.au/climate/current/soi2.shtml) and were handled using the Troup convention, whereby index values are multiplied by 10.

We used the Pearson correlation coefficient to describe the extent to which variation in each water quality variable at each river monitoring site was associated with fluctuations in the SOI (Fig. 1, Step 1). We calculated the value of this correlation coefficient for each site and variable in four steps. First, we selected all the observations of the water quality variable in the entire 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 seasonal changes for both the SOI and all water quality variables at each site. Third, because we wanted to assess only the extent to which fluctuations in the water quality variable and SOI were associated, we removed 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.

Although the SOI is cyclic, there is a tendency for SOI values to exhibit either an increasing or decreasing trend when a time window is defined by arbitrary starting and ending dates. We quantified the trend in the SOI for each time window (w) as δSOI_w (Fig. 1, Step 2). We evaluated δSOI_w by linear regression of monthly SOI values against their respective dates for rolling windows of 10- and 20-year duration starting in 1990 and incrementing by 1 year to a final window ending in 2017. This resulted in 19 and 9 time windows of durations of 10 and 20 years respectively.

For the same rolling windows of 10- and 20-year duration defined above, we characterised the direction and strength of the monotonic trend for each water quality variable at each monitoring site (Fig. 1, Step 3). Trends were assessed using Kendall’s tau (τ) statistic, which was calculated using the Mann–Kendall trend test or the seasonal Kendall test (Hirsch et al. 1982) for cases where there was a statistically significant difference in the observations grouped by their sampling month (Kruskal–Wallis test α ≤ 0.05). The two procedures are robust to non-normal data distributions, missing values, outliers and non-linear trends, and are widely used to evaluate the magnitude and significance of monotonic trends in environmental variables. Both tests are based on the Mann–Kendall S statistic, which is the difference between the number of positive and negative changes in the variable between all pairs of observations. When pairs of observations included censored values, the changes were evaluated robustly using the methods described by Helsel (2011). We converted S to Kendall’s τ by dividing it by the total number of pairs of observations to provide a standardised measure of trend direction and strength. Kendall’s τ takes values between –1 and +1, with negative and positive values indicating that the observations decreased and increased through time respectively. Kendall’s τ for each site, water quality variable, duration and time window is denoted as τ_w.

We restricted the allowable proportion of missing values to ensure that the observations consistently and adequately represented both the seasons and the entire time window. Site and variable combinations were retained for trend analysis provided 80% of sample intervals were represented by observations.

The values of water quality variables in rivers often vary in response to flow fluctuations (Hirsch et al. 1991). When trend analyses are conducted on river water quality data, the observations are often adjusted to account for the effects of flow (flow adjustment) to decrease extraneous variation and increase statistical power (i.e. to increase the confidence of trend assessments; Hirsch et al. 1982). We did not flow-adjust the water quality data because our analyses concerned relationships between trend direction and strength across multiple sites, and were not dependent on establishing the highest possible level of confidence in the direction of individual trends.

Strong relationships are consistently observed in New Zealand between water quality state and the proportion of catchment occupied by grazed grassland (i.e. used for pastoral farming; e.g. Larned et al. 2016; Julian et al. 2017). The area of grazed grassland has not changed substantially in most catchments over the past 30 years (Ministry for the Environment and Statistics New Zealand 2018). However, over the same period there have been significant changes in the spatial distribution and intensity of the four main types of pastoral livestock, namely sheep, dairy cows, beef cows and deer (MacLeod and Moller 2006; Ministry for the Environment and Statistics New Zealand 2018). Where changes in catchment area occupied by grazed grassland have occurred, they are generally associated with commensurate changes in area of plantation forests, with grazed grassland in some areas being replaced by forest and vice versa (New Zealand Government 2020b). In addition, within each of the four main types of pastoral farming (sheep, dairy, beef and deer farming) and plantation forestry, there have been differences in land use practices, such as the types and rate of adoption of mitigation measures (Monaghan et al. 2021).

For this study we used data describing the spatiotemporal variation in the four main types of pastoral livestock and in the two relevant land cover categories, namely grazed grassland and plantation forest. However, we did not have data describing changes in specific land use practices. Therefore, we derived two sets of catchment land use indicators for each monitoring site and each time window describing: (1) the change in land use intensity (δLU; Fig. 1, Step 4); and (2) the mean land use intensity in the catchment (μLU; Fig. 1, Step 4). The first set of indicators are measures of changes in the intensity of the four types of pastoral livestock and plantation forestry. The five change indicators are direct measures of change in land use intensity that have potentially driven water quality trends. The second set of indicators are measures of the mean intensity of the four types of pastoral livestock and plantation forestry. The mean indicators are indirect measures of potential changes in land use practices that are associated each of the five types of land use; they provide proxy measures of a catchment’s exposure to any land use practices that are particular to each type of land use.

The area occupied by different land cover categories is mapped nationally in New Zealand using the Land Use and Carbon Analysis System (LUCAS; Newsome et al. 2018), which includes 12 land cover categories derived from satellite imagery corresponding to four dates: 1990, 2008, 2012 and 2016. Land cover categories for the 1990 map were determined from 30-m spatial resolution Landsat-4 and Landsat-5 satellite imagery and the 2008, 2012 and 2016 maps were prepared using higher-resolution imagery from the SPOT-5 satellite (Ministry for the Environment 2012). Land cover classes were mapped at a minimum mapping area of 1 ha and include the categories grazed grassland, plantation forest planted before 1990 and plantation forest planted after 1989.

The number of animals in the four stock type categories is periodically surveyed on all livestock farms by Statistics New Zealand as part of the Agricultural Production Census (APC). We used the highest resolution versions of APC data that are publicly available, which are associated with a spatial coverage comprising 960 hexagonal grid cells (35 000 ha) that cover all New Zealand (https://statisticsnz.shinyapps.io/livestock_numbers/). For each grid cell, the number of animals in each stock type category was available for 5 census years (1994, 2002, 2007, 2012, 2017).

We combined the land cover and stock numbers data with a digital representation of New Zealand’s surface water drainage network to provide a spatially continuous representation of catchment agricultural land use for the entire country. The digital network represents New Zealand’s rivers as 590 000 segments (delineated by upstream and downstream confluences) and their catchments and is contained in a geographic information system (Snelder and Biggs 2002). We converted the stock numbers associated with the hexagonal grid and each census year into the density of animals of each stock type in the catchment of every network segment in four steps. First, all land cover categorised by LUCAS as grazed grassland (combining high- and low-producing subcategories) was intersected with the hexagonal grid. Second, the animals of each stock type in each hexagonal grid cell were evenly distributed over the land categorised as grazed grassland within the cell and the animal numbers were then converted to densities (number of animals per hectare of grazed grassland). The animal numbers before 2008 were distributed to the high- and low-producing grasslands associated with the 2008 version of LUCAS and the animal numbers in 2012 and 2017 were distributed to the high- and low-producing grasslands associated with the 2012 and 2016 versions of LUCAS respectively. Third, the combined hexagonal grid cells and grazed grassland spatial coverages were intersected with the digital network catchment boundaries. Finally, the number of animals of each stock type in each catchment was calculated as the sum of each area of grazed grassland within each catchment multiplied by the corresponding animal densities.

The animal densities were converted to stock unit (SU) equivalents to provide a measure of land use intensity that is comparable across the different stock types (Di and Cameron 2002; McDowell and Wilcock 2008; Ledgard et al. 2011). The SU equivalent is a commonly used measure of metabolic demand by livestock in New Zealand (Parker 1998; Trafford and Trafford 2011). The density of animals of each type in each catchment was converted to stocking intensity (SU ha^–1) by multiplying by their SU equivalent. The conversion accounted for the increase in metabolic demand of animals of each stock type over time due to increases in animal weight and levels of production. For example, the production of milk solids per dairy cow has risen in New Zealand from ~260 kg year^–1 in 1990 to 370 kg year^–1 in 2017 (Livestock Improvement Corporation Limited and DairyNZ Limited 2018). Our calculation of land use intensity accounted for this increasing animal size and productivity by adjusting SU values over time according to the values presented in Table 1.

Table 1.  Stock unit equivalent values assumed for sheep, beef, dairy and deer livestock classes between 1980 and 2017

Because changes in areas of grazed grassland since 1990 were generally associated with changes (either increases or decreases) in plantation forest, plantation forest land cover was included as a land use indicator. We calculated the proportion of catchment area occupied by plantation forest in the catchment of every network segment in two steps. First, for each version of LUCAS (1990, 2008, 2012 and 2016), we intersected the areas categorised as plantation forest before 1990 and after 1989 with the digital network catchment boundaries. Second, the area of plantation forest planted before 1990 and after 1989 in each year and each catchment was expressed as a proportion of the total catchment area (%).

The indicators used to characterise changes in land use intensity were evaluated for each time window and every segment of the network as the rates of change of stocking intensity and rate of change in the proportion of catchment occupied by plantation forest. Change in stocking intensity was evaluated as the change in total SUs normalised by the catchment area, divided by the time step in years (δTotal; SU ha^–1 year^–1). We also evaluated the change in the proportion of the total SUs associated with each stock type in the catchment for the time step (δDairy, δBeef, δSheep and δDeer; % year^–1). Change in plantation forest intensity was evaluated as the change in the proportion of catchment occupied by plantation forest (δForest; % year^–1).

The mean intensity indicators were evaluated for each time window and every segment of the network as the mean stocking intensity and the mean proportion of catchment occupied by plantation forest. Mean intensity was computed for total stocking intensity (μTotal; SU ha^–1), for the proportion of total stocking intensity associated with each stock type (μDairy, μSheep, μBeef and μDeer; %) and for the proportion of catchment occupied by plantation forest (μForest; %). We used linear interpolation to estimate the change in intensity and mean intensity when the beginning or end of time windows did not coincide with a year for which data were available. For the beginning of the early 20-year assessment window (1990) we used stock numbers pertaining to 1994. For the end of the late 20- and 10-year assessment windows (2017) we used plantation forest areas pertaining to 2016. Geographic information for each monitoring site was used to associate it with a segment of the digital network, and the land use indicators for that segment were used as explanatory variables in the subsequent analyses.

For each water quality variable, time window and duration, we used regression models to explain variation in the site trends (i.e. τ_w) as a function of each site’s correlation with the SOI (r_o) and both types of land use indicator (i.e. the change in land use and the mean land use; Fig. 1, Step 5). All evaluated τ_w values were used in all analyses regardless of their significance. The land use indicators were intercorrelated because they are mutually dependent (i.e. decreases in one indicator are at least partly compensated for by increases in another). To allow the sign of the fitted model coefficients to be interpreted as the direction of the relationship between the explanatory variables and the trends, we reduced the land use indicators for all combinations of sites and time windows to a set of orthogonal explanatory variables using principal component analysis (PCA). The values of each land use indicator for each site and time window were used as input to a PCA analysis, which was based on the correlation matrix so that each indicator was given the same weight. We performed a varimax rotation on the PCA, which assists interpretability by increasing the correlation of some of the individual input variables with the resultant components. We retained the components that made a significant contribution to the explained variation in the input data (broken stick method; Jackson 1993). We used the coordinates of the sites on the retained PCA axes as explanatory variables and refer to these as synthetic land use indicators denoted by PC1_w, P21_w,… PCn_w. We interpreted the synthetic land use indicators based on the loadings of the original land use indicators on the retained PCA components.

For each water quality variable, time window and duration, we fitted ordinary least-squares linear regression models of the form:

where τ_w represents the Kendall’s τ for each site, time window and duration, and β₁, β₂,…β_n are coefficients fitted to the explanatory variables that include r_o and PC1_w, PC2_w,…PCn_w.

We inspected whether there were geographic patterns in all model residuals, the occurrence of which would suggest environmental drivers of trends that were not represented by the models. We formally assessed the strength of the geographic patterns, using the Mantel test (Mantel 1967). The Mantel statistic (Mantel’s r) is the Pearson correlation coefficient between two dissimilarity matrices and is used to quantify geographic clustering (i.e. whether there is greater similarity between sites that are geographically close than between widely separated sites). The first matrix described the dissimilarity in the model residuals between pairs of monitoring sites. The second dissimilarity matrix defined the geographic (Euclidian) distance between all pairs of sites. The significance of Mantel’s r was established by permutation based on the null hypothesis of no correlation (Legendre and Legendre 1998). We controlled for false discovery by adjusting the P-values (Benjamini and Hochberg 1995).

We assessed the performance of the models and their overall significance based on their coefficients of determination (R²) and the F-test P-values (Fig. 1, Step 6). Because we were examining multiple models (i.e. 6 variables × 9 windows of 20-year duration and 19 windows of 10-year duration), we controlled for false discovery by adjusting the P-values (Benjamini and Hochberg 1995). We assessed the direction of each explanatory variable’s relationship with τ_w (i.e. positive or negative) based on the sign of the fitted coefficient and quantified confidence in that direction based on the complement of half the fitted coefficient’s adjusted P-value (i.e. 1 – P/2; Makowski et al. 2019). Confidence in the direction of the relationships between the explanatory variables and τ_w ranged between 0.5 (as likely to be the opposite of the direction indicated by the fitted coefficient) to 1 (complete confidence). For each water quality variable and duration (10 and 20 years), we summarised confidence in the direction of the relationship between each explanatory variable and τ_w as the mean of the individual confidence values pertaining to that explanatory variable over all models (i.e. windows).

For each water quality variable and each duration, we defined the consistency of each explanatory variable’s relationship with τ_w as the degree of agreement of the directions indicated by the fitted β coefficients over all time windows (Fig. 1, Steps 7 and 8). Our definitions of agreement differed for the coefficient fitted to r_o (i.e. β₁) compared with the synthetic land use indicators. For r_o, we expected that for time windows when the SOI trend (δSOI_w) was positive, the direction of the relationship indicated by β₁ would be positive and for windows when δSOI_w was negative the direction would be negative. We measured the consistency of the contribution of r_o for each water quality variable and duration as the fraction of time windows for which this expected relationship was true. For the synthetic land use indicators, we had no prior expectation of the direction of the relationship with τ_w. Therefore, for each indictor, we measured consistency as the fraction of models for which the fitted coefficients indicated the same direction as the mode direction:

where N_p is the number of models for which the fitted β coefficients were positive and N is the total number of models (i.e. 19 and 9 for the 10- and 20-year duration time windows respectively). The measure of consistency for each of the synthetic land use indicators has a minimum value of 0.5 when there are as many positive β coefficients as negative and a maximum value of 1 when all β coefficients have the same sign. We arbitrarily classified the relationships for each explanatory variable and duration as consistent and certain if consistency was ≥0.9 and mean confidence was ≥0.8.

The SOI trend (i.e. linear trend in the SOI; δSOI_w) differed between time windows and alternated between increasing and decreasing for both time window durations (Fig. 3). The amplitude of the variation in the linear trend in the SOI was greater for the 10- than 20-year time window duration (Fig. 3).

Fig. 3.  SOI trends, defined by linear trends in the SOI, for time window durations of 10 and 20 years. Each point represents the SOI trend (δSOI_w) for the time window indicated by the end year (x-axis).

Correlations between the monthly deseasonalised and detrended water quality observations and monthly deseasonalised and detrended SOI (r_o) varied considerably in strength and direction across sites and variables, ranging from –0.95 to 0.88 (Table 2). The means of the absolute values of r_o were lowest for TP (0.23) and highest for TN (0.28).

Table 2.  Distributions of r_o values by water quality variable

After discarding time windows for which <80% of sample intervals and <80% of years were represented by observations sites and variables, the lowest number of sites for any window was 78 for CLAR for the 10-year window ending 1999. The number of sites analysed increased for all water quality variables over the study period (Fig. 4). For example, for DRP, the number of sites increased from 122 for the 10-year window ending 1999 to 821 for the window ending 2017. Information describing the proportions of observations that were censored and missing is provided in the Supplementary material (Fig. S1).

Fig. 4.  Summary of trend analyses. The left-hand plot indicates the number of site trends analysed for each time window and duration. The central and right-hand boxplots indicate the distribution of Kendall’s τ (τ_w) for all sites, for each water-quality variable and assessment window for the 10- (red) and 20-year (blue) durations. The central bar indicates the median 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.

Trend strength and direction (i.e. τ_w) varied across sites for all water quality variables, time windows and durations (Fig. 4). Trend directions generally oscillated with advancing time window for the analyses of 10-year duration. For example, for CLAR and the 10-year duration, the median values of τ_w were positive for windows ending 2000–05, negative for windows ending 2006–11 and then positive for windows ending 2012–17. The oscillations in trend directions were weaker for the analyses of 20-year duration. For example, there were few changes in the sign of the median values of τ_w for the 20-year duration (Fig. 4). For the 20-year duration, the median values of τ_w were negative (i.e. most sites had decreasing trends) and decreased monotonically (i.e. the number of sites with decreasing trends increased) with advancing time window for DRP and TP, whereas median values of τ_w were positive (i.e. most sites had increasing trends) for all time windows for NO₃-N.

Spatial variation in the land use indicators across the entire study time period (i.e. 1992–2017) is shown in Fig. 5 and 6. Over the 28-year study period, the mean total stocking intensity (μTotal) across the entire river network was 3.8 SU ha^–1 (interquartile range, IQR, 0–6.8 SU ha^–1) and across all river monitoring sites it was 5.3 SU ha^–1 (IQR 1.6–7.7 SU ha^–1). Over the 28-year study period, the mean change in total stocking intensity (δTotal) across the entire river network was 0.25 SU ha^–1 (IQR –0.29, 0.58 SU ha^–1) and across all monitoring sites it was 0.6 SU ha^–1 (IQR –0.46, 1.01 SU ha^–1). μTotal was greatest in Waikato, Taranaki, Manawatū and Hawkes Bay regions of the North Island and the Canterbury and Southland regions of the South Island (Fig. 5). δTotal increased the most in the Waikato and Taranaki regions of the North Island and the Canterbury and Southland regions of the South Island (Fig. 5). δTotal decreased the most in the Auckland, Manawatū–Whanganui, Hawkes Bay and Wellington regions of the North Island and the Marlborough region of the South Island (Fig. 5).

Fig. 5.  River network segments (Strahler order ≥ 4) categorised by land use indicators over the entire study period, evaluated as (A) mean total stocking intensity (μTotal; SU ha^–1), and (B) changes in stocking intensity (δTotal; 100 × SU ha^–1 year^–1); (C) mean proportion of catchment cover by plantation forest (μForest; %), and (D) changes in proportion of catchment cover by plantation forest (δForest; % year^–1 × 100). Categories are defined by six equal quantiles of the distributions of the non-zero values of the indicators over network segments.

Fig. 6.  River network segments (Strahler order ≥5) categorised by land use indicators over the entire study period. The top panel shows the mean contribution to total stocking intensity by stock type (μDairy, μBeef, μSheep, μDeer; %). The bottom panel shows changes in the contribution of each stock type to total stocking intensity (δDairy, δBeef, δSheep, δDeer; % year^–1 × 100). Categories are defined by six equal quantiles of the distributions of the non-zero values of the indicators over network segments.

Over the 28-year study period, the mean proportion of catchment occupied by plantation forestry (μForest) across both the entire river network and river monitoring sites was 1.4%. Over the study period, the mean change in plantation forestry (δForestry) across both the entire river network and river monitoring sites was 2.8%. Increases throughout the study period in the proportion of catchment occupied by plantation forestry (δForest) were greatest across the North Island and in the eastern and northern regions of the South Island (Fig. 5).

The change in sheep stocking intensity (δSheep) over the entire study period was negative across the entire country, with the largest changes in the Canterbury and Southland regions of the South Island (Fig. 6). The change in dairy cow stocking intensity (δDairy) over the entire study period was positive across much of the country, but particularly in the Waikato region of the North Island and the Canterbury, West Coast and Southland regions of the South Island.

Throughout the study period, dairy cows made the greatest contribution to total stocking intensity (μDairy) in the Waikato and Taranaki regions of the North Island and parts of the West Coast, Canterbury and Southland regions of the South Island (Fig. 6). Sheep made the greatest contribution to total stocking intensity (μSheep) over much of the eastern side of both islands and the Manawatū region in the North Island. Beef made the greatest contribution to total stocking intensity (μBeef) in the Northland and Gisborne regions of the North Island and the southern part of the West Coast of the South Island.

The first 10 components of the PCA used to summarise the land use indicators were significant and explained 97% of the total variation in the land use indicators. The number of land use indicators with high loadings (i.e. absolute values >0.2) were generally limited on the individual PCA axes, enabling meaningful interpretation of all axes (Table 3). For example, sites with high scores on PC1 had high mean dairy cow stocking intensity (μDairy), low mean sheep stocking intensity (μSheep) and high mean total stocking intensity. Similarly, sites with high scores on PC2 had increasing dairy cow (δDairy) and decreasing sheep (δSheep) stocking intensity, and sites with high scores on PC10 had increasing plantation forest cover and high mean plantation forest cover.

Table 3.  Loadings of the land use indicators on the retained PCA components
Values of loadings with an absolute value >0.2 are shown

Of the 168 fitted models (6 variables, 19 windows of 20-year duration and 9 windows of 10-year duration), 157 (93%) were significant after adjusting for false discovery rate. The variation in τ_w explained by the linear regression models was generally low (mean R² = 0.15), varied between time windows and was slightly greater for the 20-year duration windows (mean R² = 0.18) than the 10-year duration (mean R² = 0.14; Table 4). The NO₃-N models had consistently high R² values across windows for both the 10- and 20-year durations (means of 0.17 and 0.31 respectively). The TP model had the lowest R² values for the 10-year duration (mean R² = 0.11) and the NH₄-N model had the lowest R² values for the 20-year duration (mean R² = 0.09).

Table 4.  Summary of explained variation and significance for models representing each time window for durations of 10 and 20 years
‘Explained variation’ is the mean R² over all models (first value), with the minimum and maximum R² values provided in parentheses. ‘Percentage significant’ is the proportion of models that were significant at α = 0.05 after adjusting for false discovery rate

Mapped patterns in the residuals of all models were random and did not suggest the presence of additional environmental drivers of trends that were missing from the models. The mean Mantel’s r value across all models was 0 (maximum absolute value 0.14) and only 8% of Mantel’s r values were statistically significant after adjusting for false discovery rate (for more information see Fig. S3 and S4 of the Supplementary material).

Across all models, confidence in the direction of the relationships between the explanatory variables and τ_w was highest for r_o with a mean confidence of 0.88. For the 10-year duration, r_o had the highest mean confidence (0.90), but for the 20-year duration r_o had the fourth highest mean confidence (0.84). The sign of the coefficient fitted to r_o (β₁) was generally consistent with the direction of the SOI trend (δSOI_w) for each time window (Fig. 7). For the 10-year duration time windows, consistency for r_o was >0.75 for all water quality variables and, when only time windows with absolute δSOI_w values >0.75 were considered (i.e. when time windows with weak SOI trends were excluded), all variables had β₁ consistency values of unity. For the 20-year duration, when only time windows with absolute δSOI_w values >0.75 were considered, all water quality variables except CLAR had β₁ consistency values of unity. With two exceptions (CLAR and TN for the 20-year duration), mean confidence in all fitted β₁ values was >0.75.

Fig. 7.  Relationship between confidence in the sign of the coefficient fitted to r_o (β₁) and the direction of the SOI trend (δSOI_w). The grey boxes indicate windows for which the sign of the fitted β₁ coefficient was consistent with δSOI_w and confidence was >0.9.

More relationships between the synthetic land use indicators and τ_w were consistent and certain for the 20-year duration (26 of 60 coefficients) than for the 10-year duration (4 of 60 coefficients; Fig. 8). For the 10-year duration there were more relationships between the synthetic land use indicators and τ_w with high mean confidence (i.e. >0.75) but low consistency (29) than for the 20-year duration (8).

Fig. 8.  Summary of the confidence in the direction of the relationships between the synthetic land use indicators and site τ_w values for time window durations of 10 years (left column) and 20 years (right column). The individual panels represent each water quality variable and, within each panel, each end year (x-axis) represents a model pertaining to a window and each row represents a synthetic land use indicator. The colour of the cells indicates the direction and level of confidence in the fitted relationship. Blue and red shades represent negative and positive relationships respectively. Colour intensity indicates the level of confidence. The alphabetic codes shown in the key indicate the following confidence quantities: a ≥ 0.99; 0.99 < b ≥ 0.95; 0.95 < c ≥ 0.9; 0.9 < d ≥ 0.75; e < 0.75. The solid horizontal rectangles indicate the indicators that made contributions to the models of each time window duration that were consistent (consistency ≥0.9) and certain (mean confidence ≥0.75).

For all water quality variables, three or more of the synthetic land use indicators made consistent and certain contributions to the 20-year duration models and, for NO₃-N, two and six explanatory variables made consistent and certain contributions to the 10- and 20-year duration models respectively (Table 5). Within water quality variables, the direction of the consistent and certain associations between the synthetic land use indicators and water quality trends differed. For example, for the 20-year duration, CLAR had negative and positive associations with PC1 and PC2 indicating degradation and improvement respectively. Between water quality variables, the direction of the consistent and certain associations between the synthetic land use indicators differed. For example, for the 20-year duration, NO₃-N and TP had positive and negative associations with PC1 indicating degradation and improvement respectively.

Table 5.  Summary of consistent and certain associations between the synthetic land use indicators and water quality trends
‘Impact’ indicates whether the direction of the indicator’s contribution to the trend represents degradation or improvement in water quality

Although river water quality trends are routinely analysed and reported in New Zealand, there has been minimal effort to date to rigorously attribute these trends to environmental drivers (e.g. Hamill and McBride 2003). Accurate attribution of water quality degradation to drivers is a prerequisite for effective water management (Ryberg et al. 2018). In the case of agricultural land use, reducing adverse effects on water quality without unnecessarily prohibiting agricultural practices requires accurate identification of the land uses and practices that cause water quality degradation (Ryberg et al. 2018). To the extent possible with the available data, we attributed water quality trends to the combination of an indicator of climate variability (i.e. r_o) and multiple aspects of pastoral land use. Attribution was based on statistical models that included multiple alternative drivers and consideration of the consistency, certainty and physical plausibility of the associations.

Our model results indicate that SOI trends are associated with trends in the six water quality variables at the 10-year timescale and, to a lesser degree, at the 20-year timescale. The lack of consistent and certain associations between land use indicators and water quality trends at the 10-year timescale indicates that land use signals were generally swamped by the noise of climate variability (Choquette et al. 2019). However, at the 20-year timescale, land use signals were more discernible and a wide range of consistent, certain and physically plausible land use–water quality trend associations were apparent. Confidence in these conclusions was provided by the consistency of the contribution of the land use indicators to the explanation of between-site differences in trends across time windows of 20-year duration. This confidence was increased by the fact that each window corresponded to a unique set of climate conditions (Fig. 3) and a unique set of sites (Fig. 4). Here we consider a subset of these associations and discuss their possible explanations.

The strongest relationships between water quality trends and land use indicators were for NO₃-N. Independent estimates indicate that between 1990 and 2012, nitrogen leaching from agricultural land increased at an average rate of 1.2% year^–1 (Ministry for the Environment and Statistics New Zealand 2019). The estimated increase in leaching rates is strongly associated with the replacement of the formerly dominant sheep industry by dairy farming and with increasing stocking intensity in general (Dymond et al. 2013). These observations are consistent with our attribution of degrading trends in NO₃-N to increasing dairy cow and decreasing sheep stocking intensity (PC2) and to increasing total stocking intensity (PC4) for both the 10- and 20-year durations. Our results are also consistent with catchment-scale studies of the effects of conversion to dairy farming and increasing total stocking rates on river nitrogen concentrations (Monaghan et al. 2007, 2009; Wilcock et al. 2013; Wright-Stow and Wilcock 2017).

We also attributed degrading NO₃-N trends to increasing plantation forest cover (PC3). A plausible explanation for this association relates to the effect of forestry operations on diffuse nitrogen loss rates. Catchment-scale studies indicate that NO₃-N and TN concentrations increase in response to forest planting, thinning and harvesting operations (Bäumler and Zech 1999; Gravelle et al. 2009; Hughes and Quinn 2019).

Improving trends in river DRP and TP concentrations over the past two decades have been reported from several recent national-scale studies (e.g. Larned et al. 2016; Ministry for the Environment and Statistics New Zealand 2017). These improvements may be attributable to several factors; McDowell et al. (2019) suggested that the most probable driver is the growing use of mitigation measures to reduce the loss of phosphorus from agricultural land (e.g. shifting from high- to low-solubility fertilisers). We attributed improving DRP trends to high mean dairy cow and low mean sheep stocking intensity (PC1 in Table 5). A plausible explanation for this may be that the implementation of mitigation measures to reduce phosphorus loss has been growing rapidly on dairy farms, which have had the greatest pressure to improve land use management practices under the Dairying and Clean Streams Accord and its successors (Bewsell et al. 2007; Scarsbrook and Melland 2015; Monaghan et al. 2021). The reasons for degrading DRP trends in association with high mean beef cow and low mean dairy cow stocking intensity (PC8) and degrading TP trends in association with increasing total stocking intensity (PC4) are unclear, but may reflect slower implementation of mitigation measures on non-dairy farms.

A plausible explanation for improving trends in river NH₄-N in association with increasing dairy cow stocking intensity (PC2) is the improved management of farm dairy effluent, particularly as dairy farms have intensified and have often been required to upgrade their infrastructure. Over the past 25 years, there has been a widespread shift from discharging dairy effluent directly to streams and drains to land application. This management change appears to have reduced effluent-derived inputs of NH₄-N to rivers (Monaghan et al. 2010). By contrast, degrading NH₄-N trends associated with increasing total stocking intensity (PC4) may reflect changes in production and practices that have overtaken any benefits associated with the adoption of mitigating strategies over the past two decades.

Two results of this study are indicative of the predictable effects that climate variability has on the direction and strength of water quality trends for time windows of 10- and 20-year duration. First, there was generally high confidence in the coefficients fitted to r_o for all water quality variables and both time durations (Fig. 7). Second, for models corresponding to individual time windows, the signs of the coefficients fitted to r_o were generally consistent with the signs of the corresponding SOI trend (δSOI_w; Fig. 7). The more even distribution of the signs of the coefficients fitted to r_o for durations of 10 compared with 20 years was consistent with the oscillations in strength and direction of water quality trends and the direction of δSOI_w over the shorter timescale. These results are consistent with relationships between SOI and water quality trends reported in Snelder et al. (2021).

Although the contribution of r_o to the models was generally consistent and certain, its inclusion as an explanatory variable should not be construed as comprehensively accounting for the effect of climate variability on water quality trends. ENSO strength has a strong effect on New Zealand’s climate, but it accounts for less than 25% of the year-to-year variance in seasonal rainfall and temperature at most New Zealand measurement sites (Salinger and Mullan 1999). In future investigations, alternative indicators of climate variability to SOI, such as measures of rainfall and temperature, could be assessed as in terms of the variation they explain in water quality trend direction and strength. We anticipate that measures of rainfall and temperature that best explain trend direction and strength will differ between water quality variables due to differences in how the associated contaminants are mobilised, transported, stored and transformed within catchments.

We did not investigate the mechanisms underlying the observed correlations between water quality observations at each site and the SOI (i.e. r_o). In a previous investigation, Scarsbrook et al. (2003) identified geographic patterns in site r_o values for some water quality variables that are consistent with the known impact of ENSO variation on rainfall and stream flow. These patterns suggested that the associations between water quality observations and the SOI at each site, as well as between r_o and water quality trends across sites, are due, in part, to the influence of rainfall on the mobilisation, storage and transport of contaminants. As a corollary, the confounding effect of climate variability may be reduced by apportioning water quality trends into components attributable to changes in stream flow (i.e. discharge trend) and to changes in the concentration–discharge relationship. The recently developed Weighted Regressions on Time, Discharge and Season (WRTDS; Hirsch et al. 2015) performs this type of apportionment and may allow water quality trends to be more accurately attributed to land use (Choquette et al. 2019).

The statistical models used in this study explained a small proportion of the variability in site trends (i.e. τ_w; Table 4). This is an expected result for at least four reasons. First, the SOI is an imprecise proxy for the role of climate variability on water quality. Second, variation in water quality responses to climate variability is mediated by multiple biogeochemical and hydrological processes, including contaminant mobilisation, transport, storage, attenuation, hydrological connectivity and dilution (Gascuel-Odoux et al. 2010; Mellander et al. 2018). Differences in these processes between catchments will produce inconsistent water quality responses between sites to the same climate forcing and weakens our ability to detect the contribution of agricultural land use to water quality trends (Diamantini et al. 2018; Mellander et al. 2018). Third, there has been significant promotion and uptake of mitigation measures to reduce the adverse effects of pastoral agriculture on water quality in New Zealand over the past three decades (Monaghan et al. 2021). These measures and others, such as improved management of point source discharges, are not explicitly represented by our models and are likely to have had spatially variable effects, further contributing to the unexplained variation. Fourth, the resolution of our land use indicators was limited so that land use drivers were represented by the mean of, and changes to, catchment land use. However, catchment averages may not represent the most influential drivers of water quality changes: land use at specific locations may be more influential than catchment averages. Increased explanation of between-site trends may be possible if future research is able to find more precise descriptors of both the climatic and land use drivers of water quality.

Our study indicates that climate variation strongly contributes to trends in water quality and, at the 10-year timescale, almost completely overwhelmed the ability to detect the signal of land use in water quality trends. At the 20-year timescale, our study was able to detect the signal of agricultural land use and land use changes on water quality trends in New Zealand. Our rigorous approach to attributing water quality trends to land use shows that associations are complex and depend on both the water quality variable and the specific land uses and intensity levels. These results are somewhat at odds with the contention that the expansion and intensification of dairy farming has negatively affected water quality over the past 25 years (Foote et al. 2015). Of the 30 consistent and certain associations between agricultural land use indicators and water quality trends, nine were associated with water quality improvement. There was a clear signal of degrading TN and NO₃-N associated with the intensification of dairy farming (PC2). However, dairy intensification was also associated with improving CLAR and NH₄-N and high mean dairy (PC1) was associated with improving DRP. Water quality state is strongly associated with the proportion of catchment area occupied by agricultural land and therefore improving trends associated with agriculture may represent slow rates of improvement from a poor state. However, the association between current water quality state and agricultural land use is not evidence that ongoing changes in that land use are driving water quality degradation. Our study shows that considerable analytical effort is required to extract an understanding of the drivers of water quality trends. The effort to better link water quality changes to land use will enable more precise and effective management actions to be prescribed to arrest degradation in the future.

The authors declare that they have no conflicts of interest.

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