River water quality changes in New Zealand over 26 years: response to land use intensity

Relationships between land use and water quality are complex with interdependencies, feedbacks, and legacy effects. Most river water quality studies have assessed catchment land use as areal coverage, but here, we hypothesize and test whether land use intensity – the inputs (fertilizer, livestock) and activities (vegetation removal) of land use – is a better predictor of environmental impact. We use New Zealand (NZ) as a case study because it has had one of the highest rates of agricultural land intensification globally over recent decades. We interpreted water quality state and trends for the 26 years from 1989 to 2014 in the National Rivers Water Quality Network (NRWQN) – consisting of 77 sites on 35 mostly large river systems. To characterize land use intensity, we analyzed spatial and temporal changes in livestock density and land disturbance (i.e., bare soil resulting from vegetation loss by either grazing or forest harvesting) at the catchment scale, as well as fertilizer inputs at the national scale. Using simple multivariate statistical analyses across the 77 catchments, we found that median visual water clarity was best predicted inversely by areal coverage of intensively managed pastures. The primary predictor for all four nutrient variables (TN, NOx , TP, DRP), however, was cattle density, with plantation forest coverage as the secondary predictor variable. While land disturbance was not itself a strong predictor of water quality, it did help explain outliers of land use–water quality relationships. From 1990 to 2014, visual clarity significantly improved in 35 out of 77 (34/77) catchments, which we attribute mainly to increased dairy cattle exclusion from rivers (despite dairy expansion) and the considerable decrease in sheep numbers across the NZ landscape, from 58 million sheep in 1990 to 31 million in 2012. Nutrient concentrations increased in many of NZ’s rivers with dissolved oxidized nitrogen significantly increasing in 27/77 catchments, which we largely attribute to increased cattle density and legacy nutrients that have built up on intensively managed grasslands and plantation forests since the 1950s and are slowly leaking to the rivers. Despite recent improvements in water quality for some NZ rivers, these legacy nutrients and continued agricultural intensification are expected to pose broad-scale environmental problems for decades to come.

Abstract. Relationships between land use and water quality are complex with interdependencies, feedbacks, and legacy effects. Most river water quality studies have assessed catchment land use as areal coverage, but here, we hypothesize and test whether land use intensity -the inputs (fertilizer, livestock) and activities (vegetation removal) of land useis a better predictor of environmental impact. We use New Zealand (NZ) as a case study because it has had one of the highest rates of agricultural land intensification globally over recent decades. We interpreted water quality state and trends for the 26 years from 1989 to 2014 in the National Rivers Water Quality Network (NRWQN) -consisting of 77 sites on 35 mostly large river systems. To characterize land use intensity, we analyzed spatial and temporal changes in livestock density and land disturbance (i.e., bare soil resulting from vegetation loss by either grazing or forest harvesting) at the catchment scale, as well as fertilizer inputs at the national scale. Using simple multivariate statistical analyses across the 77 catchments, we found that median visual water clarity was best predicted inversely by areal coverage of intensively managed pastures. The primary predictor for all four nutrient variables (TN, NO x , TP, DRP), however, was cattle density, with plantation forest coverage as the secondary predictor variable. While land disturbance was not itself a strong predictor of water quality, it did help explain outliers of land use-water quality relationships. From 1990 to 2014, visual clarity significantly improved in 35 out of 77 (34/77) catchments, which we attribute mainly to increased dairy cattle exclusion from rivers (despite dairy expansion) and the considerable decrease in sheep numbers across the NZ landscape, from 58 million sheep in 1990 to 31 million in 2012. Nutrient concentrations increased in many of NZ's rivers with dissolved oxidized nitrogen significantly increasing in 27/77 catchments, which we largely attribute to increased cattle density and legacy nutrients that have built up on intensively managed grasslands and plantation forests since the 1950s and are slowly leaking to the rivers. Despite recent improvements in water quality for some NZ rivers, these legacy nutrients and continued agricultural intensification are expected to pose broad-scale environmental problems for decades to come.

Introduction
River water quality reflects multiple activities and processes within its catchment, including geomorphic processes, vegetation characteristics, climate, and anthropogenic land uses (Brierley, 2010). Relationships between water quality and these catchment characteristics are not straightforward because all of these factors interact over both space and time. For example, if intensive livestock grazing occurs on steep slopes, surface runoff and consequently river turbidity is expected to be greater than if grazing occurs on flatter areas; in other respects, if fertilizers are heavily applied to sandy soils with high drainage density, rivers will likely become eutrophied over a period of decades due to legacy nutrients slowly leaking to the rivers through groundwater (McDowell et al., 2008). The influence of land use on water quality has also been shown to vary among different climates (Larned et al., 2004). With all of the various types of intensive land uses that have occurred across diverse landscapes over hundreds of years, rivers with degraded water quality are now widespread.
Historically, water quality in rivers was managed to meet minimally acceptable standards or maximum pollutant load limits (Baron et al., 2002;Boesch, 2002;Howard-Williams et al., 2010). However, in the last decade, a greater emphasis has been placed on maximizing the ecosystem services provided by healthy rivers, which is driving efforts to further improve water quality (Brauman et al., 2007;Davies-Colley, 2013). Early efforts in developed countries to improve water quality focused on point-source pollution, particularly wastewater discharges from factories and treatment plants (Campbell et al., 2004). While the broad-scale reduction in point-source pollution elevated many water quality variables above minimal standards, most rivers globally still have water quality impairments due to diffuse pollution from fine sediments, nutrients, and other contaminants (Vorosmarty et al., 2010). Although considerable effort has been directed at monitoring and reducing diffuse pollution with some success, the legacy of pollutants from various land uses remains (Boesch, 2002;Kronvang et al., 2008;Zobrist and Reichert, 2006). Agricultural land uses are by far the greatest contributors of diffuse pollution globally (Foley et al., 2005;Vitousek et al., 1997); however, the "intangible" sources of diffuse pollution make it difficult to assign cause-and-effect relationships between land use and water quality (Campbell et al., 2004).
Many studies have used theoretical or numerical models to examine relationships between land use and water quality because of the lack of consistent water quality monitoring over long periods (bracketing land use change). While modeling approaches can be useful for catchments where much is known about its landscape, modeling may not work well for large, heterogeneous catchments because land-water relationships are complex with interdependencies, feedbacks, and legacy effects. Empirical studies can shed light on some of these complexities, but they are only useful for their particular catchments and may have limited generality or transferability. Comparisons of many diverse catchments is probably most useful to advance understanding of broad-scale landwater relationships (Zobrist and Reichert, 2006).
One of the most comprehensive empirical multi-catchment studies to date on land use-water quality relationships has been the study by Varanka and Luoto (2012) of 32 boreal rivers in Finland. They analyzed five water quality variables over 10 years as a function of a suite of physiographic, climate, and land use variables. A similar study was conducted on many of the same rivers in Finland, but with a more sophisticated temporal analysis (Ekholm et al., 2015). In a study of 11 Swiss watersheds, Zobrist and Reichert (2006) analyzed export coefficients of six water quality variables from biweekly, flow proportional, composite samples over a 24-year period within the context of land use.
All of these studies, and most catchment land use studies, assessed land use (or land use change) as areal coverage. However, land use intensity -the inputs (e.g., fertilizer, livestock) and activities (e.g., vegetation removal) of land usecould be a better predictor of environmental impact for being a more direct measure of impact than land use alone (Blüthgen et al., 2012;Ramankutty et al., 2006). Unfortunately, our understanding of the patterns, processes, and impacts of land use intensity is inadequate because of (1) its complex, multidimensional interactions with other landscape variables, and (2) the lack of appropriate datasets across broad spatiotemporal scales Erb et al., 2016). New Zealand (NZ) provides a valuable test bed for the patterns, processes, and impacts of land use intensity because over the past 3 decades pasture area has decreased but livestock densities and fertilizer inputs have increased (MacLeod and Moller, 2006;StatsNZ, 2015). Like Finland and Switzerland, NZ has an extensive long-term river water quality monitoring network, which has allowed for many studies on river water quality state and trends (Smith et al., 1996(Smith et al., , 1997Scarsbrook et al., 2003;Scarsbrook, 2006;) and effects of land use areal coverage (Davies-Colley, 2013;Larned et al., 2004Larned et al., , 2016. However, this dataset has not been assessed as regards changes in land use intensity that have occurred over the same period. Here, we investigate long-term relationships between land use intensity, geomorphic processes, and river water quality in NZ -which provides a particularly valuable case study because (1) it has had one of the highest rates of agricultural land intensification over recent decades (OECD/FAO, 2015) and thus serves as a potential indicator for countries that are also increasing agricultural intensity; (2) it has a long, consistent, and comprehensive national water quality dataset; and (3) it is physiographically diverse. We examined monthly data for a suite of water quality variables that extend over a 26-year period for 77 diverse catchments. We then compared these states and trends of river water quality to landscape data that characterized the catchments' geomorphology, soil properties, and hydro-climatology as well as temporal changes in land use areal coverage and land use intensity, specifically livestock density and land disturbance, defined here as bare soil resulting from vegetation loss. Altogether, these analyses reveal coincident spatiotemporal patterns in land use intensity and water quality over a quarter of a century. Most of our analyses were performed at the catchment scale, which integrates the spatiotemporal changes that are reflected in our water quality measurements and is the most appropriate scale to manage diffuse pollution (Howard-Williams et al., 2010).

Study area
New Zealand is a small island nation (∼ 268 000 km 2 ) located between the South Pacific Ocean to the east and the Tasman Sea to the west. Its two main islands, North Island and South Island, are located between 34 and 47 • S latitude. Being located on the active boundary between the Australian and Pacific plates, NZ's geology and geomorphology are very diverse, including active volcanoes, karst regions, a range of high-fold mountains (the Southern Alps), large coastal plains, and rolling hills across both hard and soft rocks. Being stretched latitudinally, with nowhere more than about 150 km from the sea, between two major ocean waters combined with its topographic variability, NZ also has a diverse climate with regional extremes, including sub-tropical in the far north, temperate in the central North Island, extremely wet on the western side of the Southern Alps (up to 10 m annually), and semi-arid in the rain shadow to the east of the Southern Alps.
New Zealand is the last major habitable landmass to be settled by humans. Eastern Polynesians first arrived around 1300 AD (Wilmshurst et al., 2008). Europeans first arrived in the late-1700s, but large-scale settlement did not begin until the 1840s. Broad-scale agriculture spread shortly after and has been intensifying since. While we address land use changes at the national scale in this study, our water quality analyses focus on 77 diverse catchments across NZ (Fig. 1).

Water quality data
Water quality data were obtained from NZ's National Rivers Water Quality Network (NRWQN), which is operated and maintained by the National Institute of Water & Atmospheric Research (NIWA). This network represents one of the world's most comprehensive river water quality datasets: 13 water quality and 2 biomonitoring variables have been measured monthly (via in situ measurements and grab samples), with supporting flow estimation, from 1989 to 2014 at 77 sites, whose catchments cumulatively drain approximately half of NZ's land surface (Davies-Colley et al., 2011). Further, this dataset has been operationally stable throughout its history, which allows us to calculate trends over this period. For this study, we focused on 11 water quality variables and their coincident flow (Table 1). We did not analyze ammoniacal nitrogen (NH 4 ) because early NH 4 samples were biased high by laboratory contamination (Davies-Colley et al., 2011).
All water quality variables, except water temperature (T w ), were flow normalized (for each site separately) in JMP ® Pro (v 11.2.1) with local polynomial regression (LOESS) using a quadratic fit, a tri-cube weighting function, a smoothing window (alpha) of 0.67, and a four-pass robustness to minimize the weights of outliers (Cleveland and Devlin, 1988), where flow-adjusted value = raw value -LOESS value + median value. With LOESS, there is no assumption about the water quality variable's relationship with flow. For example, although visual clarity usually decreases systematically with increasing flow (Smith et al., 1997), algae blooms at low flows can sometimes reduce clarity. LOESS also allowed us to examine relative water quality changes over long periods.
We assessed water quality states and trends with ANZECC (2000) guidelines, which are the 20th percentile of the first decade of the NRWQN record for reference. These guidelines are "trigger values" that when exceeded trigger a management response to protect ecosystem health (Hart et al., 1999). Although these trigger values are not effectsbased standards (which would be difficult to define for the wide variety of NZ ecosystems), they do provide a useful reference for comparing water quality states and trends. Upland and lowland catchments, distinguished by the 150 m elevation threshold, have different guidelines that take into account that lowland rivers are typically more turbid and nutrient rich.

Physiographic data
Water quality metrics and trends were compared to a suite of landscape variables (Table 2). Catchment morphometrics (area, slope, ruggedness) were obtained from a 30 m digital elevation model (DEM) that we rescaled (in order to align with other gridded spatial datasets) from the 25 m DEM produced by Landcare Research (LCR). This 25 m DEM was interpolated from 20 m contours of the national TOPOBASE www.hydrol-earth-syst-sci.net/21/1149/2017/ Change in stock unit density (SUD 2012(SUD −1990 Difference between SUD in 2012 and 1990 (SU ha −1 ) Statistics NZ (territorial authority) digital topographic dataset supplied by Land Information NZ (LINZ; scale: 1 : 50 000). Catchment area (A) is the drainage area (in km 2 ) above the NRWQN station, derived using Arc Hydro tools in ArcGIS 9.3.1 in combination with the River Environment Classification (REC, v2.0), the national hydrography dataset derived from a 30 m hydrologically correct DEM (Snelder et al., 2010). Mean catchment slope (S c ) was derived from the same software package, using a 3 × 3 cell window. We defined ruggedness (R r ) as the standard deviation of the 30 m slope grid for each catchment (Grohmann et al., 2011). Drainage density (D d ) was calculated from the ratio of the total length of REC streams to catchment area (in km km −2 ). Soil data were obtained from the 1 : 63 360 Fundamental Soils Layers (FSL), which is maintained by LCR. Methods and data descriptions for this soils database are described in Webb and Wilson (1995) and Newsome et al. (2008). Catchment-scale soil variables (mean value across catchment) that we included in our analysis for being expected to be related to water quality were soil depth (Z s ), percent of catchment dominated by silty and clayey surface soils (SC%), soil pH (pH s ), cation exchange capacity (CEC), organic matter percentage (OM%), and phosphate retention (P ret ). Phosphate retention is a measure (in %) of the amount of phosphate that is removed from solution by the soil via sorption (Saunders, 1965). Thus, soils with high P ret have low P availability for plant growth.
Median annual precipitation (MAP), median annual temperature (MAT), and median annual sunshine (MAS) averaged across each catchment was obtained from NIWA's National Climate Database, which contained 5 km gridded daily weather data (Tait and Turner, 2005). Our values for these three variables represent the median annual precipitation (total mm yr −1 ), temperature (mean • C), and sunshine (hours yr −1 ) for the period 1981-2010. Relative water storage (RWS) was calculated as the proportion of the annual catchment water yield (i.e., total volume of water leaving the catchment in a year) stored in lakes and reservoirs. Reservoir/lake storage was obtained from the Freshwater Ecosystems of NZ (FENZ) database, described in Snelder (2006). The last hydro-climatological variable we included in our analyses was the median discharge (Q 50 ), which was calculated from the NRWQN "flow stamping" at times of water quality sampling from 1989 to 2014.

Land use and intensity data
There are two national land use datasets for NZ. The Land Use and Carbon Analysis System (LUCAS) was developed by the NZ Ministry for the Environment (MfE, 2012) for reporting and accounting of carbon fluxes and greenhouse gas emissions, as required by the United Nations Framework on Climate Change and the Kyoto Protocol. Accordingly, LU-CAS uses 1990 as its reference year and maps land use in 12 classes for 2008 and 2012. The Land Cover Database (LCDB) was developed by LCR, with contributions from MfE, Department of Conservation, Ministry for Primary Industries, and Regional Councils (LCR, 2015). LCDB contains 35 land use classes for 1996, 2001, 2008, and 2012. Both datasets use a minimum mapping area of 1 ha, and use many of the same data and methods to map land use. There are however, some key differences in their class designations and classifications that are important to our analyses: (1) LU-CAS includes Manuka/Kanuka as forest, whereas LCDB designates Manuka/Kanuka as shrub; (2) LUCAS lumps all post-1989 forests into one class, whereas LCDB differentiates between indigenous and plantation forests; (3) LUCAS uses a conservative approach to map high-producing grasslands, whereas LCDB uses phenological information to provide more accurate estimations of high-producing grassland. Because of our focus on (water quality-impacting) plantation forests and high-producing grasslands, we used the LCDB (v4.1) for the midpoint year 2001 for our spatial and statistical analyses. We used LUCAS only to quantify long-term changes from 1990 to 2012, before the LCDB was initiated in 1996. Table 3 describes the land use classes we used in this research, which classes are included from both datasets, and the national comparison between LUCAS and LCDB for 2012.
There are numerous metrics for land use intensity . At the catchment scale, we used livestock density as a metric for all grasslands; and we used land disturbance, defined here as bare soil resulting from vegetation loss, as a metric for high-producing grasslands and plantation forests. We also used national-scale annual fertilizer data  from StatsNZ (2015) to compare longterm trends of river nutrient concentrations to nutrient inputs. Livestock numbers for dairy cattle, beef cattle, sheep, and deer (at 1 ha resolution) for each catchment were derived from maps provided by Ausseil et al. (2013), which is representative for the year 2011. To assess total livestock impact, we multiplied each livestock type by its AgriBase stock unit (SU) coefficient: sheep = 0.95 SU, deer = 1.9 SU, beef cattle = 5.3 SU, and dairy cattle = 6.65 SU (Woods et al., 2006). The total SU for each catchment was then normalized by total catchment area, expressed as stock unit density (SUD) in SU ha −1 .
Changes in SUD from 1990 to 2012 (SUD 2012(SUD −1990 ) were assessed using district-level data from StatsNZ (2015) on total numbers of sheep, deer, beef cattle, and dairy cattle. These livestock numbers were then aggregated for each catchment and multiplied by their respective SU coefficient. Stock unit densities were then compared between 1990 and 2012 to assess change in livestock intensity in each catchment. For Whakatane and Kawerau districts, 1993 was used because 1990 data were unavailable.
Land disturbance (i.e., bare soil resulting from vegetation loss) was quantified for all high-producing grasslands (D HG ) and plantation forests (D PF ), as well as the whole catchment (D C ) for the period 2000-2013. The methods for calculating and validating disturbance are described in de . Briefly, MODIS BRDF corrected reflectance data (MCD43A4) at 463 m spatial resolution and 8-day temporal resolution was used to calculate Tasseled Cap brightness, greenness, and wetness based on the coefficients following Lobser and Cohen (2007). These indices consist of linear combinations of all seven MODIS reflectance bands to represent general image brightness which is comparable to albedo, image greenness which is comparable to the better known vegetation indices such as NDVI and EVI, and image wetness which is linked to the amount of water captured in the vegetation, most comparable to normalized difference water indices. Missing pixels were ignored. We then calculated the mean and standard deviation of each tasseled cap index for each combination of land cover class (LCR, 2015) and climatic region for each 8-day time period. We then used these measures to standardize the calculated tasseled cap indices. To determine how disturbed each pixel was at any point in time, we then calculated the forest and grassland disturbances. The forest disturbance index is calculated as the standardized brightness minus the standardized greenness and wetness. The idea is that disturbed forests appear brighter and less green and less wet than undisturbed forests. The grassland index is the negative sum of all indices, indicating that disturbed grasslands appear darker, less green, and less wet than undisturbed grasslands. MODIS disturbance data were visually validated against 7500 random pixels from Landsat imagery and corresponding 15 high-resolution Orbview-3 and Ikonos images. The overall accuracy of the disturbance index based on Landsat data was 98 %.

Statistical methods
We used non-parametric Spearman rank correlation coefficients (r s ) to look at relationships between variables because many of the relationships were curvilinear. Statistical significance was taken to be an alpha of 0.05. Bivariate comparisons between all variables (Tables 1-3) were performed to explore for associations and identify correlated variables before later multivariate analyses. Median values (from the 26year monthly time series) for water quality variables at each site were used when compared to physiographic and land use variables of their corresponding catchment. Stepwise regression was then used to rank order the relative contributions of multiple landscape variables associated with each major water quality variable. Stepwise regression was used because it accounts for correlations among the independent landscape variables. The order of variables in the stepwise regression model and the sign of their coefficient (proportional [+] vs. inverse [−]) provides an objective measure of the contribution of each landscape variable to river water quality. The level of entry into the model was set to p = 0.05. All the above statistical analyses were performed in JMP ® Pro (v 11.2.1). Temporal trends in flow-normalized water quality  and disturbance (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) data were assessed with the seasonal Kendall (SK) test, which was corrected for temporal autocorrelation using the rkt R package; missing values were ignored. We also calculated the SK slope estimators (SKSEs) using the same R package. Because some NR-WQN sites had multiple measurements in some months, a few records (no more than five) were removed from each site in order to ensure 12 monthly values for each year for the SKSE test. There were also occasional missing values for some variables throughout the time series, particularly in the early years. Of particular note, there were no TN values for 1994 as a result of contamination by leaking ammonia refrigerant during storage of frozen subsamples. HV1 did not have data for 18 months from 2012 to 2014.
In order to make trend comparisons among sites and derive an estimate of percent change per year, we normalized SKSE values by dividing them by the raw data median to give the relative SKSE (RSKSE) in percent change per year (Smith et al., 1996). Given that water temperature (T w ) uses an arbitrary scale in • C, we only report SKSE values for this variable. We also used the trend categories of Scarsbrook (2006): (1) no significant trend -the null hypothesis for the SK test was not rejected (p > 0.05); (2) significant increase/decrease -the null hypothesis for the SK test was rejected (p < 0.05); and (3) "meaningful" increase/decrease -the trend was significant and the magnitude of the trend (RSKSE) was greater than 1 % per year. A 1 % change per year translates to slightly more than 10 % change per decade (due to compounding), a rate of change that is easily detectable and observable.

Physiographic characteristics
The 77 NRWQN catchments were physiographically diverse in terms of morphometric, soil, and hydro-climatological variables (Table 4; Table S1 in the Supplement). Most notable with regards to its direct influence on runoff and water quality was median annual precipitation (MAP), which ranged from 533 to 7044 mm yr −1 . When combined with the wide range of catchment areas (A), median discharge (Q 50 ) varied over 3 orders of magnitude, from 0.4 to 515 m 3 s −1 , and annual water yield from 103 to 3475 mm yr −1 . In terms of soil, about a quarter of the catchments had very sandy surface soils (SC% < 10) and a quarter had fine-textured soils (SC% > 70). Phosphate retention (P ret ), an important variable for fertilizer management and consequently water qual-ity, was particularly high (> 57 %; 10th percentile) for seven catchments in the central North Island. Several physiographic variables (Table 2) displayed strong latitudinal trends from north to south and many were strongly correlated (p < 0.001; Fig. S1 in the Supplement). In consideration of these relationships and perceived importance for water quality (sensu Varanka and Luoto, 2012), we used the following subset of minimally correlated physiographic variables for subsequent multivariate analyses: catchment slope (S c ), silt-clay percentage (SC%), phosphate retention (P ret ), and median flow (Q 50 ).

Land use areal coverage and temporal changes
Land use in NZ, like physiography, varied widely, and our 77 catchments captured this diversity ( Fig. 1; Table S2). In 2001, 13 catchments were dominated by non-plantation forests (NF), while three catchments were dominated by intensively managed plantation forests (PFs); 13 catchments were dominated by shrub/grassland (SG) that was not intensively managed. The most dominant land use was grasslands that were intensively managed (high-producing grasslands; HG), covering the majority of the area for 31 catchments. Open water (OW) was the majority land use for only one catchment and relatively high (> 10 %) for two others. Barren/other (BO), which was largely bare rock, was relatively high (> 10 %) for 13 mountainous catchments. Urban (UR) coverage rarely exceeded 1 %, with only one catchment greater than 2 %. Annual cropland (AC) exceeded 1 % in 11 catchments, but never exceeded 8 %. Vegetated wetland (VW) and perennial cropland (PC) were minimal in all catchments, each rarely exceeding 1 %.
In general, NF, SG, and BO areas dominated mountainous catchments with high S c and low Z s ; while HG dominated most lowland catchments with low S c , high Z s , and high pH s . Like HG, PF mostly occurred on flat areas (r s = −0.48 with S c ) with thick soils (0.35 with Z s ) that were less acidic (0.31 with pH s ). Given the relative dominance of catchment land use, relationships with physiographic variables, and potential effects on water quality in NZ rivers (Davies-Colley, 2013;Howard-Williams et al., 2010), the land use variables used for subsequent multivariate analyses were NF, SG, HG, PF, and OW.
Land use areal coverage did not change much from 1990 to 2012 across NZ (Fig. 2) or in many catchments (Table S2). The greatest change was a 13.4 % increase in PF in GS1, which was almost entirely accounted for by a 13 % decrease in SG; 13 other catchments experienced small increases (3.0-6.6 %) in PF, accounted for by decreases in SG or HG or both. HM3 and HM4 had the greatest increases in HG at 3.4 and 2.0 %, respectively. HG for the other 75 catchments remained virtually unchanged (< 0.4 %) or decreased. WH3 had the greatest decrease in HG at −4.8 %. Land use areal coverage change in other catchments was negligible.

Land use intensity and temporal changes
Changes in total stock unit density between 1990 and 2012 (SUD 2012(SUD −1990 ) were also minor with only two catchments changing more than 1.6 SU ha −1 over this period (Table S3). Temporal changes in SUD 2012−1990 for 56 of the 77 catchments were within the range of −1.0 to 1.0 SU ha −1 . Although land use areal coverage and total livestock densities changed little over the period 1990-2012, livestock types changed considerably for many catchments (Table S3) and across NZ (Fig. 2). The general pattern was dairy cattle replacing sheep. The number of dairy cattle from 1990 to 2012 increased in 72 catchments, with a mean increase of 0.6 SU ha −1 for all catchments, whereas the number of sheep decreased in all 77 catchments (mean = −0.9 SU ha −1 ). Deer and beef cattle numbers changed little: 0.0 and −0.2 SU ha −1 , respectively.
When 2011 livestock densities were compared with physiographic variables, the strongest relationships were found with combined SUD of dairy and beef cattle (hereafter SUD cattle ; Fig. S2). SUD cattle decreased strongly with increasing slope, S c (r s = −0.79), but increased with Z s (0.43), pH s (0.32), and P ret (0.27). SUD cattle also increased with MAT (0.68) and MAS (0.42), but decreased with MAP (−0.34). Thus, the highest cattle densities were found in catchments, such as WA3 (with the highest SUD cattle at 15.7 SU ha −1 ), that were relatively flat, warm, sunny, and dry, with deep soils that had relatively high pH and high P retention. HG had similar, but less strong, correlations with these same physiographic variables.
Catchment disturbance (D C ) varied widely over both space and time between 2000 and 2013 (Table S4). The max- imum amount of D C at one time was 35.7 % for WN3 on 7 April 2003, almost entirely due to bare pastures. D C exceeded 15 % on six other occasions (264 days in total) in this catchment. In general, the North Island (Fig. 3) had a greater extent and intensity of disturbance than the South Island (Fig. 4). The most intense disturbances occurred as a result of plantation forest harvests, and these disturbances were on average visible for about 1.5 years up to about 4 years, with exceptions lasting more than 6 years. Indeed, D C was strongly correlated to PF coverage (r s = 0.51). The catchment with the highest median D C (10.5 %) was RO3, which had 69.8 % of its catchment in PF and 17.7 % in HG; 14 other catchments had D C above 5 %, and two-thirds of these were dominated by either PF or HG.
We also analyzed disturbance of plantation forests (D PF ) and high-producing grasslands (D HG ) separately for each catchment. For catchments with at least 21.4 km 2 (100 MODIS pixels, for the sake of statistical robustness) of plantation forest, the mean (±SD) D PF (from 2000 to 2013) was 10.6 ± 5.6 %. The catchments with the highest D PF were those with low mean annual precipitation, MAP (r s = −0.42). There were no significant relationships between D PF and any of the other physiographic variables. For catchments with at least 21.4 km 2 of high-producing grass- lands, the mean (±SD) D HG was 6.0 ± 6.4 %. The catchments with the highest D HG were those with low mean annual sunshine (MAS; r s = −0.25), low mean annual temperature (MAT; −0.30), high catchment slope (S c ; 0.25), and high ruggedness (R r ; 0.31). The six catchments with the highest D HG (> 15 %) all had low phosphate retention (P ret ; < 32 %). While it is assumed that greater densities of livestock lead to greater pasture disturbance, we did not find a proportional relationship between SUD and D HG among catchments. In fact, the highest median D HG was found for catchments with low SUD (r s = −0.45). Over time, however, we observed a fairly strong trend (r s = 0.50) of lower D HG with decreasing SUD (−SUD 2012(−SUD −1990 ). In all there were seven catchments with significant or meaningful decreases in D HG from 2000 to 2013 (assessed with SKSE), all of which had a negative SUD 2012−1990 .

Water quality characteristics and trends
Median monthly values of water quality variables for the 77 catchments ranged widely (Table 5; Table S5). Some rivers had exceptional water quality all around, while others had either current issues with multiple variables or worsening temporal trends (assessed with SKSE from 1989 to 2014; Table 6). Because of the dependence of water quality on flow, we first assessed temporal trends in Q. Only two catchments had significant increases in Q, with one also being meaningful. Three catchments had significant decreases in Q and five others also had meaningful decreases in Q. Water temperatures (T w ) were not particularly high for any of the catchments; however, 21 rivers had significant increases in T w , possibly the signature of climate change. Because of its strong latitudinal trend (stronger than any land use effect), T w was not analyzed further. Dissolved oxygen (DO) was close to 100 % for most catchments, but was particularly low (< 90 %) for two catchments: one affected by periurban activities (AK2) and one affected by discharge from a large pulp mill (RO2). Temporal trends in DO from 1989 to 2014 were relatively minor (RSKSE < 0.5 % yr −1 ), except for RO2 which had a significant increase attributable to progressive improvements in treatment of organic waste from its large pulp mill. Conductivity (COND) was relatively low (< 115 µS cm −1 ) for all South Island catchments and varied considerably for the North Island (54-528 µS cm −1 ). Most catchments (52/77) experienced significant or meaningful increases in COND from 1989 to 2014. Water pH (pH w ) was neutral to alkaline for all rivers, which have been described as calcium-sodium bicarbonate waters by Close and Davies-Colley (1990), and only displayed minor changes over the 26-year study period.
Median visual water clarity (CLAR) was exceptionally high (> 5 m) for seven catchments and very low (< 1 m) for 22 catchments. Since 1989, CLAR has improved in almost half of the rivers, and has worsened in four rivers (Table 6;  Table S5). Water turbidity (TURB) was strongly inversely proportional to CLAR (r s = −0.97) and generally followed opposite trends of CLAR. Colored dissolved organic matter (CDOM) was low for most of the rivers, with only five catchments greater than 2.0 m −1 ; 19 of the catchments have experienced significant or meaningful decreases in CDOM since 1989, possibly due to the loss of wetlands across NZ. Only one catchment had a meaningful increase in CDOM.
Total nitrogen (TN) was relatively high (> 455 mg m −3 ) for almost a third of the catchments, with the vast majority (17/23) of these being lowland catchments. Most of these catchments also had relatively high NO x ; 33 catchments had significant or meaningful increases in TN from 1989 to 2014, while only five had significant or meaningful decreases in TN (Table 6). NO x not only had a similar number of increasing temporal trends but also had meaningful decreases for 12 catchments. Total phosphorus (TP) followed a similar geographical pattern as TN; 18 of the 23 catchments with relatively high TP (> 30 mg m −3 ) were lowland catchments. Most of the catchments with relatively high TP (18/23) also had relatively high dissolved reactive phosphorus (DRP) (> 9.5 mg m −3 ); 17 catchments had meaningful increases in DRP, compared to only three with meaningful decreases. There was more of a balance in temporal trends of TP, with eight meaningful increases and seven meaningful decreases.
In addition to the expected correlations between CLAR and TURB, and among the nitrogen and phosphorus constituents, several other significant relationships existed among the water quality variables (Fig. S3). Taking into consideration this broad multicollinearity, we focused our multivariate analyses on several key water quality variables, particularly those that experienced the most changes from 1989 to 2014 (Table 6): CLAR, TN, NO x , TP, and DRP.

Water quality relationships with physiography, land use, livestock density, and disturbance
CLAR generally decreased with A (−0.37; all following parentheses in this section are r s unless specified). Except for TURB (0.32), no other water quality variables had significant relationships with catchment area. Several water quality variables correlated with catchment slope (S c ), including TN (−0.72), TP (−0.63), and DRP (−0.65), meaning N and P concentrations were relatively high in lowland (low slope) catchments. DRP (0.65) and TP (0.61) were directly proportional to mean annual temperature (MAT), but this association probably arises because the highest phosphorus values occurred mainly in lowland catchments and some of the northernmost catchments, with temperature being strongly correlated with altitude and latitude. DRP also had a signif- icant relationship with soil phosphate retention, P ret (0.35).
No other strong physiographic relationships emerged from our analyses. The strongest relationships between water quality and land use areal coverage (Table 7) included HG, which had strong positive relationships with several water quality variables except CLAR, which decreased as HG increased. The lessermanaged SG had generally opposite relationships with water quality, but note that SG did not have significant relationships with TURB or CLAR. NF followed the same trends as SG, but had fewer significant relationships with water quality. PF, on the other hand, followed the same trends as HG, with poorer water quality being associated with greater coverage of PF; although correlations were not as strong as HG. CDOM, DRP, and all N constituents had significant negative correlations with OW, meaning that water quality improved with greater OW coverage, plausibly due to entrapment of fine sediment and nutrients.
Water quality was significantly correlated with all SUD metrics (Table 7; Fig. S4), except deer (SUD de ), which only had relatively weak relationships with TN and NO x . The nutrients and CDOM had the strongest correlations with SUD cattle , which includes both dairy and beef cattle. COND, CLAR, and TURB had the strongest (slightly) correlations with SUD be . Overall, degraded water quality was strongly associated with high livestock densities, even stronger than areal coverage of HG.
No significant correlations between water quality and total catchment disturbance (D C ) were found; however, there were significant associations when disturbance was isolated by high-producing grasslands (D HG ) and plantation forest (D PF ; Table 7). Unexpectedly, CLAR and TURB were not correlated to D HG , and surprisingly, the rest of the water quality variables had a significant inverse relationship with D HG . Conversely, CLAR was the only water quality variable correlated to plantation forest disturbance, D PF (r s = −0.27). Some interesting results emerged when temporal trends in water quality (via SKSE) were assessed for catchments with high disturbance. Of the 15 catchments with D c greater than 5 %, six had meaningful increases in TURB; while only one had a meaningful decrease in TURB. Most of these 15 catch-ments also experienced significant increases in TN (9 catchments; 7/9 also meaningful) and NO x (10 catchments; 8/10 also meaningful). Interestingly, TP and DRP significantly increased in only two of these highly disturbed catchments.

Multivariate water quality relationships
In order to build on the above correlation analyses, the water quality variables of CLAR, TN, NO x , TP, and DRP were each assessed in a multivariate stepwise regression, using the following 10 physiographic and land use independent variables: S c , SC%, P ret , Q 50 , NF, SG, HG, PF, OW, and SUD cattle (Table 8). The residual plots for all five water quality variables met the assumptions of normality and linearity, but displayed heteroscedasticity with a wide scatter for high values. CLAR was correlated to −HG, followed by +OW, −Q 50 , and −PF, where signs represent whether the relationship is positive (+) or inverse (−). Thus, water clarity was predictably lower for larger rivers that drain larger areas of high-producing grasslands and/or plantation forests, but improved with increased open-water coverage (Fig. 5).
The combined stock unit density for beef and dairy cattle (SUD cattle ) was the primary predictor for all four nutrient variables, with TN, TP, and DRP also being proportional to PF coverage (Table 8). Dissolved oxidized nitrogen (NO x ) was not proportional to PF, or any other independent variable in the stepwise regression. Coverage of HG and silt-clay surface soils (SC%) were also proportional factors for TN. Whether intensity or areal coverage, land use was the primary and secondary predictor for all five water quality variables (Fig. 5).

River water quality states and trends
We characterized water quality states and trends for 77 river sites across NZ using a wide range of flows and water quality conditions for each site, including some small floods. We acknowledge that our analyses did not fully capture large floods due to their short durations, unlikelihood of occurring during Figure 5. Multivariate relationships between major water quality variables (median value for each site) and land use variables. For each plot, the primary explanatory variable from the stepwise regression (Table 8) is the x axis, with bubble color representing the secondary explanatory variable. Note that oxidized nitrogen (NO x ) did not have a secondary explanatory variable. Selected catchments discussed in the text are labeled. the preset monthly sampling, and the fact that we relied on grab samples. These episodic floods are particularly important for water quality of downstream waters such as lakes and estuaries (Stamm et al., 2014). The uncertainty surrounding our lack of flood samples could have been mitigated by composite samples or supplemental flood samples; however, our 26 years of monthly samples for each site (n = 312) did allow us to confidently report median conditions and temporal trends in water quality (Moosmann et al., 2005). There was a wide range of water quality across NZ rivers (Table 5), with drastic differences between upland catchments and the more intensively managed lowland catchments. Overall, lowland rivers had considerably lower CLAR and higher TURB, TN, NO x , TP, and DRP. Only two (alpine glacial flour-affected) upland rivers were below the ANZECC CLAR guideline of 0.6 m, whereas 17 lowland rivers were below the ANZECC guideline of 0.8 m. Similarly, 13 lowland catchments exceeded the ANZECC TN guideline of 614 mg m −3 , but only eight upland catchments exceeded the much lower guideline of 295 mg m −3 . Almost three-quarters of these catchments (15/21) also exceeded the NO x guideline of 444 mg m −3 (lowland) and 167 mg m −3 (upland). There were a similar number of sites exceeding ANZECC guidelines for TP (33/26 mg m −3 for lowland/upland) and DRP (10/9 mg m −3 for lowland/upland), each with at least 20 and most of these were corresponding.
Our results on the state and trends of the 77 NRWQN catchments generally accord with earlier NRWQN studies (e.g., ) and a recent publication by Larned et al. (2016), which analyzed water quality states and trends for 461 NZ river sites for the period 2004-2013.
Based on ANZECC (2000) trigger values, we have organized the catchments into four classes (Fig. 6): (I) clean river with high visual water clarity (CLAR) and low dissolved inorganic nutrients (DIN), (II) sediment-impacted river with low CLAR and low DIN, (III) nutrient-impacted river with high CLAR and high DIN, and (IV) sediment-and nutrientimpacted river with low CLAR and high DIN. Note that the term "sediment impacted" is a connotation for total suspended solids (TSS), which includes organic matter as well. In agriculture-dominated catchments, both mineral sediment and particulate organic matter can greatly increase TSS (Julian et al., 2008). We use CLAR as a preferred metric for suspended matter because TSS are not routinely measured in the NRWQN (or other monitoring networks) while CLAR correlates strongly to TSS (r = −0.92), and better than TURB (r = 0.87) . Further, CDOM in NZ rivers is low with minimal impact on CLAR. We use NO x as our preferred metric for DIN because it is least affected by suspended sediment and soil properties (compared to DRP). However, catchments that exceed ANZECC guidelines for DRP are indicated in Fig. 6 by gray-filled markers.
When this classification is combined with the SKSE trend analyses (Table 6), we obtain a clear picture of the current and potential state of NZ rivers (Fig. 6). Before individual rivers are discussed, we first point out key differences between the upland and lowland catchments, which will later be placed within the context of physiography and land use intensity. Most obvious, and consistent with the findings of Larned et al. (2004), was that lowland rivers were much more degraded, particularly by sediment. More than a third of the lowland catchments were either class II or IV (17/44), whereas only two upland catchments were class II. None of the upland catchments were class IV, and more than twothirds were clean rivers (class I). Both types had a similar number of nutrient-impacted rivers (class III). Particularly concerning is that almost half of the lowland rivers (19/44) are currently experiencing meaningful increases (> 1 % per year) in NO x , DRP, or both. The other striking trend is that many of the lowland rivers are becoming clearer, with 18/44 experiencing meaningful increases in CLAR -which, plausibly, has been attributed to increasing riparian fencing to exclude cattle from channels (Davies-Colley, 2013;Larned et al., 2016).
While clearer rivers are seen as an improvement in water quality; when combined with increasing nutrients, warmer water, and lower flows, the perfect recipe for toxic algae blooms is created (Dodds and Welch, 2000;Hilton et al., 2006). Only recently has the widespread problem of toxic algae blooms in NZ rivers been evidenced (Wood et al., 2015;McAllister et al., 2016), and our results indicate that this problem could worsen given the increasing trends we found in water temperatures, inorganic nutrients, and most influential in our opinion, water clarity. Nutrient enrichment and global warming receive the most attention when it comes to degraded water quality, but rivers have increasingly become light limited (Hilton et al., 2006;Julian et al., 2013) such that when clarity improves in warm, nutrient-rich rivers, algae can proliferate. Particularly problematic for NZ is that its lowland catchments, which are warmer, have much greater DRP and NO x , and have longer water residence times, are the ones becoming appreciably clearer (Fig. 6). If droughts become more frequent and intense in NZ, toxic algae blooms are also likely to become more frequent, more widespread, and more problematic. However, this algae response is complex and depends on a number of interacting factors such that the apparent potential for increasing algal nuisance might not necessarily be realized in some rivers (Dodds and Welch, 2000;Hilton et al., 2006).

The role of physiography in dictating land use intensity across New Zealand
While physiography did not emerge as a significant independent variable in the multivariate analyses (except TN with SC%), physiography is important because it largely controls the location and intensity of agricultural land uses. The greatest coverages of HG and the highest densities of cattle (SUD cattle ), the two primary explanatory variables for all five major water quality variables (Table 8), were both found predominantly in flat areas with deep soils located in warm, sunny, and relatively dry climates. Livestock in NZ depend almost exclusively on pasture grasses and thus their productivity is maximized when pasture productivity is maximized. The very large cattle are not well suited for steep slopes, particularly dairy cattle, which can weigh more than 500 kg. Deep soils are important because they absorb and hold more water for plant uptake, and are not as susceptible to waterlogging, especially in wetter climates. Year-round and intense grazing is best supported by warm and sunny climates where pasture grasses are highly productive and recover quickly following intense grazing such as strip/rotational grazing, which is common in NZ dairy farms. Another soil property we found to be positively correlated to SUD cattle was phosphate retention (P ret ). The highest dairy cow densities were found on Allophanic volcanic soils with high P ret , likely because these soils respond favorably to P fertilizer and thus can be managed more intensively. However, soils with high P ret require more P fertilizer, and thus generally have higher export of DRP to rivers. Our finding of a significant positive correlation between these two variables is consistent with this interpretation. Further, we found that high-producing pastures with high P ret had the lowest disturbance (D HG ), indicating that these intensively managed pastures recover quickly following grazing. In a more comprehensive study of land disturbance across the North Island Figure 6. River water quality classes for upland (a) and lowland (b) catchments in New Zealand: (I) clean river with high visual water clarity (CLAR) and low dissolved inorganic nutrients (DIN); (II) sediment-impacted river with low CLAR and low DIN; (III) nutrientimpacted river with high CLAR and high DIN; and (IV) sediment-and nutrient-impacted river with low CLAR and high DIN. Classes are organized by ANZECC (2000) trigger values for water clarity (x axis) and NO x (y axis). Catchments that exceed ANZECC guidelines for DRP are indicated by gray-filled markers. Arrows indicate direction of trend over the 26 years inclusive from 1989 if significant (dashed) or meaningful (solid). No arrow means the trend was not significant. of NZ, de Beurs et al. (2016) also found that Allophanic soils had the least disturbance among all soil orders. Where high livestock densities occur in less than ideal conditions, land disturbance is likely. Our catchment-scale analyses limit our interpretation of specific situations, but based on our results, field observations, and previous remote sensing analyses, pasture disturbance in NZ will likely be highest during droughts on steep, south-facing slopes with thin soils being heavily grazed by sheep. Under these conditions, grasses will be grazed down to bare soil and recover very slowly.
Plantation forests (PFs) in NZ also correlated with thick soils with relatively high P ret on flat areas, particularly the pumice soils of the central North Island. The porous nature of the pumice soils allows them to efficiently hold and regulate nutrients, water, and air while being well-drained and resistant to compaction and flooding. Under these conditions, radiata pine (the dominant PF species in NZ) grows rapidly (mean harvest cycle of 28 years) and can be harvested yearround. Since 1990, however, many of the PF additions have occurred on steeper slopes in response to carbon credit incentives, greater economic demand for wood products (PCE, 2013), and the need for soil erosion control on steep pasture susceptible to land sliding (Parkyn et al., 2006).

Land use intensity and water quality in New
Zealand rivers

High-producing pastures and livestock densities
HG coverage was the primary explanatory variable for visual clarity (CLAR; Table 8, Fig. 5). CLAR in NZ rivers is mostly influenced by mineral and organic particulates . Livestock reduce visual clarity in multiple ways, especially in NZ where high densities of multiple types of livestock tread year-round on relatively steep slopes with highly erodible soils vegetated by shallowroot, introduced grasses that are susceptible to destabilization (McDowell et al., 2008). The year-round treading is particularly important because most NZ regions during winter are very wet with short days, which increases soil disturbance (pugging and compaction) and slows recovery times. Where livestock have direct access to rivers, their trampling of riverbanks and instream disturbance is often the main contributor to reduced CLAR (Trimble and Mendel, 1995;McDowell et al., 2008). The lowland flatter areas in NZ have high HG coverage and high cattle stock densities (SUD cattle ). These lowlands also have high drainage densities -often increased by artificial drainage. The influence of HG on CLAR is thus exacerbated by this interaction of high SUD cattle and artificial drainage. Interestingly, SUD cattle was not an explanatory variable for CLAR in the stepwise regression, which is likely a result of two factors. First, HG and SUD cattle are highly correlated, and stepwise regression does not include secondary variables that are explaining the same proportion of variance as the primary independent variable. Second, we found that CLAR has actually improved in catchments where SUD cattle is high and/or has increased (Fig. 6), which we noted earlier could be a result of increased riparian fencing. In 2003, NZ implemented the Dairying and Clean Streams Accord, which has led to the exclusion of dairy cattle from 87 % (as of 2012) of perennial rivers greater than 1 m in width (Bewsell et al., 2007;Howard-Williams et al., 2010;Gunn and Rutherford, 2013). By excluding (dairy) cattle from channels and riparian zones, the contribution of riverbank and bed erosion to degraded CLAR has likely been mitigated and reduced over time (Trimble and Mendel, 1995;Hughes and Quinn, 2014). Indeed, CLAR has been significantly and meaningfully improving in many of NZ's rivers (Table 6), even those with increasing SUD cattle , albeit from a fairly degraded condition.
Another potential explanation for improved water clarity at numerous sites is the considerable decrease in sheep density across the NZ landscape. NZ had 57.65 million sheep in 1990. By 2012, that number had been reduced by almost half, to 31.19 million (StatsNZ, 2015). Although cattle are larger and have a greater treading impact per animal, the much greater number of sheep means that SUD may be broadly comparable as regards environmental impact. Another difference is that sheep are generally placed on steeper, less stable slopes in NZ, where headwater stream channels are located. Where there are breaks in slope (even small ones), sheep create tracks of bare soil with their hooves and hillside scars with their bodies (for scratching and shelter), both of which can enhance soil erosion (Evans, 1997). Further, cattle (using their tongues) leave approximately half the grass height on the pasture after grazing, whereas sheep (using their teeth) graze approximately 80 % of grass height (down to bare soil in dire conditions), leaving it exposed to erosion (Woodward, 1998). Considering all these factors, sheep can have a greater impact on sediment runoff into rivers, and consequently visual clarity, than suggested by their aversion to water vs. cattle's attraction to water. Although not isolated in our analyses, the particulate fractions of TN and TP have likely been affected by similar processes as CLAR and may follow the same temporal trends .
While HG was also strongly correlated to river nutrient concentrations (Table 7), the primary explanatory variable for all four major nutrient metrics (Table 8, Fig. 5) was land use intensity as measured by livestock density of beef and dairy cattle (SUD cattle ). The difference between these two explanatory variables may seem trivial; however, the distinction is important if we want to understand future trends and effectiveness of water quality management strategies. As we demonstrated, the area of land used for HG has not changed much since 1990 (Fig. 2). In fact, it has decreased or stayed virtually the same in all but two of the 77 catchments. Yet, nutrient concentrations have been increasing in many of the rivers (Table 6), which we attribute to (1) increasing numbers of cattle (mostly dairy) on both HG and SG, and (2) legacy nutrients being slowly delivered to the rivers in groundwater.
From 1990 to 2012, NZ approximately doubled its number of dairy cattle, exceeding 6.4 million. (StatsNZ, 2015). This enormous addition to a country that is only 268 000 km 2 in area, has been accompanied by more than 1.426 million tons of P -based fertilizers and 335 000 tons of N-based fertilizers annually (1990StatsNZ, 2015). Of the nutrients consumed by lactating dairy cows, approximately 66 % of P and 79 % of N are returned to the landscape in the form of urine and feces . This results in about 940 000 tons of P -based and 260 000 tons of N-based diffuse pollution, which is an underestimation because cloverrye grass dairy pastures also receive large inputs from fixed atmospheric N (Ledgard, 2001). Some of these nutrients will be transported to rivers during subsequent storms, but a majority will remain (building up) in the landscape to be slowly added to rivers over decadal timescales (Howard-Williams et al., 2010).

Plantation forests
All water quality variables were significantly correlated to plantation forest coverage (PF ; Table 7), with a negative relationship with CLAR but positive for all other variables. From the stepwise regression, PF emerged as an explanatory variable for all major water quality variables except NO x (Table 8), suggesting that its dominant impact on river water quality was from surface runoff. Plantation forestry activities can add a considerable amount of sediment and nutrient pollution to rivers, especially during and immediately following harvesting (Fahey et al., 2003;Croke and Hairsine, 2006;Davis, 2005). This harvesting period of maximum soil disturbance usually lasts about 2 years (Fahey et al., 2003), but the land cover may remain sparsely vegetated and susceptible to erosion for several years (but usually not more than 5 years; de . The greatest PF impact on sediment runoff, and thus potentially CLAR, is usually from road sidecast/runoff, shallow landslides, and channel scouring/gullying (Fahey et al., 2003;Motha et al., 2003;Fransen et al., 2001).
Rivers receive a pulse of nutrients during the forest harvest, but fertilizers are also applied at time of re-planting and sometimes routinely to enhance growth (Davis, 2005). Radiata pine in the pumice soils of the central North Island, the dominant area of PF in NZ, are particularly responsive to both N and P fertilizers and thus likely receive ample supplements. Like pasture fertilizers, some of these nutrients may be delivered to rivers during intense precipitation, but there is also a legacy of nutrients left behind. Fertilizers have been applied to plantation forests in NZ since the 1950s, with an intense period of application in the 1970s (Davis, 2005). While fertilization rates (tons ha −1 yr −1 ) have decreased since 1980, the amount of NO x leaving catchments mostly covered in PF has significantly and "meaningfully" increased since 1989. None of these catchments had more than 17.7 % HG, none had major increases in HG (< 0.3 %), none had major increases in SUD cattle (< 0.7 SU ha −1 ), and none had a significant increase in D PF . What the catchments did have in common were all had gravelly/sandy pumice soils (< 4.5 SC%) and all were intensively managed as reported by Davis (2005) and as indicated by high D C (> 6.8 %). The extended periods of non-vegetated land due to weed control also increases the amount of nutrients delivered to rivers over the long term (Davis, 2005).

Land disturbance and water quality
So far, we have discussed how land use, livestock densities, and fertilizer inputs affect water quality, with a focus on sediment and nutrient runoff. When land is disturbed (i.e., bare soil), sediment/nutrient mobilization can be enhanced. The most intense and longest lasting disturbances occurred during plantation forest harvests. Following harvest, we found that the land remained disturbed for 1-6 years, with a mean of 1.5 years. The overall mean and median D PF among all catchments was 10 %, which means that plantation forestry leaves large areas of disturbed land at any one time. When this bare land is exposed to intense precipitation, large quantities of sediment and nutrients can be mobilized into the rivers. This process has been documented for numerous catchments across NZ (Basher et al., 2011;Hicks et al., 2000;Phillips et al., 2005). Because these disturbances only last a few years, they typically do not show up as temporal trends (via SKSE); however, it is possible that they produce enough readily available sediment to impact water quality for longer periods (Kamarinas et al., 2016).
The coincidence of rainstorms on disturbed pasture could have the same effect on sediment/nutrient runoff if the pasture is connected to the stream network via steep slopes or adjacent channels/canals (Dymond et al., 2010;Kamarinas et al., 2016). Pastures become disturbed from overgrazing, strip grazing, pugging/soil compaction, tilling/reseeding, cropping/harvesting, or landsliding on steep slopes. Given the high intensity of grazing management in NZ, all of these are common. While D HG was lower than D PF on average, D HG had a higher maximum (Table 4). Spatiotemporal patterns in disturbance between these two land uses were also different . D PF covered large areas and lasted years at a time, whereas D HG had two patterns: (1) one related to dairy cattle strip grazing, which were short lived due to quick recovery times of grasses in fertilized soils; and (2) more widespread and longer continuous disturbances occurring on steeper slopes grazed by sheep and beef cattle, particularly following drought periods. Because our disturbance analyses had a spatial resolution of 463 m, we likely missed some paddock-scale disturbances. Future work could use Landsat imagery (30 m resolution) to assess disturbance (sensu de . All six catchments with meaningful increases in D HG had large increases in dairy cattle density 1990-2012 (Tables S3, S4). Not surprisingly, all six catchments suffered impacts to water quality. Five of the six had meaningful increases in DRP and three had meaningful increases in NO x and TN. One had a meaningful increase in TURB and three had significant reductions in DO. One of these catchments, in particular, may provide a glimpse into NZ's future if agricultural intensification continues. The Waingongoro River catchment (WA3) is covered almost entirely by HG (91.2 %), with practically all of this land being used for intensive strip grazing. The SUD da was 15.0 SU ha −1 in 1990 and increased to 15.4 SU ha −1 by 2012. The D HG from 2000 to 2013 had a strong increasing trend of 9.8 % yr −1 RSKSE, associated with the intensification of dairy operations (Wilcock et al., 2009). The result of all this intensification was that WA3 had meaningful increases in TP and DRP. The reason TN and NO x did not display significant trends here is because of the extreme monthly variability in river nitrogen concentrations, possibly due to livestock rotations, fertilizer applications, and precipitation events. Noteworthy is that these significant trends of increasing SUD da , D HG , and nutrients are occurring not only in lowland catchments on the North Island (WA3, HV2) but also in upland catchments of the North Island (RO6), as well as both lowland (TK1) and upland (CH3, TK2) catchments on the South Island.
While disturbance was not itself a strong predictor of water quality, it did help explain outliers of land use-water quality relationships. For example, streams with high DRP (> 20 mg m −3 ; 10th percentile) had one of two dominant land uses, either PF (RO2, RO3) or HG (HM5, WA3, WA9, HM4, HM2). The one exception was RO4, which had relatively low coverage of PF (11.2 %) and HG (2.9 %). In fact, RO4 is dominated by NF (79.1 %). Upon closer examination, we found that the small areas of PF and HG in RO4 were disturbed frequently. Further, most of the disturbed forestry occurred on steep slopes and most of the disturbed pastures (practically all sheep and beef) occurred on hilly terrain adjacent to stream channels. Our high temporal-resolution analyses of disturbance showed that even though this catchment is mostly indigenous forest, intense disturbances on small proportions of developed land can have a considerable impact on water quality. RO4 is also experiencing significant increases in TURB and TP, as well as a significant decrease in Q. Another outlier example was RO3, which was the only non-HGdominated catchment with high NO x (634 mg m −3 ). RO3 was dominated by PF (69.8 %), but it had the highest median disturbance (10.5 %) of all catchments. This catchment also exceeds ANZECC guidelines for DRP and has experienced meaningful increases in TURB, TN, and NO x .
We believe that land disturbance and consequently river water quality will continue to worsen in some NZ catchments based on the following. More plantation forests were planted 1993-1997 (3810 km 2 ) than any other 5-year period in NZ history (NZFFA, 2014). With a 28-year mean age of harvest, NZ will experience its greatest coverage and intensity of forest disturbance around 2025. When combined with drought and intense storms, the potential for nutrient and sediment mobilization is high, especially given that approximately 45 % of these plantings occurred on high-producing grasslands (NZFFA, 2014) where many of the legacy nutrients will be exported to rivers during forest harvest (Davis, 2014). If carbon prices continue to stay low, there will be a high likelihood that many of the harvested forests will be converted to pasture, adding even more nutrients to NZ rivers (PCE, 2013). Given that the central government created a national policy goal of nearly doubling the export to gross domestic product ratio by the year 2025 (MBIE, 2015), NZ is likely to see continued increases in livestock density, fertilizer inputs, and supplemental feed to support these extra livestock, all of which will add even more pressure and risks of eutrophication on NZ's rivers.

Conclusions
This study had the overall goal of describing how changes in land use intensity impact river water quality across broad scales and over long periods. To address this goal we used a combination of "brute force" statistical analyses (in terms of hundreds of analyses using a suite of physiographic, land use, and water quality data for 77 catchments over 26 years) and careful examination (using multi-resolution data to find patterns and relationships among these variables). This goal was ambitious and we likely missed some relationships and details of water quality changes. However, we found empirical evidence for several key relationships among land use intensity, geomorphic processes, and water quality, which we now place into a broader perspective.
The greatest negative impact on river water quality in NZ in recent decades has been high-producing pastures that require large amounts of fertilizer to support high densities of livestock. While this finding has been previously published (Davies-Colley, 2013;Howard-Williams et al., 2010;and references within), our results and supporting information show that the relationship between high-producing pastures and water quality is complicated, being dependent on livestock type/density, disturbance regime, and physiography, particularly soil type. Dairy cattle receive much of the blame for degraded water quality because of their high nutrient requirements (Howard-Williams et al., 2010), but beef cattle can also strongly degrade water quality due to comparable required inputs and grazing on steeper land with a higher potential for runoff (McDowell et al., 2008). Further, pasture designations/boundaries are becoming increasingly blurred by modern cattle management, with greater movements of dairy and beef cattle among pastures, greater use of highproducing pastures for beef, over-wintering of dairy cattle on beef pastures, and cross-breeding (Morris, 2013). While riparian fencing has plausibly improved the clarity of NZ rivers, the removal of millions of sheep from steep slopes has also likely played a role that should be investigated further.
New Zealand is the global leading exporter of whole milk powder, butter, and sheep products, and NZ's prominence in these industries is likely to continue over the next decade (OECD/FAO, 2015). In this most recent environmental review by the Organisation for Economic Co-operation and Development, NZ, had the highest percent increase (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) in agricultural production out of 29 OECD countries, the highest percent increase in N fertilizer use, and the second highest increase in P fertilizer use. This agricultural intensification over our study period is reflected in overall nutrient enrichment of NZ rivers. If cattle continue to be added at the rates we documented, additional fertilizers and supplemental feed will be needed. Even if best management practices are adopted to reduce nutrient export to rivers, there is already a half-century legacy of nutrients distributed across the NZ landscape that will continue to leak to the rivers (Larned et al., 2016). Indeed, the full impact of agricultural intensification on river water quality will not be fully appreciated for another several decades (Howard-Williams et al., 2010;Vant and Smith, 2004). Having an extensive national network like the NRWQN to document and study these water quality changes will be important.

Data availability
The land disturbance data characterized in Figs. 3 and 4 are available on the following website: http://tethys.dges.ou.edu/NZ_disturbance/# (de Beurs, 2017). Users can use a web map to identify their area of interest and view/download the disturbance time series for that particular area. The NRWQN dataset is detailed in Davies-Colley et al. (2011), hosted by NIWA, and is publicly available for non-commercial use at https://niwa. co.nz/freshwater/water-quality-monitoring-and-advice/ national-river-water-quality-network-nrwqn. We note that a large number of NRWQN site operations (at time of writing) are being transferred to regional government authorities, with NIWA monitoring continuing at only about 30 reference sites. Data summaries for all 77 catchments can be found in the Supplement Tables.
The Supplement related to this article is available online at doi:10.5194/hess-21-1149-2017-supplement.