Management of groundwater quality is assisted by an understanding of reference conditions, which describe the concentration ranges expected for key substances in the absence of human impact. This study evaluates reference conditions for NO₃–N in New Zealand groundwater based on three complementary methods: hierarchical cluster analysis, relationships to groundwater age, and regression against a measure of land-use impact. The three methods result in very similar national-scale estimates of reference conditions for NO₃–N concentration in oxic, minimally impacted groundwater, with the 80th, 90th and 95th percentiles found to be 1.65 ± 0.12, 1.97 ± 0.14 and 2.32 ± 0.14 mg/l, respectively (weighted average ± 95% confidence level), in good general agreement with previous studies from New Zealand and overseas. Anoxic groundwaters were treated separately for definition of reference conditions, with the 80th and 90th percentiles of NO₃–N found to be 0.04 ± 0.01 and 0.16 ± 0.01, respectively (the 95th percentile could not be estimated reliably). For both oxic and anoxic groundwater, where a site-specific investigation has not been conducted to estimate reference conditions at a local scale, we suggest that the 80th percentile is an appropriate national-scale default threshold, to match the thresholds used for surface waters under the Australian and New Zealand Guidelines for Fresh and Marine Water Quality.

Introduction

Elevated nitrate-nitrogen (NO₃–N) concentration in groundwater is a worldwide issue (Abascal et al. 2022). Excess NO₃–N in groundwater is derived from a variety of current and legacy sources including agriculture, wastewater, stormwater and industry (Wang et al. 2013; Abascal et al. 2022). Elevated NO₃–N concentrations can affect groundwater ecosystems, compromise the use of groundwater as a potable water supply and, through groundwater-surface water interactions, affect the quality of surface water and dependent ecosystems in receiving environments (Carpenter et al. 1998; Marmonier et al. 2018).

In New Zealand, NO₃–N inputs into freshwater systems are primarily from livestock urine (Parfitt et al. 2012). Urine patches contain 600–1000 kg N ha⁻¹(Monaghan et al. 2007). These loads are well above the nitrogen requirements for pasture growth (cf. 200 kg N ha⁻¹), which leaves the remainder susceptible to leaching (Di and Cameron 2007). Elevated NO₃–N concentrations and increasing temporal trends have been reported for several long-term groundwater monitoring sites across New Zealand (Daughney and Wall 2007; Daughney and Randall 2009; Morgenstern and Daughney 2012; Moreau et al. 2016; StatsNZ 2020), a pattern also evident in surface waters (StatsNZ 2022).

Understanding of reference conditions for NO₃–N in groundwater is vital for selecting appropriate management approaches and objectives (Shand et al. 2007). Reference conditions specify the expected range of NO₃–N concentrations in groundwater in the absence of human impacts. The terms baseline and background are also used (Edmunds et al. 2003; Reimann et al. 2005), but here we use the term reference conditions due to its common use in freshwater policy (Hess et al. 2020). Understanding of reference conditions enables partitioning of observed NO₃–N concentration into the fraction that is anthropogenic vs. natural. Without such understanding, it is difficult to quantify the magnitude of human impacts on groundwater quality (Brydie et al. 2014). In practice it can be scientifically challenging to determine reference conditions representative of the complete absence of human impacts and so some studies instead develop estimates for minimally disturbed conditions (Huo et al. 2012; McDowell et al. 2013; Swanson et al. 2020), a convention also taken in the present study. Regardless of whether defined as a complete absence or a minimal level of human disturbance, understanding of reference conditions enables policy-makers, land managers and communities to set limits or targets for NO₃–N concentrations in groundwater that are realistic and achievable.

Reference conditions are described as a distribution, not a single number, to encapsulate natural range and variation (Reimann et al. 2005; Shand et al. 2007; Edmunds and Shand 2008). Statistically derived numeric thresholds are often selected from the distribution to assist comparison to concentrations that have been measured in water quality investigations. For example, the Australian and New Zealand Guidelines for Fresh and Marine Water Quality define Default Guideline Values for many water quality parameters in surface water, which are based on the 20th and 80th percentiles of the distribution of concentrations inferred to occur under reference conditions (ANZG 2018). Other thresholds have also been recommended, including the 95th percentile (Edmunds et al. 1987) or the 97.7th percentile (Langmuir 1997). The 90th percentile has also been proposed as an appropriate threshold for small datasets (fewer than 60 measurements) or where human impacts cannot be unequivocally excluded (Hinsby et al. 2008). Note that reference condition thresholds are not necessarily management limits or targets that needs to be met, but rather their exceedance may simply be taken as a prompt to consider whether further investigation should be undertaken to determine whether aquatic ecosystems are sufficiently protected.

Reference conditions are often defined separately for specific water bodies or types of water bodies. For example, the range of natural concentrations of NO₃–N in groundwater may vary according to aquifer lithology, redox conditions, flow-paths (including aquifer depth) and groundwater age. Unless there is a geological source of NO₃–N, its concentration under reference conditions is likely to be low under anoxic (reducing) conditions where nitrate is removed through microbiologically-facilitated denitrification, especially in older groundwater. In contrast, detectable levels of NO₃–N are expected for oxic groundwaters, even in the absence of human impacts (Madison and Brunett 1985; Daughney and Reeves 2005). Accordingly, different definitions of reference conditions (and their corresponding percentile-based thresholds) may be needed for different aquifers, parts of aquifers, or aquifer types (Shand et al. 2007).

A variety of methods are available to estimate reference conditions for freshwater quality. Ideally, reference conditions could be determined from measurements made at monitoring sites located in pristine, unimpacted locations (Shand et al. 2007; ANZG 2018). However, this approach is often hampered by a lack of monitoring sites in pristine areas (Dodds and Oakes 2004). This is the case in New Zealand (Daughney et al. 2012; McDowell et al. 2013), especially for groundwater.

In the absence of monitoring sites in pristine locations, another option is to use historical data from any or all available sites regardless of their location, but to estimate reference conditions based only on samples collected in pre-industrial time periods (Edmunds et al. 2003; Limbrick 2003; Griffoen et al. 2008; Hinsby et al. 2008). However, use of historical data may not be possible because older measurements may be of poor quality or not available at all (Shand et al. 2007). As an alternative to reliance on historical measurements, water dating techniques can be used to estimate groundwater quality prior to anthropogenic activity (Morgenstern and Daughney 2012; Huang et al. 2012). Isotope ratios such as ¹⁸O/¹⁶O and ¹⁵N/¹⁴N can also be used with the aid of modelling to estimate the source of NO₃–N and a mass balance of N, allowing the estimation of reference conditions for aquifers with mixed inputs (Minet et al. 2017).

In the absence of data from unimpacted monitoring sites or time periods, various statistical techniques have been used to estimate reference conditions for the chemical composition of groundwater. In early work, reference conditions were determined as those concentrations lower than two standard deviations from the mean of all measured concentrations (Sinclair 1974; Nolan and Hitt 2003). As the mean was vulnerable to skewed data, subsequent refinement changed this to lower than two median absolute deviations from the median of all observed concentrations (Reimann et al. 2005). More recent work has focused on cumulative probability plots (Reimann et al. 2005; Panno et al. 2006; Koh et al. 2009). For example, using log-scale cumulative probability plots, changes in linear slopes are inferred to differentiate NO₃–N concentrations that are indicative of reference conditions versus human impact, or to indicate different processes controlling groundwater quality (Coetsiers et al. 2009; Gemitzi 2012; Rahman et al. 2020). A related approach developed through the EU BRIDGE programme uses cumulative probability plots in concert with pre-selection methods to identify and exclude monitoring results that may show evidence of human impact, for example, based on parameters such as Cl concentration (Griffoen et al. 2008; Hinsby et al. 2008). Iterative outlier removal and gaussian mixture models can also be used to distinguish groundwater quality measurements that represent unimpacted versus impacted conditions (Nakić et al. 2007; Kim et al. 2015; Manu et al. 2022).

For streams, rivers and lakes, the reference conditions have also been inferred by isolating an anthropogenic factor such as the proportion of catchment area under intensive agriculture (Dodds and Oakes 2004; Abell et al. 2019). In surface waters, examples exist where variation is captured by classification systems (Snelder and Biggs 2002), resulting in several definitions of reference conditions. By grouping sites with similar characteristics and plotting contaminant concentrations at sites against, for example, upstream area in intensive agriculture, a regression can be fitted for each environmental class. Here the intercept represents reference conditions, because in this case it is predicting the chemical composition of freshwater in the absence of any upstream intensive agriculture (McDowell et al. 2013).

Each of the above-listed methods for evaluating reference conditions comes with advantages and disadvantages. Statistical methods for estimating reference conditions are fast, cheap, and can be relatively simple, but often suffer from higher uncertainty – especially if based on few measurements over large areas or complex geologies (Cruz and Andrade 2015). Moreover, these statistical methods do not necessarily provide unequivocal insight into the specific geochemical processes or anthropogenic drivers that are controlling the concentrations of NO₃–N. In comparison, the use of isotopes and mass balance modelling is relatively expensive and time consuming but can produce robust definitions for reference conditions. Techniques that rely on environmental classifications and/or comparison of observed NO₃–N concentrations to the levels of anthropogenic stressors, such as area of intensive agriculture, require information on the source areas (capture zones) of recharge, but this information is often lacking for groundwater monitoring sites. Of note for the context of the present investigation, relatively few previous studies have applied more than one method for evaluation of reference conditions for the chemical composition of groundwater, meaning that the strengths and weaknesses of the different techniques are difficult to compare and assess robustly.

The objective of this study is to estimate reference conditions and corresponding thresholds for NO₃–N concentrations in New Zealand groundwater. This is undertaken using national-scale datasets and comparison of results from three methods. The first method is hierarchical cluster analysis (HCA), a multivariate statistical method applied by Daughney and Reeves (2005). The second method involves a comparison of NO₃–N concentrations to groundwater residence times derived by fitting lumped parameter models to measurements of age tracers such as tritium, as demonstrated by Morgenstern and Daughney (2012). The final method employs a comparison of NO₃–N concentrations to a metric of land use intensity in the area around each monitoring site, like the approach that has been applied to determine reference conditions for chemical indicators in New Zealand surface waters (McDowell et al. 2013). While all three methods have been applied previously in New Zealand, they have not employed data from the same sites and time periods, nor have they used consistent percentile-based thresholds that enable meaningful comparison of the inferred reference conditions.

Methods

Calculations were performed with R version 4.2.2, using the specific functions and packages detailed hereafter.

Groundwater quality data

This study used two groundwater quality datasets (). The first dataset was sourced from the New Zealand National Groundwater Monitoring programme (NGMP), a long-term research and monitoring programme that includes approximately 100 monitoring sites located across the country (Supplementary Material Table 1). NGMP sites are selected to encompass a range of land uses, aquifer confinement and lithology. Further details on the NGMP and its sites are provided by Daughney and Reeves (2005) and Daughney et al. (2012). The second dataset included measurements from 944 state-of-the-environment sites monitored by regional authorities and compiled by Land Air Water Aotearoa (LAWA). These state-of-the-environment monitoring sites are selected for a range of purposes according to the different monitoring objectives and network designs devised by the individual regional authorities.

Figure 1. Locations of monitoring sites in the NGMP dataset (left panel) and LAWA dataset (right panel). Symbology corresponds to hydrochemical facies identified by HCA.

Groundwater sampling and analytical methods are generally comparable between the NGMP dataset and the LAWA dataset (Daughney et al. 2012). Most sites are sampled quarterly, though some regional authority sites are monitored monthly or annually. Samples are collected according to a standard protocol (New Zealand Ministry for the Environment 2006). NGMP samples are analysed in the field for electrical conductivity, pH and temperature, and analysed in the laboratory for the concentrations of major ions (Na, K, Ca, Mg, HCO₃, Cl, SO₄) along with Br, F, Fe, Mn, SiO₂ and selected forms of nitrogen and phosphorus (typically at least NO₃–N, NH₄–N and PO₄–P). While many of these same parameters are also monitored by regional authorities (see Daughney and Randall 2009), only NO₃–N, NH₄–N, PO₄–P, Cl and electrical conductivity are routinely made publicly available through the LAWA website (see Data Availability).

For the NGMP and LAWA datasets, this study used measurements made during two non-overlapping ten-year time periods: 1 January 2000 to 31 December 2009; and 1 January 2010 to 31 December 2019. Ten-year time periods were selected to provide sufficient data for analysis within each time window, particularly for estimation of upper and lower percentiles in the distributions of NO₃–N concentration (Helsel et al. 2020). The start and end dates for these time windows are arbitrary but were selected such that the second period used in this study approximated the most recent intervals used for groundwater quality trend reporting by Statistics New Zealand (1999–2018; StatsNZ 2020) and LAWA (2012–2021; LAWA 2022).

For both time periods, the NGMP and the LAWA datasets contain censored analytical results (i.e. results reported as being below a detection limit), which is typical for water quality datasets (Helsel et al. 2020). The level of censoring varies by site, by parameter and over time. For example, censored results for Ca concentrations are very uncommon in the NGMP and LAWA datasets, whereas censored results account for a total of 20% and 11% of all NO₃–N results in these two datasets, respectively, reflecting the very low concentrations that occur in some groundwaters (Langmuir 1997). Other parameters including F, Fe, Mn, NH₄–N and PO₄–P also commonly display censored results for a non-trivial proportion of sites and samples. For these parameters, the censoring threshold can change over time, e.g. due to improved analytical methods. To illustrate, across the NGMP and LAWA datasets there are more than ten different censoring thresholds for NO₃–N, ranging from <0.001 to <1. The NGMP and LAWA datasets were therefore prepared for subsequent analysis using two steps appropriate for censored values as described below.

The first step in preparing the groundwater quality datasets for subsequent analysis was to calculate a median value for each parameter at each site, based on its available time-series data. This was undertaken separately for the two ten-year periods used in this study by applying the ros function of the NADApackage to perform regression on order statistics. This approach uses log-transformation and Weibull plotting positions of uncensored values to replace the censored values with numeric estimates that are then used in the calculation of the site-specific median for a given parameter (Helsel et al. 2020). This approach is suitable for datasets with multiple censoring thresholds and censoring levels of up to 70%–80% (Helsel and Cohn 1988). In this study, the ros function was not applied to any site/parameter dataset that had a censoring level of greater than 70%, and instead the median was recorded as a censored value at the highest censoring threshold for the relevant parameter and site. As a hypothetical example, a site having measured concentrations of <0.005, <0.003 and <0.001 in three different samples would have its median recorded as <0.005. At the end of this data preparation step, the time-series datasets had been reduced to two s × p arrays (one array for the medians calculated for each ten-year period), where s and p were the number of sites and number of parameters in the dataset, respectively, and in which a portion of the results remained as censored values.

The second step in the data preparation was to replace the censored median values that remained after the first step. This was undertaken on a per-parameter basis whereby, across all sites, any remaining censored values were replaced with the highest censoring threshold for that parameter, and all uncensored values less than the same censoring threshold were similarly replaced (Helsel et al. 2020). For example, sites having median values of 0.22, 0.108, <0.03, 0.019 and <0.01 for a particular parameter would be replaced by 0.22, 0.108, 0.03, 0.03 and 0.03. The reason for this approach is that several of the methods used to estimate reference conditions cannot function with censored values, and it is not possible to determine which of these last three values is largest or smallest, hence they were treated as equivalent.

The impacts of these approaches for dealing with censored data are discussed below. Note that while the above method description refers to calculation of median values on a site/parameter basis, the same methods were also used to calculate the 80th percentiles used as input for estimation of reference conditions via regression against a measure of land-use impact.

Results and discussion

Hierarchical cluster analysis

HCA conducted with the NGMP dataset identified three hydrochemical facies of relevance to the present study (). At the highest HCA separation threshold, the NGMP sites were partitioned into Cluster 1 and Cluster 2, representing oxic and anoxic groundwaters, respectively. Groundwater sites assigned to Cluster 1 were typified by concentrations of NO₃–N that were above the analytical detection-limit and concentrations of NH₄–N, Fe and Mn that were near or below their analytical detection limits, whereas sites assigned to Cluster 2 displayed the opposite pattern. The shifts in the concentrations of these redox-sensitive substances reflect their expected behaviour in response to aquifer microbial processes (Chapelle 2000) and natural water-rock interaction (Langmuir 1997). At a lower HCA separation threshold, the oxic groundwaters were subdivided into Cluster 1A and Cluster 1B, interpreted to represent groundwaters that were impacted or minimally impacted by human activity, respectively. We infer that the degree of human impact is a primary driver that differentiates Cluster 1A from Cluster 1B, as evidenced by the former cluster’s typically higher concentrations of NO₃–N, sometimes also accompanied by elevated concentrations of other substances such as K, Na and/or Cl. Further lowering the HCA separation threshold allows additional subclusters to be identified, which within Cluster 1B are inferred to be driven by aquifer lithology, i.e. clastic or carbonate versus volcanic or volcaniclastic (Daughney and Reeves 2005); however, these subclusters are not strongly differentiated by their NO₃–N concentrations so are not discussed in this paper.

Figure 2. Dendrogram for the NGMP dataset based on HCA conducted with site-specific median concentrations of Br, Ca, Cl, F, Fe, HCO₃, K, Mg, Mn, Na, NH₄–N, NO₃–N, PO₄–P, SiO₂ and SO₄ for samples collected in the period 2000–2009 (top panel) or 2010–2019 (bottom panel). Terminus of each vertical line along x-axis represents one NGMP site. Labels correspond to hydrochemical facies described in the text and number of sites (n) assigned to them.

These three hydrochemical facies represented by Clusters 1A, 1B and 2 in the NGMP dataset are very similar to those previously identified using the same HCA methods by Daughney and Reeves (2005) and later by Daughney et al. (2012). Repeated identification of these same main hydrochemical facies over different time periods likely occurs because groundwater chemistry is relatively static or otherwise changing slowly at most NGMP sites (Daughney and Reeves 2006; Moreau and Daughney 2021). For the sites where groundwater quality has changed more significantly over time, it tends to shift from one to another of the three hydrochemical facies described in the preceding paragraph, rather than towards some previously unrecognised hydrochemical facies. This suggests that water-rock interaction, aquifer redox processes and human impacts are overarching and enduring controls on groundwater quality in New Zealand.

Using the three main hydrochemical facies identified by HCA with the NGMP data collected up to the end of 2003, Daughney and Reeves (2005) interpreted that reference conditions for NO₃–N in oxic, minimally impacted groundwater in New Zealand could be defined based on the sites assigned to Cluster 1B. Noting that reference conditions must be defined as a range to reflect natural variation, Daughney and Reeves (2005) proposed that NO₃–N concentrations above the 75th percentile or 95th percentile observed for Cluster 1B were ‘probably’ and ‘almost certainly’ indicative of human impact, respectively, for which the corresponding NO₃–N thresholds were 1.6 and 3.5 mg/l (the 80th percentile was 1.81 mg/l, though this particular threshold was not originally reported). Data from NGMP sites that were assigned to Cluster 1A were inferred to show a significant degree of human impact and so were not used for the estimation of reference conditions for oxic groundwaters. Data from NGMP sites assigned to Cluster 2 were typified by reduced groundwater with low or below-detection concentrations of NO₃–N so were likewise excluded from the estimation of reference conditions for oxic groundwater. However, noting that reduced groundwaters can occur naturally, thresholds were reported separately for NO₃–N in anoxic groundwater (Cluster 2), namely as 0.02 and 0.05 mg/l for the 75th and 95th percentiles, respectively.

Following the same approach using the NGMP dataset in the present study and the 80th percentile as an example, we calculated the corresponding NO₃–N concentration in oxic minimally impacted groundwater (Cluster 1B) to be 1.70 ± 0.23 mg/l if based on samples collected in the decade 2000–2009, or 1.44 ± 0.24 mg/l if based on samples collected in the decade 2010–2019 (see for other selected percentiles). For comparison, the 80th and 90th percentiles for NO₃–N concentration in anoxic groundwater (Cluster 2) were found to be 0.04 ± 0.01 and 0.18 ± 0.02 mg/l (2000–2009 period), or 0.03 ± 0.01 and 0.13 ± 0.02 mg/l (2010–2019 period), respectively (95th percentiles could not be estimated robustly). That the selected percentiles in NO₃–N concentration for each hydrochemical facies changes relatively little between the two ten-year periods used in the present study, or in comparison to the earlier study of Daughney and Reeves (2005), is to be expected if these facies do in fact represent minimally impacted groundwaters, which would therefore not be affected by changing loadings of NO₃–N arising from temporal changes in use of land or water in the area around the monitoring sites.

 

Table 1. Selected percentiles ± 95% confidence intervals for NO₃–N concentration (mg/l) in oxic, minimally impacted groundwater as determined in this study using data from two time periods for the NGMP and LAWA datasets.

The same HCA methods cannot be applied to the LAWA dataset because it only includes analyses of NO₃–N, NH₄–N, PO₄–P, Cl, and electrical conductivity. Therefore, we used the LDA model to predict each site’s hydrochemical facies based only on these available parameters. The LDA model was initially calibrated using the NGMP dataset and was able to correctly apportion 90% of the NGMP sites to their originally assigned clusters (1A, 1B, 2). Next, we confirmed that the LDA model was able to achieve >80% accuracy in predicting the hydrochemical facies assignments previously reported by Daughney and Randall (2009) for a national dataset of over 1000 regional council groundwater monitoring sites. Using the validated LDA model applied to the LAWA sites, we calculated that the 80th percentile for NO₃–N concentration in oxic, minimally impacted groundwater (Cluster 1B, n = 278) is 1.65 ± 0.28 mg/l if based on samples collected in the decade 2000–2009, and 1.68 ± 0.31 mg/l if based on samples collected in the decade 2010–2019 (). These results are in good agreement with the 80th percentile values determined from the NGMP dataset for these same time periods. For anoxic groundwater (Cluster 2) the 80th percentile NO₃–N concentration was found to be 0.03 ± 0.01 mg/l for the period 2000–2009 and 0.04 ± 0.01 mg/l for the period 2010–2019, whereas the 90th percentiles were 0.05 ± 0.02 mg/l for both time periods (results not tabulated). We caution that requirement to use an LDA model to predict cluster assignments means that the percentile values derived from the LAWA dataset should not be used beyond providing general validation of the results from the NGMP dataset.

We note certain advantages and limitations of HCA for estimating reference conditions for groundwater quality. One advantage is that the identification of minimally impacted groundwater is based on several parameters (e.g. Na, K, Cl) in addition to NO₃–N. This means that groundwater that is only marginally degraded based on its NO₃–N concentration can still be identified as impacted based on other indicators, thereby allowing for added sensitivity that a univariate method may miss. A key limitation of the HCA method is that it does not provide any definitive information about the actual processes that differentiate the clusters; instead, these drivers can only be inferred, which may be difficult if there are several interacting processes involved. For example, the distinction between Clusters 1A and 1B may be due not only to the degree of human impact, but also to age of the groundwater or the relative proportions of groundwater recharge that a site receives from rainfall passage through the soil zone compared to seepage from rivers (Daughney and Reeves 2005). Some of these limitations are addressed with the method applied in the following section.

Relationship to groundwater age

Groundwater age distribution parameters for the NGMP sites are provided in Supplementary Material Table 1, based on fitting of LPMs to site-specific time-series measurements of several age tracers, mainly tritium, SF₆ and chlorofluorocarbons. For most sites, the mean groundwater age inferred in this study is comparable to that originally reported by Daughney et al. (2010). However, significant differences in the LPM-derived mean groundwater age do occur for some sites, typically because the additional age tracer data collected since 2010 now allow for constraint of more complex age distributions, e.g. binary mixing models due to convergence of old and young water flow paths at the sampled discharge sites. There is also a small minority of NGMP sites at which the groundwater age distribution is interpreted to have shifted over time, for example due to changes in rates of recharge or abstraction that could affect flows of groundwater through some parts of some aquifer systems.

displays the relationships between mean groundwater age and the median NO₃–N concentrations for the NGMP sites. The highest concentrations of NO₃–N are found in sites with young, oxic groundwaters that show evidence of human impact based on HCA (Cluster 1A). Groundwaters that are oxic but minimally impacted (Cluster 1B) typically have lower but detectable concentrations of NO₃–N, even for sites with mean groundwater age >100 years. Regardless of mean age, anoxic groundwaters tend to have very low concentrations of NO₃–N. These general relationships between mean groundwater age and NO₃–N concentration have been previously reported for the NGMP sites (Morgenstern and Daughney 2012).

Figure 3. Median NO₃–N concentration for the periods 2000–2009 (open symbols) and 2010–2019 (closed symbols) versus inferred mean groundwater age for sites in the NGMP dataset. Symbol colours represent hydrochemical facies. Black lines are derived from the constraint line method: solid hyperbola is the fitted constraint line; solid and dashed horizontal lines indicate the asymptote of the constraint line and the 80th percentile for oxic groundwater, respectively. Blue dotted line is derived from the quantile regression method, indicating the 80th quantile for all Cluster 1B groundwaters.

Estimation of reference conditions for NO₃–N via its relationship to mean groundwater age was first reported by Morgenstern and Daughney (2012) using a visual graph fitting approach. A threshold of >0.2 mg/l was identified for the NGMP dataset as separating reference conditions for unimpacted oxic groundwater versus groundwater affected by ‘low intensity land-use’. We note that this threshold was defined based on only three NGMP sites. A second threshold of 2.5 mg/l NO₃–N was identified to differentiate oxic groundwaters impacted by ‘low-intensity’ versus ‘high intensity’ land-use. No thresholds were reported for NO₃–N in anoxic groundwater.

In the present study, we estimated reference conditions for NO₃–N using a constraint line approach (Hao and Wu 2017). The best-fitting constraint line is a hyperbolic equation with a maximum NO₃–N concentration that approaches an asymptote of 2.2 mg/l as mean groundwater age approaches 250 years (), which is towards the upper end of values reported for the NGMP dataset and the effective limit of interpretability of LPMs based on the age tracers used in this study (Morgenstern and Taylor 2009; Morgenstern and Daughney 2012). Once the constraint line had been fitted, we determined selected percentiles for NO₃–N concentration in oxic minimally impacted groundwater using all oxic groundwaters (Clusters 1A and 1B) that plot below the horizontal asymptote of the constraint line. Again, using the 80th percentile as an example, the threshold for oxic groundwater is determined to be 1.79 ± 0.32 mg/l for the 2000–2009 period, and 1.79 ± 0.35 mg/l for the 2010–2019 period (). We note that there is a possible bias in this approach because the calculation includes some sites that are inferred to show evidence of human impact based on HCA (i.e. sites assigned to Cluster 1A). This constraint line approach cannot be used to estimate percentiles in the NO₃–N concentrations for anoxic groundwater.

For comparison we used the NGMP dataset for a second approach, which ignores the constraint line and instead applies quantile regression using only those sites assigned to Cluster 1B, i.e. oxic and minimally impacted according to HCA. The slope of the quantile regression line is not significantly different from zero (p > 0.05) for either period, so the intercept is taken as the threshold concentration of NO₃–N in oxic minimally impacted groundwater. Again, using the 80th quantile as an example, the regression line’s intercept was determined to be 1.61 ± 0.60 for data from the 2000–2009 period, or 1.61 ± 0.53 for data from the 2010–2019 period (). This approach can also be applied to anoxic groundwater (Cluster 2), which yielded estimates of the 80th and 90th percentiles in NO₃–N concentration of 0.05 ± 0.02 and 0.30 ± 0.07 (2000–2009 period), or 0.05 ± 0.01 and 0.35 ± 0.15 (2010–2019 period), respectively (results not tabulated); 95th percentiles could not be estimated robustly.

Due to the lack of an available compilation of their groundwater age distributions, the LAWA dataset cannot presently be used to validate the results from the NGMP dataset. However, some regional and catchment-scale investigations have been undertaken and reveal generally similar relationships between groundwater quality and groundwater age as observed in the NGMP dataset (e.g. Morgenstern et al. 2015, 2018, 2019). For example, concentrations of parameters such as Na, HCO₃, SiO₂, Fe, Mn and PO₄–P have generally been observed to increase with groundwater age, reflecting likely origin from natural water-rock interaction, whereas concentrations of parameters such as NO₃–N, Cl and SO₄ have generally been observed to be higher in younger groundwaters, suggesting that they may be sourced from human activities and/or subject to reactions that tend to cause their concentrations to decrease over time (e.g. denitrification, sulphate reduction).

We note certain advantages and limitations of estimating reference conditions for groundwater quality based on its relationship to groundwater age. One advantage is that groundwater age can be helpfully compared to the history of land use in the vicinity of the monitoring site, which provides constraint on the level of human impact and therefore the NO₃–N concentration that might be expected. We point out that the LPMs applied in this study follow the assumption that the groundwater age distribution at each site is constant over the long term, despite seasonal fluctuations that may exist at some sites (e.g. Toews et al. 2016). Constancy of groundwater age distribution is a reasonable assumption for New Zealand groundwater systems, which experience temperate climates and are not generally highly exploited, but we acknowledge that this assumption may be violated if this method is applied elsewhere. A disadvantage of the groundwater age approach is that, as a univariate method, the relationships to groundwater quality are explored only one parameter at a time. However, as shown in this study, this limitation can be overcome by categorising the monitoring sites using a multivariate hydrochemical method such as HCA. Another disadvantage of the groundwater age method is that, like HCA, it does not provide any definitive information about the actual land-use intensity that occurred in the vicinity of a well at the time of recharge. This limitation is addressed with the method applied in the following section.

Regression against a measure of land-use impact

displays plots of the 80th percentile in the measured NO₃–N concentration in groundwater at each site compared to the average NO₃–N leaching rate in a circular area around that site (similar plots were generated for other percentiles, not shown). The results were not significantly different for circular area radii from 0.5 to 5.0 km (p > 0.1), so all following results are based on a radius of 2.0 km. We are unable to accurately quantify the uncertainties in the average NO₃–N leaching rates depicted in because such uncertainties were not reported by Ausseil and Manderson (2018). These data were generated using a Nutrient Budget model (Overseer), whose uncertainty based on changing climate and soil data is estimated to be 27 ± 9% (Tavernet and Harper 2022), and we acknowledge that modelled NO₃–N leached from land can be highly variable within and between farms according to how they are managed (Parliamentary Commissioner for the Environment 2018).

Figure 4. 80th percentile in NO₃–N concentration for the period 2010–2019 versus average modelled NO₃–N leached from livestock in 2017 within a 2 km radius of each site (note differences in axis scales), with corresponding linear regression lines and 95% confidence limits (shaded areas): A, All NGMP sites plotted together as a single group. B, NGMP sites segregated by hydrochemical cluster with individual regression lines; symbol size is scaled to the fraction of groundwater that has age less than 5 years. C, All LAWA sites plotted together as a single group. D, LAWA sites segregated by hydrochemical clusters predicted by LDA.

Firstly, we considered all the NGMP sites together in a single group (A). Linear regression produces a slope that is significantly different from zero (p < 0.01). The regression line has an intercept of 1.42 ± 0.53 mg/l, which is interpreted to be the 80th percentile for NO₃–N concentration in groundwater in the absence of any NO₃–N leached from livestock ().

We performed a second evaluation in which the NGMP sites were segregated into hydrochemical facies, and each site’s weighting in the regression was scaled according to its fraction of groundwater aged 5 years or less (B). This grouping of the NGMP sites according to their hydrochemical facies helps to understand any observed relationships between NO₃–N leaching rate and the observed NO₃–N concentrations in groundwater, and the age-weighting maximises the regression influence for those sites with residence time distributions most closely corresponding to the 2017 year used for the NO₃–N leaching map. Note that sites assigned to Cluster 2 are anoxic and hence denitrification is likely to have removed almost all NO₃–N from the groundwater, which explains the lack of discernible relationship to mapped NO₃–N leached from livestock. Sites assigned to Clusters 1A and 1B are oxic, so are more likely to show relationships between NO₃–N concentration and the mapped NO₃–N leached from livestock. Based on their hydrochemistry, sites assigned to Cluster 1A are inferred to show evidence of human/agricultural impact and so, following the approach of McDowell et al. (2013), the 80th percentile for NO₃–N concentration in minimally impacted groundwater can be inferred from the intercept of the regression line (1.87 ± 1.63 mg/l). By comparison, sites assigned to Cluster 1B are interpreted to have little or no evidence of human impact. The Cluster 1B regression line has an intercept of 0.99 ± 0.45 mg/l and a slope that is not significantly different from zero (p > 0.1). In other words, even though NO₃–N leaching from livestock is mapped in the vicinity of some Cluster 1B sites, the impacts of this leaching are not observable as elevated NO₃–N concentrations in the groundwater. Given that Cluster 1B groundwaters are predominantly oxic and hence would have low likelihood for the occurrence of denitrification, we suggest that a lack of relationship between NO₃–N concentration in the groundwater and mapped NO₃–N leaching from livestock could indicate that these sites are recharged predominantly from river seepage rather than rainfall passage through the soil zone, and/or recharged from an area outside the 2.0 km radius used in this study. In either case, the similarity of the intercepts for the Cluster 1A and Cluster 1B regression lines provides independent but very similar estimates for the 80th percentile for NO₃–N concentration in oxic, minimally impacted groundwater in New Zealand. Finally, the regression line for Cluster 2 has an intercept of 0.05 ± 0.45 mg/l and a slope that is not significantly different from zero (p > 0.1), providing an estimate of the 80th percentile NO₃–N concentration that would be expected for anoxic groundwaters. Similar approaches can be used to derive estimates for other percentiles of NO₃–N concentration in oxic and anoxic groundwater under reference conditions ().

We repeated the analyses described above using the LAWA dataset. Again, using the 80th percentile as an example, linear regression based on all LAWA sites as a single group (C) yielded an intercept of 2.23 ± 0.43 mg/l (), providing a similar estimate as derived from the NGMP dataset for of the 80th percentile for NO₃–N concentration in groundwater in the absence of any NO₃–N leached from livestock. As noted above, including all the LAWA sites in the regression does not account for differences in the main hydrochemical drivers such as redox condition or potential degree of human impact. We therefore developed separate regression lines by categorising the LAWA sites according to their hydrochemical facies as based on the LDA model (D). LPM parameters are not presently compiled for the LAWA sites, so it is not possible to weight each site’s influence in the regressions according to its groundwater age distribution. For the 80th percentile, the resulting regression for Cluster 1A sites yields an intercept of 6.73 ± 0.77 mg/l, which is markedly different from the result obtained with the same method applied to the NGMP sites (1.87 ± 1.63) and indeed from all other methods and datasets employed in this study (). Further, some LAWA sites assigned to Cluster 1B have markedly elevated NO₃–N concentrations. We conclude that the uncertainties of the LDA model (potentially causing some sites to be assigned to the incorrect cluster) and the lack of groundwater age information hampers the utility of this approach when applied to the LAWA dataset.

We identify certain advantages and limitations of estimating reference conditions for groundwater quality based on regression against a measure of land-use impact. One clear advantage is that land-use impacts are quantified for each site, whereas the other methods employed in this study can only infer the extent of land-use impacts from indirect proxies. A challenge is to account for each site’s capture zone, groundwater transit time, and any hydrochemical processes, such as denitrification, any of which could confound the relationship between NO₃–N leaching rates and the observed concentrations of NO₃–N in groundwater. As shown in this study, these challenges can be reasonably overcome through combined use of a circular area surrogate for the capture zone, LPMs to describe groundwater age distribution, and HCA to assign sites to distinct hydrochemical facies.

Conclusions

This study has illustrated that reference conditions for NO₃–N in groundwater are helpfully estimated through the application of a range of complementary methods, datasets, and monitoring time periods. At the national scale in New Zealand, we estimate that the 80th percentile in NO₃–N concentration is 1.65 ± 0.12 for oxic, minimally impacted groundwater and 0.04 ± 0.01 for anoxic groundwater, using data collected through the National Groundwater Monitoring programme over the period 1 January 2000 to 31 December 2019. In keeping with the approach used for New Zealand surface waters (ANZG 2018), we suggest that the 80th percentile should be used as national default threshold for reference conditions, for comparison to the NO₃–N concentrations observed in environmental monitoring programmes, except where site-specific investigations have indicated that different threshold values should be used. While this study has focussed on NO₃–N, the same approaches could be used for estimation of reference conditions for other groundwater quality variables.

Acknowledgements

We thank regional authority groundwater scientists for collection of groundwater samples over many years to generate the data used in this study. Funding for analysis of groundwater samples was provided by regional authorities and the New Zealand Ministry of Business, Innovation and Employment as part of the Strategic Science Investment Infrastructure Funds (National Groundwater Monitoring Programme, contract C05X1701). Some funding for this work was received from the Our Land and Water National Science Challenge (Ministry of Business, Innovation and Employment contract C10X1507). We thank Associate Editor Nick Mortimer and two anonymous reviewers for suggestions that greatly improved this manuscript.

Data availability statement

Groundwater quality data from the NGMP are available from GNS Science at https://ggw.gns.cri.nz. Groundwater age interpretations for the NGMP sites are provided in the Supplementary Material for this publication. Groundwater quality data from the LAWA data are available from www.lawa.org.nz.

Disclosure statement

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

References

• Abascal E, Gómez-Coma L, Ortiz I, Ortiz A. 2022. Global diagnosis of nitrate pollution in groundwater and review of removal technologies. Science of the Total Environment. 810:152233. ISSN 0048-9697; https://doi.org/10.1016/j.scitotenv.2021.152233. 
• Abell JM, Özkundakci D, Hamilton DP, van Dam-Bates P, McDowell RW. 2019. Quantifying the extent of anthropogenic eutrophication of lakes at a national scale in New Zealand. Environmental Science & Technology. 53(16):9439–9452. 
• Anayah FM, Almasri MN. 2009. Trends and occurrences of nitrate in the groundwater of the West Bank, Palestine. Appl Geogr. 29:588–601. 
• Angel JC, Peterson EW. 2015. Nitrates in karst systems: Comparing impacted systems to a relatively unimpacted system. J Geogr and Geol. 7:65–76. 
• ANZG. 2018. Australian and New Zealand guidelines for fresh and marine water quality. Australian and New Zealand Governments and Australian State and Territory Governments. Canberra ACT, Australia. www.waterquality.gov.au/anz-guidelines. 
• Ausseil A-G, Manderson A. 2018. Spatial nitrate leaching extent: an update for environmental reporting. Maanaki Whenua Landcare Research Contract Report LC3245. 
• Back W. 1961. Techniques for mapping of hydrochemical facies. United States Geological Survey Professional Paper 424-D: 380-382. 
• Back W. 1966. Hydrochemical facies and groundwater flow patterns in northern part of Atlantic Coastal Plain. United States Geological Survey Professional Paper 498-A. 
• Brydie J, Jones D, Jones JP, Perkins E, Rock L, Taylor E. 2014. Assessment of baseline groundwater physical and geochemical properties for the quest carbon capture and storage project, Alberta, Canada. Energy Procedia. 63:4010–4018. ISSN 1876-6102; https://doi.org/10.1016/j.egypro.2014.11.431. 
• Carpenter SR, Caraco NF, Correll DL, Howarth RW, Sharpley AN, Smith VH. 1998. Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecological Applications. 8(3):559–568. 
• Chapelle FH. 2000. The significance of microbial processes in hydrogeology and geochemistry. Hydrogeol J. 8:41–46. 
• Cloutier V, Lefebrve R, Therrien R, Savard MM. 2008. Multivariate statistical analysis of geochemical data as indicative of the hydrogeochemical evolution of groundwater in a sedimentary rock aquifer system. J Hydrol. 353:294–313. 
• Coetsiers M, Blaser P, Martens K, Walraevens K. 2009. Natural background levels and threshold values for groundwater in fluvial Pleistocene and Tertiary marine aquifers in Flanders, Belgium. Environmental Geology. 57(5):1155–1168. 
• Costa JL, Massone H, Martínez D, Suero EE, Vidal CM, Bedmar F. 2002. Nitrate contamination of a rural aquifer and accumulation in the unsaturated zone. Agricult Water Manag. 57:33–47. 
• Cruz JV, Andrade C. 2015. Natural background groundwater composition in the Azores archipelago (Portugal): a hydrogeochemical study and threshold value determination. Science of The Total Environment. 520:127–135. 
• Daughney CJ, Morgenstern U, van der Raaij R, Reeves RR. 2010. Discriminant analysis for estimation of groundwater age from hydrochemistry and well construction: application to New Zealand aquifers. Hydrogeol J. 18:417–428. 
• Daughney CJ, Raiber M, Moreau-Fournier M, Morgenstern U, van der Raaij R. 2012. Use of hierarchical cluster analysis to assess the representativeness of a baseline groundwater quality monitoring network: comparison of New Zealand’s national and regional groundwater monitoring programs. Hydrogeol J. 20:185–200. 
• Daughney CJ, Randall M. 2009. National groundwater quality indicators update: state and trends 1995-2008. GNS Science Consultancy Report 2009/145. GNS Science, Wellington, New Zealand. https://environment.govt.nz/publications/national-groundwater-quality-indicators-update-state-and-trends-1995-2008/. 
• Daughney CJ, Reeves RR. 2005. Definition of hydrochemical facies in the New Zealand National Groundwater Monitoring Programme. J Hydrol. (New Zealand). 44:105–130. 
• Daughney CJ, Reeves RR. 2006. Analysis of temporal trends in New Zealand’s groundwater quality based on data from the National Groundwater Monitoring Programme. J Hydrol. (New Zealand). 45(1):41–62. 
• Daughney CJ, Wall M. 2007. Groundwater quality in New Zealand: state and trends 1995-2006. GNS Science Consultancy Report 2007/23. 65p. 
• Davis M. 2014. Nitrogen leaching losses from forests in New Zealand. NZ J Forestry Sci. 44:2. 
• De Stefano L, Martínez-Santos P, Villarroya F, Chico D, Martínez-Cortina L. 2013. Easier said than done? The establishment of baseline groundwater conditions for the implementation of the water framework directive in Spain. Water Resour Manage. 27:2691–2707. doi:10.1007/s11269-013-0311-6. 
• Di HJ, Cameron KC. 2007. Nitrate leaching losses and pasture yields as affected by different rates of animal urine nitrogen returns and application of a nitrification inhibitor—a lysimeter study. Nutrient Cycling in Agroecosystems. 79(3):281–290. 
• Dodds WK, Oakes RM. 2004. A technique for establishing reference nutrient concentrations across watersheds affected by humans. Limnology and Oceanography: Methods. 2(10):333–341. 
• Dollery R, Li S, Dickenson NM. 2019. Nutrient-enriched soils and native N-fixing plants in New Zealand. J Plant Nutr Soil Sci. 182:104–110. 
• Edmunds WM, Cook JM, Darling WG, Kinniburgh DG, Miles DL, Bath AH, Morgan-Jones M, Andrews JN. 1987. Baseline geochemical conditions in the chalk aquifer, Berkshire, U.K.: a basis for groundwater quality management. Appl Geochem. 2:251–274. 
• Edmunds WM, Ma J, Aeschbach-Hertig W, Kipfer R, Darbyshire DPF. 2006. Groundwater recharge history and hydrogeochemical evolution in the Minqin Basin, North West China. Appl Geochem. 21:2148–2170. 
• Edmunds WM, Shand P, editors. 2008. Natural groundwater quality. Oxford, UK: Blackwell Publishing. 
• Edmunds WM, Shand P, Hart P, Ward RS. 2003. The natural (baseline) quality of groundwater: a UK pilot study. Sci Tot Env. 310:23–35. 
• Farnham IM, Stetzenbach KJ, Singh AK, Johannesson KH. 2002. Deciphering groundwater flow systems in Oasis Valley, Nevada, using trace element geochemistry, multivariate statistics, and geographical information system. Chemometrics and Intellegent Laboratory Systems. 60:265–281. 
• Fleming SW, Quilty EJ. 2006. Aquifer responses to El Niño–Southern Oscillation. Southwest British Columbia. Groundwater. 44:595–599. 
• Freeze RA, Cherry JA. 1979. Groundwater. New Jersey: Prentice Hall. 
• Gemitzi A. 2012. Evaluating the anthropogenic impacts on groundwaters; a methodology based on the determination of natural background levels and threshold values. Environmental Earth Sciences. 67(8):2223–2237. 
• Griffoen J, Passier HF, Klein J. 2008. Comparison of selection methods to deduce natural background levels for groundwater units. Env Sci Tech. 42:4863–4869. 
• Guggenmos MR, Daughney CJ, Jackson BM, Morgenstern U. 2011. Regional-scale identification of groundwater-surface water interaction using hydrochemistry and multivariate statistical methods, Wairarapa Valley, New Zealand. Hydrol Earth Syst Sci. 15:3383–3398. 
• Güler C, Thyne GD, McCray JE, Turner AK. 2002. Evaluation of graphical and multivariate statistical methods for classification of water chemistry data. Hydrogeol J. 10:455–474. 
• Hao R, Wu J. 2017. Relationship between paired ecosystem services in the grassland and agro-pastoral transitional zone of China using the constraint line method. Agriculture, Ecosystems & Environment. 240:171–181. 
• Helsel DR, Cohn TA. 1988. Estimation of descriptive statistics for multiply censored water quality data. Water Resources Research. 24:1997–2004. 
• Helsel DR, Hirsch RM, Ryberg KR, Archfield SA, Gilroy EJ. 2020. Statistical methods in water resources. U.S. Geological Survey Techniques and Methods, book 4, chap. A3, 458 p. https://doi.org/10.3133/tm4a3. 
• Hess S, Alve E, Andersen TJ, Joranger T. 2020. Defining ecological reference conditions in naturally stressed environments – how difficult is it? Marine Environmental Research. 156:104885. 
• Hinsby K, Condesso de Melo MT, Dahl M. 2008. European case studies supporting the derivation of natural background levels and groundwater threshold values for the protection of dependent ecosystems and human health. Sci Tot Environ. 401:1–20. 
• Hope BK, Clarkson JR. 2014. A strategy for using weight-of-evidence methods in ecological risk assessments. Human and Ecol Risk Assessment. 20:290–315. 
• Huang T, Pang Z, Yuan L. 2012. Nitrate in groundwater and the unsaturated zone in (semi)arid northern China: baseline and factors controlling its transport and fate. Env Earth Sci. 70:145–156. 
• Huo SL, Zan FY, Chen Q, Xi BD, Su J, Ji DF, Xu QG. 2012. Determining reference conditions for nutrients, chlorphyll a and Secchi depth in Yungui Plateau ecoregion lakes, China. Water Env J. 26:324–334. 
• Johnson TD, Belitz K. 2009. Assigning land use to supply wells for the statistical characterization of regional groundwater quality: correlating urban land use and VOC occurrence. J Hydrol. 370:100–108. 
• Jurgens BC, Böhlke JK, Eberts SM. 2012. TracerLPM (Version 1): an Excel® workbook for interpreting groundwater age distributions from environmental tracer data. U.S. Geological Survey Techniques and Methods Report 4-F3, 60p. 
• Kim KH, Yun ST, Kim HK, Kim JW. 2015. Determination of natural backgrounds and thresholds of nitrate in South Korean groundwater using model-based statistical approaches. J Geochem Expl. 148:196–205. 
• Koh D-C, Chae G-T, Yoon Y-Y, Kang B-R, Koh G-W, Park K-H. 2009. Baseline geochemical characteristics of groundwater in the mountainous area of Jeju Island, South Korea: implications for degree of mineralization and nitrate contamination. J Hydrol. 376:81–93. 
• Land Air Water Aotearoa. 2022. Groundwater quality. https://www.lawa.org.nz/explore-data/groundwater-quality/#/tb-national. 
• Langmuir D. 1997. Aqueous environmental geochemistry. New Jersey: Prentice-Hall. 
• Limbrick KJ. 2003. Baseline nitrate concentration in groundwater of the chalk in South Dorset, UK. Sci Tot Env. 314-316:89–98. 
• Linkov I, Loney D, Cormier S, Satterstrom FK, Bridges T. 2009. Weight-of-evidence evaluation in environmental assessment: review of qualitative and quantitative approaches. Sci Tot Env. 407:5199–5205. 
• Madison RJ, Brunett JO. 1985. Overview of the occurrence of nitrate in groundwater of the United States. In: National Water Summary 1984 – Hydrologic events, selected water quality trends, and groundwater resources. United States Geological Survey Water-Supply Paper 2275. Washington, DC: United States Government Printing Office; p. 93–105. 
• Maloszewski P, Zuber A. 1996. Lumped parameter models for the interpretation of environmental tracer data. Report IAEA-TECDOC--910. Vienna, Austria: International Atomic Energy Agency. 
• Manu E, Afrifa GY, Ansah-Narh T, Sam R, Loh YSA. 2022. Estimation of natural background and source identification of nitrate-nitrogen in groundwater in parts of the Bono, Ahafo and Bono East regions of Ghana. Groundwater for Sustainable Development. 16:100696. 
• Marmonier P, Maazouzi C, Baran N, Blanchet S, Ritter A, Saplairoles M, Dole-Olivier MJ, Galassi DMP, Eme D, Doledec S, Piscart C. 2018. Ecology-based evaluation of groundwater ecosystems under intensive agriculture: a combination of community analysis and sentinel exposure. Sci Tot Env. 613-614:1353–1366. 
• Martindale H, van der Raaij R, Daughney CJ, Morgenstern U, Singh R, Jha N, Hadfield J. 2019. Assessment of excess N2 for quantifying actual denitrification in New Zealand groundwater systems. J Hydrol. NZ. 58:1–17. 
• McDowell RW, Snelder TH, Cox N, Booker DJ, Wilcock RJ. 2013. Establishment of reference or baseline conditions of chemical indicators in New Zealand streams and rivers relative to present conditions. Marine & Freshwater Research. 64(5):387–400. 
• Minet EP, Goodhue R, Meier-Augenstein W, Kalin RM, Fenton O, Richards KG, Coxon CE. 2017. Combining stable isotopes with contamination indicators: a method for improved investigation of nitrate sources and dynamics in aquifers with mixed nitrogen inputs. Water Research. 124:85–96. 
• Monaghan RM, Hedley MJ, Di HJ, McDowell RW, Cameron KC, Ledgard SF. 2007. Nutrient management in New Zealand pastures— recent developments and future issues. New Zealand Journal of Agricultural Research. 50:181–201. 
• Moreau M. 2021. Long-term trend analysis of groundwater quality at Te Waikoropupū Springs – Tākaka. Wairakei, NZ: GNS Science Consultancy Report 2021/06. 39 p. 
• Moreau M, Daughney CJ. 2021. Defining natural baselines for rates of change in New Zealand’s groundwater quality: dealing with incomplete or disparate datasets, accounting for impacted sites, and merging into state of the-environment reporting. Sci Tot Env. 755:143292. 
• Moreau M, Riedi MA, Aurisch K. 2016. Update of National Groundwater Quality Indicators: state and trends 2005–2014. GNS Science Consultancy Report 2015/107. 42 p. 
• Morgan CO, Winner MD Jr. 1962. Hydrochemical facies in the 400 foot and 600 foot sands of the Baton Rouge area, Louisiana. United States Geological Survey Professional Paper 450-B: B120-121. 
• Morgenstern U. 2008. Lake Taupo catchment groundwater age distribution and implications for future land-use impacts. GNS Science Consultancy Report 2007/301. 38 p. 
• Morgenstern U, Begg JG, van der Raaij R, Moreau M, Martindale H, Daughney CJ, Franzblau RE, Stewart MK, Knowling MJ, Toews MW, et al. 2018. Heretaunga Plains aquifers: groundwater dynamics, source and hydrochemical processes as inferred from age, chemistry, and stable isotope tracer data. GNS Science Report 2017/33. 
• Morgenstern U, Daughney CJ. 2012. Groundwater age for identification of baseline groundwater quality and impacts of land-use intensification – The National Groundwater Monitoring Programme of New Zealand. Journal of Hydrology. 456-457:79–93. 
• Morgenstern U, Daughney CJ, Leonard G, Gordon D, Donath FM, Reeves R. 2015. Using groundwater age and hydrochemistry to understand sources and dynamics of nutrient contamination through the catchment into Lake Rotorua, New Zealand. Hydrol Earth Syst Sci. 19:803–822. 
• Morgenstern U, Davidson P, Townsend D, White P, van der Raaij R, Stewart M, Moreau M, Daughney C. 2019. From rain through river catchment to aquifer: the flow of water through the Wairau hydrologic system. GNS Science Report 2019/63. 
• Morgenstern U, Reeves R, Daughney C, Cameron S, Gordon D. 2004. Groundwater age and chemistry, and future nutrient load for selected Rotorua Lakes catchments. Institute of Geological & Nuclear Sciences science report 2004/31. 73p. 
• Morgenstern U, Taylor CB. 2009. Ultra low-level tritium measurement using electrolytic enrichment and LSC. Isot Environ Health Stud. 45:96–117. 
• Musgrove M, Opsahl SP, Mahler BJ, Herrington C, Sample TL, Banta JR. 2016. Source, variability, and transformation of nitrate in a regional karst aquifer: Edwards aquifer, central Texas. Sci Tot Env. 568:457–469. 
• Nakić Z, Posavec K, Bačani A. 2007. A visual basic spreadsheet macro for geochemical background analysis. Groundwater. 45:642–647. 
• Nejatijahromi Z, Nassery HR, Hosono T, Nakhaei M, Alijani F, Okumura A. 2019. Groundwater nitrate contamination in an area using urban wastewaters for agricultural irrigation under arid climate condition, southeast of Tehran, Iran. Agr Water Management. 221:397–414. 
• New Zealand Ministry for the Environment. 2006. A national protocol for state of the environment groundwater sampling in New Zealand. Report ME781, Ministry for the Environment, Wellington, New Zealand. https://environment.govt.nz/publications/a-national-protocol-for-state-of-the-environment-groundwater-sampling-in-new-zealand/. 
• New Zealand Ministry for the Environment. 2022. National policy statement for freshwater management 2020 – amended December 2022. https://environment.govt.nz/assets/publications/National-Policy-Statement-for-Freshwater-Management-2020.pdf. 
• Nolan BT, Hitt KJ. 2003. Nutrients in shallow ground waters beneath relatively undeveloped areas in the conterminous United States. 21 p. 
• Panno SV, Kelly WR, Martinsek AT, Hackley KC. 2006. Estimating background and threshold nitrate concentrations using probability graphs. Groundwater. 44(5):697–709. 
• Parfitt RL, Stevenson BA, Dymond JR, Schipper LA, Baisden WT, Ballantine DJ. 2012. Nitrogen inputs and outputs for New Zealand from 1990 to 2010 at national and regional scales. NZ J. Agricultural Research. 55:241–262. 
• Parliamentary Commissioner for the Environment. 2018. Overseer and regulatory oversight: models, uncertainty and cleaning up our waterways. Parliamentary Commission for the Environment, Wellington. 
• Pauwels H, Muller D, Griffoen J, Hinsby K, Melo T, Brower R. 2007. BRIDGE Background cRiteria for the IDentification of Groundwater thrEsholds. Final Activity Report. [accessed 2022 Dec]. https://cordis.europa.eu/project/id/6538/reporting. 
• Radu E, Balaet R, Vliegenthart F, Schipper P. 2010. Derivation of threshold values for groundwater in Romania, in order to distinguish point & diffuse pollution from natural background levels. Environ Eng Res. 15:85–91. 
• Rahman A, Tiwari KK, Mondal NC. 2020. Assessment of hydrochemical backgrounds and threshold values of groundwater in a part of desert area, Rajasthan, India. Environmental Pollution. 266:115150. 
• Raiber M, White PA, Daughney CJ, Tschritter C, Davidson P, Bainbridge S. 2012. Three-dimensional geological modelling and multivariate statistical analysis of water chemistry data to analyse and visualise aquifer structure and groundwater composition in the Wairau Plain, Marlborough District, New Zealand. J Hydrol. 436-437:13–34. 
• Reimann C, Filzmoser P, Garrett RG. 2005. Background and threshold: critical comparison of methods of determination. Science of The Total Environment. 346(1):1–16. 
• Seaber GE. 1962. Cation hydrochemical facies of groundwater in the Englishtown Formation, New Jersey. United States Geological Survey Professional Paper 450-B: B124–126. 
• Shand P, Edmunds WM. 2008. The baseline inorganic chemistry of European groundwaters. In: Edmunds WM, Shand P, editors. Natural groundwater quality. New York: John Wiley & Sons; p. 22–58. 
• Shand P, Edmunds WM, Lawrence AR, Smedley PL, Burke S. 2007. The natural (baseline) quality of groundwater in England and Wales. British Geological Survey Research Report No. RR/07/06. 
• Sinclair AJ. 1974. Selection of threshold values in geochemical data using probability graphs. Journal of Geochemical Exploration. 3(2):129–149. 
• Snelder TH, Biggs BJF. 2002. Multi-scale river environment classification for water resources management. Journal of the American Water Resources Association. 38:1225–1239. 
• Snelder TH, Larned ST, Fraser C, De Malmanche S. 2021. Effect of climate variability on water quality trends in New Zealand rivers. Marine and Freshwater Res. 73:20–34. 
• StatsNZ. 2019. Nitrate leaching from livestock – published April 2019. [accessed 2022 April 1]. https://www.stats.govt.nz/indicators/nitrate-leaching-from-livestock. 
• StatsNZ. 2020. Groundwater quality. [accessed 2022 April 1]. https://www.stats.govt.nz/indicators/groundwater-quality. 
• StatsNZ. 2022. River water quality – Nitrogen – published April 2022. [accessed 2022 July 1]. https://www.stats.govt.nz/indicators/river-water-quality-nitrogen. 
• Stevens G. 2010. Groundwater quality in Tasman District. Tasman District Council, Reference R10003. 
• Swanson SK, Graham GE, Hart DJ. 2020. Using reference springs to describe expected flow, temperature, and chemistry conditions for geologically related groups of springs. Env Eng Geosci. 26:331–344. 
• Tavernet J-P, Harper J. 2022. Uncertainty analysis – propagation of input climate data and soil uncertainties on nitrate-loss estimated by the Overseer science model. [accessed 2023 May 22]. https://assets.ctfassets.net/bo1h2c9cbxaf/4HprUvSG1fiB7eeoRJgrxW/40a39b19ac30c3d49dcfdc818c62a436/Uncertainty_analysis_report.pdf. 
• Toews MW, Daughney CJ, Cornaton FJ, Morgenstern U, Evison RD, Jackson BM, Petrus K, Mzila D. 2016. Numerical simulation of transient groundwater age distributions assisting land and water management in the Middle Wairarapa Valley, New Zealand. Water Resour Res. 52:9430–9451. 
• Toews MW, Donath F. 2015. Capture zone delineation of community supply wells and State of the Environment monitoring wells in the Greater Wellington Region. GNS Science Report 2015/06. 
• United States Environmental Protection Agency. 2016. Weight of evidence in ecological assessment. Washington, DC: Risk Assessment Forum. Report no. EPA/100/R-16/001. 
• Urresti-Estala B, Carrasco-Cantos F, Vadillo-Pérez I, Jiménez-Gavilán P. 2013. Determination of background levels on water quality of groundwater bodies: A methodological proposal applied to a Mediterranean River basin (Guadalhorce River, Málaga, southern Spain). Journal of Environmental Management. 117:121–130. 
• Utom AU, Werben U, Leven C, Müller C, Knöller K, Vogt C, Dietrich P. 2020. Groundwater nitrification and denitrification are not always strictly aerobic and anaerobic processes, respectively: an assessment of dual-nitrate isotopic and chemical evidence in a stratified alluvial aquifer. Biogeochemistry. 147:211–223. 
• Wang L, Butcher AS, Stuart ME, Gooddy DC, Bloomfield JP. 2013. The nitrate time bomb: a numerical way to investigate nitrate storage and lag time in the unsaturated zone. Environmental Geochemistry and Health. 35:667–681. 
• Wheeler DC, Nolan BT, Flory AR, DellaVale CT, Ward MH. 2015. Modeling groundwater nitrate concentrations in private wells in Iowa. Sci Tot Env. 536:481–488. 
• Wilson SR, Close ME, Abraham P, Sarris TS, Banasiak L, Stenger R, Hadfield J. 2020. Achieving unbiased predictions of national-scale groundwater redox conditions via data oversampling and statistical learning. Sci Tot Env. 705:135877. 
• Zuber A, Witczak S, Rózański K, Śliwka I, Opoka M, Mochalski P, Kuc T, Karlikowska J, Kania J, Jackowicz-Korczyński M, Duliński M. 2005. Groundwater dating with 3H and SF6 in relation to mixing patterns, transport modelling and hydrochemistry. Hydrol Proc. 19:2247–2275.