The implications of lag times between nitrate leaching losses and riverine loads for water quality policy

Understanding the lag time between land management and impacts on riverine nitrate–nitrogen (N) loads is critical to understand when action to mitigate nitrate–N leaching losses from the soil profile may start improving water quality. These lags occur due to leaching of nitrate–N through the subsurface (soil and groundwater). Actions to mitigate nitrate–N losses have been mandated in New Zealand policy to start showing improvements in water quality within five years. We estimated annual rates of nitrate–N leaching and annual nitrate–N loads for 77 river catchments from 1990 to 2018. Lag times between these losses and riverine loads were determined for 34 catchments but could not be determined in other catchments because they exhibited little change in nitrate–N leaching losses or loads. Lag times varied from 1 to 12 years according to factors like catchment size (Strahler stream order and altitude) and slope. For eight catchments where additional isotope and modelling data were available, the mean transit time for surface water at baseflow to pass through the catchment was on average 2.1 years less than, and never greater than, the mean lag time for nitrate–N, inferring our lag time estimates were robust. The median lag time for nitrate–N across the 34 catchments was 4.5 years, meaning that nearly half of these catchments wouldn’t exhibit decreases in nitrate–N because of practice change within the five years outlined in policy.

. Summary statistics for biophysical parameters describing the 77 catchments of the National River Water Quality Network. Geology was added as a categorical variable assigned to either sedimentary, metamorphic or volcanic 74 . 1 SEM standard error of the mean. 2 Calculated according to the method of Pelletier and Andréassian 55 .  81 The lag time between nitrate-N leaching losses and river loads was calculable for 49 of the Network catchments using at least one of two methods, a cross-correlation function and cumulative generalised additive models (cumulative GAMs). Lags could not be calculated for 28 catchments. However, of the 49 catchments, five were excluded from further analysis as they were deemed as 'impacted' by point sources or had highly modified flows caused by hydroelectricity generation, and 10 were excluded as they were under MDC and therefore unlikely to show any changes in anthropogenic inputs. Although these exclusions may have biased our estimates of median lag times, we have no evidence to say that their distribution of lag times would be any different to those included in further analysis. Indeed, a Mann-Whitney test of the coefficients of variation for nitrate-N leaching losses and loads showed no difference between median coefficient in those catchments where calculations could or could not establish lag times. This suggests that it was not the magnitude of change but when changes were occurring that was causing the detection (or not) of lag times. Supporting this, a Mann-Whitney test of the annual median slope for leaching losses over time showed those catchments with a lag time (1.445 tonnes yr −1 ) were significantly greater than those without lag times (0.0026 tonnes yr −1 ). The near zero increase in mean annual nitrate-N leaching losses is to be expected in catchments with no significant agricultural activity.
Significant lag times in nitrate-N loads were estimated for 34 catchments, 18 and 31 via the cross-correlation and cumulative GAM approaches, respectively; 12 catchments exhibited significant lag times using both approaches ( Table 2, Fig. 3). Lag times using the cross-correlation approach were related to those using the cumulative GAM approach (Cross-correlation = 1.23 × Cumulative GAM − 3.2, R 2 = 0.41, P < 0.019). The mean and median lag time across all 34 catchments were 5.1 and 4.5 years, respectively.
Regression analysis indicated that lag times between nitrate-N leaching losses and loads (of individual techniques or a mean where both returned a significant result) could be predicted (adjusted R 2 = 0.47; Supplementary  Table S2) with increasing mean catchment altitude (metres above sea level), stream order at the catchment outlet (Strahler), and deceasing mean slope (degrees), population density (people km −2 ) and geologic class (going from sedimentary to volcanic): Amongst these parameters, the most variance was explained by mean altitude, mean slope, stream order, and evapotranspiration, with smaller contributions from geologic class, and population density. Including other parameters such as land use class, area, baseflow fraction, and rainfall resulted in small improvements in predictive power and an overfitted model.
Hydrologic comparison of lag times. Mean transit time data for water (or hydrologic lag) were available for 14 catchments, eight of which had significant lag times estimated by either the cross-correlation or cumulative GAM approaches. Mean lag times in these eight catchments were all greater than or within the range of mean transit times for baseflow (Table 2), which is likely to be a combination of young, shallow groundwater and older interflow. On average, nitrate-N lag times in these nine catchments were 2.1 years longer than the midpoint of mean transit times for baseflow.
(1) Lag = −9.1 + 0.0085Mean Altitude + 0.592Stream order −1.097Mean Slope + 0.227Evapotranspiration − 0.068Populaton density − 0.784Geology class Figure 1. Range of the cumulative sum of nitrate-N leached (tonnes) over 1990-2018 for each livestock class, the sum of livestock classes and the load in the river. Also given is the range of nitrate-N loads retained in the river relative to the sum leached. Boxes describe the 25th, 50th and 75th percentiles, whiskers are the 5 and 95th percentiles. Outliers are indicated by black dots.   Table 2. Mean yield (± standard error) of nitrate-N and lag times (years) calculated using cross-correlation (using pre-whitened data) and cumulative generalised additive models (GAMs) and the resulting mean lag time after filtering out sites that were either impacted by hydroelectric schemes or were recently disturbed causing uncharacteristically high sediment loads, or under minimally disturbed conditions (MDC), i.e. > 90% native forest, mountains or scrub. The mean hydrological transit time for baseflow is also given for reference. Lag times in parentheses were for sites that were either impacted or under MDC and excluded from further analysis. 1 A lag time either could not be calculated or no data exists.   Zealand boundary data sourced from https:// www. stats. govt. nz/ while catchment boundaries were sourced from the River Environment Classification database (https:// niwa. co. nz/ fresh water-and-estua ries/ manag ement-tools/ river-envir onment-class ifica tion-0). www.nature.com/scientificreports/ N-retention did not correlate to lag times and was only partly explained by an increase in area and mean altitude, presumably as leached nitrate-N was subjected to a longer and more tortuous route to the catchment outlet leading to more removal by processes like denitrification 20 . The absence of predictors like land use may reflect the fact that we included these as a static point in time, whereas they would be changing during the study. However, as different catchments would be changing at different times, no one metric to capture that change could be included in the model.
A negative N-retention indicates that more nitrate-N was leaving the catchment than was estimated to be leached. Negative N-retention percentages were largely confined to those in MDC. The negative N-retention from these catchments is likely caused by low agricultural activity and not including N inputs via atmospheric deposition or erosion, which in agricultural catchments would be obscured by N inputs from animals 21 . Parfitt et al. 22 estimated losses via erosion and leaching of N in catchments under native forestry in New Zealand, i.e. MDC, to be about 3 kg N ha −1 yr −1 , whereas atmospheric inputs in rainfall were 1.5 kg N ha −1 yr −1 except in catchments with high rainfall where inputs can be > 5 kg N ha −1 yr −1 . On average, MDC catchments had a 350 m greater mean altitude and about a 50% greater rainfall than other catchments (see Supplementary data).
Amongst catchments that were impacted by agricultural activity, the range of N-retention (− 2 to 99%) was broader than found in other studies. As mentioned above, our predictor variables were largely unable to explain variation in N-retention. Other studies have focused on the interaction of catchment characteristics such as biogeochemistry 23 , and catchment management such as point sources, agricultural intensification 17,24 and flowpaths 8,25 . In contrast to these studies, our catchments were more diverse in altitude, geology, and land use intensity (Table 1), and had little point source contributions 26 . For instance, the catchments studied by Dupas et al. 17 had a rainfall range of 195-689 mm yr −1 and agricultural land use not less than 30%; ours ranged from 828-3696 mm and 0-99% (Table 1). Such diversity should generate greater predictive power, if spread linearly over the range of values. However, the large size of our catchments (Table 1) coupled with a large intra-catchment diversity in altitude, soil type, and climate implies that mean catchment variables may not be representative of the catchment processes that control N-retention, or that the interaction of these variables causes a lag in their effect. Indeed, long lag times of up to 34 years (mean of 5.5 years) caused by a diversity of flow paths was labelled as a cause of attenuation of nitrate-N loads in the large Grand River catchment (6800 km 2 ) of southwestern Ontario 12 .

Lag times.
A handful of studies have used a data driven approach to establish lag times between nitrate-N inputs (or leaching losses) and nitrate-N loads in rivers. For instance, Van Meter and Basu 12 found lag times varied from 12-34 years with longer lag times in the lower watershed corresponding to catchments that were dominated by groundwater flows. Similarly, transit time distributions (effectively lag times in the way they were assessed) of 2-14 years were calculated by Dupas et al. 17 with longer lag times associated with catchments with granite and mixed-lithology than calculated for catchments with schist lithology. In a mixed forest/agricultural catchment, Ehrhardt et al. 24 found lag times in nitrate-N ranging from 7-22 years depending on stream location and season as either shallow, proximal (shorter lag times; e.g. through riparian soils and tile drains) or deeper (longer lag times; e.g. aquifer-driven) flow-paths contributed the nitrate load.
Our data showed that a combination of mean altitude, mean slope, evapotranspiration, and stream order strongly influenced lag times, with weaker influence from geologic class and population density. An increase in lag time is expected as stream order increases 27 because deeper flow-paths carrying older water will likely contribute more of the total stream flow. However, in catchments of comparative size but with greater mean altitude and slope, steeper slopes will promote deeper vertical infiltration resulting in a wider range of flow paths with different ages than in flatter areas, lengthening lag times 20,28 . Of the other variables, increasing lag times would be expected with decreasing population density which would act as a surrogate for increasing stream order and catchment size. Likewise, certain geologic classes can increase lag times, especially in catchments with porous bedrock like chalk 29 . Although limestone and chalk geology is rare in New Zealand, long lag times (70-110 years) have been noted for groundwater in catchments with porous volcanic geology 7 .
Nitrate-N lag times are likely longer than hydrologic lag (viz. transit) times, owing to the presence of biogeochemical lags that result in N being stored in the soil 30 . However, data is emerging to suggest that hydrologic lags dominate overall time lags in nitrate loads. This has been attributed to not only the mobile nature of nitrate-N but also to periods of sustained N inputs 24 or of diminished biogeochemical N retention 4 -effectively surpassing the ability of soils and flow paths to remove added N. Dupas et al. 17 found lag times that were equivalent to hydrologic times but noted that a biogeochemical lag was still likely which would lengthen the tail of nitrate-N delivery.
We made no attempt to differentiate between hydrologic and biogeochemical lags. However, owing to the presence of isotope and modelling data for nine of the catchments with nitrate-N lags we could infer that the nitrate-N lag was on average only 2.1 years longer than the hydrologic lag. However, data for hydrologic lag times were derived under baseflow conditions (except for CH4), meaning that younger water from surface runoff pathways would likely decrease the hydrologic lag. The baseflow index varied from 0.03 to 0.88 across all 77 catchments and from 0.17 to 0.50 in the eight catchments that had hydrologic lag times. Assuming water from surface runoff was at least half the age of groundwater (viz. baseflow) 31 , a mean calculated hydrologic lag time for these nine catchments would be 1.1, approximately 3.9 years younger than the mean nitrate-N lag.

Limitations of lag time calculations.
Our lag time calculations have limitations in both data and analytical approaches. Although nitrate-N leaching losses were calculated using software that has been developed and calibrated to nitrate-N losses in New Zealand 32 , the good spatial representativeness that was achieved at a farm scale during census years (varying from 3-6 years apart) reduced to trends in land use associated with a local government district level in other years. However, one benefit of calculating lag times in large catchments is that they often match district boundaries and land use change tends to be within rather than between catch- www.nature.com/scientificreports/ ments. Outside of pastoral land, our estimates of nitrate-N leaching losses may have also been hampered by the accuracy of export coefficients for other crops. However, their coverage and expansion over the last 30 years has been small ( Supplementary Fig. S2).
Analytically, there are also limitations to the utility of the cross-correlation and cumulative GAM approaches to lag times. Considering cross-correlation, one limitation is that correlation does not confirm a causal link between catchment nitrate-N leaching losses and riverine nitrate-N loads 33 . Secondly, the cross-correlation approach depends on the changes in the two series to identify significant lags. Since our catchment N leaching losses were on an annual timestep, correlations were based on only 28 data points. Therefore, longer lags identified with cross-correlation will rely on fewer data points and require caution. A finer seasonal scale 24 may have shown better results. Additionally, if gross nitrate-N leaching loss is relatively linear through time, then identifying the correlation between changes in the two series will be difficult for the catchment. Thirdly, this approach overlooks the cumulative nitrate-N leaching loss; riverine nitrate-N load is likely influenced by all the losses up to the characteristic lag time, not just nitrate-N losses at that lag.
While the cumulative GAM approach connects nitrate-N leaching losses to riverine loads in an adaptive statistical model, it too has limitations in estimating lag times. First, like the cross-correlation approach there are limited data, i.e. up to 28 annual values for each catchment. Second, our models only approximate the relationship between nitrate-N leaching losses and river loads. Our model 4 attempts to increase realism by including lagged surplus precipitation as a predictor. Ideally, leaching losses would take the form of an interaction term between subsurface hydrologic flow and nitrate-N leached from the root zone and soil profile, thereby accounting for biogeochemical processes regulating nitrate-N availability but also the timing of subsurface flows. Rather, our data only permitted us to include the separate effects of cumulative excess precipitation in 4 . Third, in some cases, it was not possible to identify the smooth term f 4 in 3,4 due to limited observations or circumstances leading to poor identifiability (e.g., weak but consistent linear trends in leaching losses). Finally, our nitrate-N leaching losses consider denitrification in the topsoil but ignore the likely removal of nitrate-N in other pathways before leaving the catchment 34 . However, since we use cumulative catchment nitrate-N leaching losses to model the variability in riverine loads, we do not require a nitrate-N mass balance. Hence, we assume that variation in nitrate-N leaching losses is predictive of changes in river loads and are most predictive at the characteristic lag time for the catchment.
Impact on local policy development. Policy to protect water quality requires landowners and managers to put in place actions to mitigate leaching losses of nitrate-N to decrease the load of nitrate-N in receiving waterbodies. In New Zealand, a recent analysis found that 43% of agricultural land was in catchments where the current load exceeded the maximum allowed 35 . Other work 36 estimated that, had farming mitigation practices over 1995-2015 not been adopted, 45% more N (largely as nitrate-N) would have been lost. Despite these efforts, the expansion of intensive land uses has increased N loads by 25% nationally 36 . Where it was assumed all actions to mitigate N losses were adopted, additional modelling showed that future N loads could decrease by about a third 37 , reducing the area still exceeding the maximum allowed to about 5% 38 . However, this modelling assumes that actions were adopted over a period of 20 years, commensurate with the mean rate of adopting agricultural practice in Australasia of about 17 years 39 . Government policy aims to show improvements in water quality metrics within five years and to make waterways healthy within a generation 4 . Although policy can enforce action to occur quickly, it is still likely that their implementation and effectiveness will take time to reach their full potential. Our work would suggest that, with a median lag time in nitrate loads of 4.5 years, targeted improvement would not be possible in nearly half of our catchments, which are representative of agricultural land use in New Zealand. It is likely that small catchments or sub-catchments of those with longer lag times would respond quicker. Improving the estimation of such changes could be aided with adjustments of the current monitoring network 40,41 .

Materials and methods
Nitrate-N leaching estimates. Annual estimates of nitrate-N leaching loss were generated nationally using the method of Dymond et al. 42 . Briefly, this method combines farm level data for livestock class (beef cattle, dairy cattle, sheep and deer) and livestock numbers, collected approximately during agricultural censuses (1994, 2002, 2007, 2012, and 2017) 43 , with modelled estimates of nitrate-N loss from the root zone and soil profile (via the model OVERSEER 44 ) for those livestock numbers and classes across 100 unique soil by climate combinations identified at level II of the Land Environments of New Zealand spatial database 45 . Annual counts of livestock types in-between census years were taken from district level data (n = 53). These data were allocated to pastoral land uses according to the Land Cover Database (1997/98, 2001/02, 2008/09, 2012/13, 2018/19) 46 and then to properties using AgriBase 47 . It was assumed that changes between census years were proportional to livestock classes and numbers as indicated across farm types within a district. For example, if dairy cattle numbers increased by 50% between two censuses, the increase was apportioned equally to all dairy farms identified by the Land Cover Database and AgriBase within that district.
Arable and horticultural land uses occupy small areas of New Zealand (c. 0.5Mha cf 11.4 M ha for cattle, sheep and deer 48 ). Few data exist or can be modelled successfully for the wide range of arable and horticultural crops and crop-rotations used year-to-year. Hence, we used a constant nitrate-N loss estimate of 30 kg N ha −1 yr −1 for arable and horticultural land based on the median of field studies 49 . Inputs from native forestry were set at 1.5 kg N ha −1 yr −122 .
Catchment concentration and discharge data. We calculated annual nitrate-N loads for 77 sites from 1989-2020 from the National River Water Quality Network run by the National Institute of Water and Atmos- www.nature.com/scientificreports/ pheric Research (NIWA) and Regional Councils in New Zealand. The Network is located on 48 of New Zealand's rivers covering a range of flow regimes, catchment characteristics, and land use. Some rivers contain more than one site. We refer to a site as a river with its own catchment. No significant point sources are included on the network. The approach used a GAM to predict daily loads from monthly nitrate-N concentration measurements and daily mean discharge, accounting for the time of year and flow regime. Daily loads were summed to annual loads. All monthly nitrate-N concentration data were sourced from NIWA. Data were also secured from NIWA for daily mean discharge (calculated from 15 min observations) at each of these rivers from 1989-2010 and for 33 rivers from 2010-2020. Discharge for the remaining 44 rivers from 2010-2020 were sourced from a combination of Regional Councils and hydroelectric power producers. A description of the rivers, methods used, and quality of the data are available elsewhere 6,50,51 .
Gaps in the stream discharge records were < 1% of all data with maximum lengths of consecutive missing data < 40 days for 63 sites but 51 to 432 days for 9 rivers. To infill gaps we imputed values using the 'GR4J' hydrological model 52 via the 'airGR' package in R 53 . This rainfall-runoff model used daily gridded precipitation and potential evapotranspiration (Penman method), sourced from the National Virtual Climate Station network 54 that were subsequently averaged across catchments with inverse distance weighted interpolation.
Discharge flow components ('baseflow' and 'quick flow') were identified at each river with the hydrograph separation technique detailed by Pelletier and Andréassian 55 . We note that some catchments are subject to dam regulation and/or glacial melt (e.g., in the Otago region), meaning that 'baseflow' for these catchments will be an arbitrary, slow component of mean daily discharge (Q); consequently we de-emphasize baseflow effects for such catchments. For the entire record, baseflow fraction was calculated as daily baseflow divided by daily Q.
Estimating catchment nitrate loads. We calculated nitrate-N loads (kg N d −1 ) for dates with grab samples as nitrate-N concentration times daily Q. Generally, concentrations were above detection limits (1 µg N L −1 ); however, three rivers had 20-46% of observations at or below detection limits (up to ~ 7% of observations for 31 other rivers). While censored values can bias statistical models at low ranges of concentration 56 , we consider this a minor problem for our objective of determining annual nitrate-N loads and long-term trends in stream nitrate-N load regimes, where the majority of load is typically delivered during storm events with concentrations orders of magnitude greater than detection limits (see Supplementary Fig. S3). Hence, we used half the detection limit for censored observations and accept the minor amount of bias this has on our load models.
To model stream nitrate-N loads, we used concentration data to fit GAMs 57 for each river. Our approach is closely related to the framework of Hirsch et al. 58 but allows more flexible model building and testing. We explored several variants of the model, but found the following to be most generally applicable: where µ is the conditional expected nitrate-N load, the link function used (g) is the log-link, α is an intercept term, f 1 -f 3 are smooth functions (see below) of the predictors, and the observed data (y) are modelled as gamma-distributed with mean µ and scale λ. Notably, the gamma distribution with the log-link: (1) accounts for the heteroskedasticity common in many water quality data, since λ can vary linearly with the magnitude of the load (µ), and (2) easily allows for predictions of nitrate-N load on the original response scale without the need for bias-correction when back-transforming predictions on the log scale 59 . Alternative distributions (e.g., log-normal) yielded poorer fits.
We chose f 3 to be a thin-plate regression spline function of time (t) to account for non-linear and variable trends in nitrate-N load over the ~ 30-year period. We modelled seasonality with a cyclic-cubic spline (f 2 ) of dayof-year (DOY; 1 to 365/366). Finally, since nitrate-N concentrations typically vary with flow (Q) but also depend on whether the flow is derived from a recent storm event (low baseflow fraction) or a longer-term recession (higher baseflow fraction), we modelled this interaction between baseflow fraction and log-Q with a bivariate tensor product smooth (f 1 ). This term comprises a thin-plate regression spline for log-Q and a cubic-regression spline for baseflow fraction plus the interaction between the two smooth functions. For choice of smoother used, we generally opted for the default thin-plate regression spline or, if the predictor was cyclical, the cyclic-cubic spline but opted for a cubic-regression spline for some terms (e.g., for baseflow fraction in the tensor product smooth) to reduce computational cost; while thin-plate regression splines are more robust 57 , the choice of a cubicregression spline for well-constrained variables such as baseflow fraction did not significantly influence the fit.
Through prior experience, literature review 56 , and exploratory analyses, we deemed this model structure sufficient and robust in capturing the dominant features of nitrate-N load in these streams. More complicated models (e.g., where the concentration-discharge relationship captured in f 1 could also itself vary with time 58 ) are possible, data permitting, but, in general, few cases warranted the more complex structure, delivering only marginal improvements in predictive performance (below) over the more parsimonious model (this is explored and output in the Figshare repository: see Data availability section). This simpler GAM fit had significant effects (using the approximate p-value for GAM smooth components 60 ) for all terms for 70 of the 77 rivers; the remaining rivers had marginally low p-values (up to p = 0.15) for either the time trend or seasonal component but, for simplicity, these minimally influential terms were kept in the model for these rivers.
The GAMs were fitted via restricted maximum likelihood (REML) with the 'mgcv' package within R 61 . We assessed each fit with various residual diagnostic plots. When trialling alternative model fits, we compared models via the Akaike information criterion (AIC), residual diagnostics, and posterior simulation behaviour (e.g., to check for over-fitting). www.nature.com/scientificreports/ Annual nitrate-N loads were calculated by predicting daily nitrate-N load across the record and summing for each year. The uncertainty in these annual loads was estimated by drawing 1000 random sets of parameters from the posterior distribution of the parameters in the original GAM, re-estimating the annual load with these simulated parameters, and summarizing these simulations with a 95% credible interval 57,61 . Calculating lag-times between nitrate leaching and catchment loads. To estimate the potential lag time between nitrate-N leaching losses and riverine nitrate-N loads, we conducted three different analyses: (1) cross-correlation, (2) cumulative GAMs, and (3) assessed mean hydrologic transit times from modelling and isotope data. Each has its strengths and weaknesses (see "Discussion"), but we synthesize all three approaches when discussing lags in nitrate-N delivery in our catchments.
Cross correlation analysis. Recent work 12,62 assessed time lags in catchment N loads by pairing time series of catchment N inputs and N outputs (usually on an annual time step) and calculating the cross-correlation function for each lag time of interest (in years). This cross-correlation is the Pearson correlation coefficient between annual catchment nitrate-N load, y t , and the preceding (lagged) annual catchment nitrate-N leaching losses, x t-k , where k is lag in years from 1 to some maximum, k max . Hence, if changes in the riverine nitrate-N load regime lags changes in catchment nitrate-N leaching losses by k = 5 years, we expect the cross-correlation to reach a maximum near k = 5. Here, our x t series is the annual nitrate-N leaching losses and our y t is the annual nitrate-N load normalized for effects of all variables except time, i.e. analogous to the 'flow-normalized' loads in Hirsch et al. 58 , which corrects the nitrate-N loads for the variability in those exogenous predictor variables.
It is necessary to 'pre-whiten' both time-series before calculating the cross-correlation 63 (though some authors argue against this 12 ). Pre-whitening ensures that potential autocorrelation and non-stationarity within the two series do not falsely produce significant cross-correlations. To pre-whiten, we first-order differenced both x t and y t and then fitted an auto-regressive integrated moving average (ARIMA) model to the differenced x t . We used the automatic procedure for fitting a suitable ARIMA developed by Hyndman and Khandakar in the 'forecast' package 64 . The resultant ARIMA model was then used to filter both differenced x t and y t prior to calculating the cross-correlation 63 . We calculated cross-correlations up to a k max of 20 years.
Cumulative GAMs. Recognizing the limitations of the cross-correlation approach, we also modelled riverine nitrate-N loads using lagged, cumulative nitrate-N leaching losses as predictors by modifying the GAM in Eq.
(2) to: We substitute out the time trend component (f 3 ) for a smooth function of cumulative nitrate-N leaching losses for some lag of k years (L k ). Our hypothesis with this change is that, with other dominant features of nitrate-N loads controlled for via Q, baseflow fraction, and DOY, the main driver behind changes in the nitrate-N loads is the change in cumulative nitrate-N leaching losses. In essence, the best lagged cumulative nitrate-N leaching losses predictor, L k , will best account for the smooth time trend in the original model fit.
We also modified this GAM to account for climatic variability where, even if there is considerable nitrate-N leaching lost from the root zone and soil profile but relatively little surplus rainfall, the potential for leached nitrate-N reaching the stream is low: where S k is the cumulative sum of daily precipitation minus potential evapotranspiration (mm) for the preceding period of k years. Hence S k provides a reasonable proxy for hydrologic conditions in the catchment for the lag k considered.
Both f 4 and f 5 were incorporated as smooth functions in the GAMs. While we used a thin-plate regression spline for f 5 , we used a shape-constrained spline for f 4 (see Pya and Wood 65 ) to avoid the possibility of negative loads. We found that unconstrained fits for the L k smooth could sometimes be unrealistic (e.g., decreasing riverine nitrate-N load despite greater L k )-this was usually due to lack of observations to inform the fit. Based off previous studies of catchment nitrate-N dynamics 66-68 , we would strongly expect riverine nitrate-N loads to only increase or level off with increasing L k . We therefore enforced this constraint by modelling f 4 as a monotonicincreasing function with the 'scam' package in R 65 .
Using the cumulative GAMs Eqs. ( 3) and (4), we test k lags for the L k predictor of 1 to 10 years. At a given lag, we fit Eqs. (3) and (4) where L k and S k are their cumulative sums for that lag. As a comparison, we also re-fit our original model Eq. (2) since the underlying data available becomes shorter as the considered lag becomes longer. We extract the effective degrees of freedom and approximate p-value for the smooth term for L k (f 4 ) in Eqs. (3) and (4) as well as all models' AIC. We then make relative comparisons between these models at each lag by subtracting the AIC of Eqs. (3) and (4) from the AIC of Eq. (2): a 'rule of thumb' here is that a model outcompetes another candidate model when its AIC is more than 2 units lower 69 . We note that lags with wellperforming models of either Eqs. (3) or (4) tended to have effective degrees of freedom of f 4 closer to a value of 1 (nearly linear): these model fits favoured simple functions of L k for predicting riverine nitrate-N loads.
For each river, we present the mean of the two lags with the lowest AIC (> 2 units from the fit of Eq.( 2)) and P value < 0.1 as the output of the cumulative GAM.
Mean transit time from isotope data. Estimates of hydrologic mean transit time for surface waters were sourced from the literature. These estimates place a minimum bound for lag time for water and hence nitrate-N for nine www.nature.com/scientificreports/ of the catchments. They were derived by a combination of different model types (e.g., binary mixing, dispersion, and exponential-piston flow) and isotope (e.g., 3 H, 18 O) data from water samples taken at baseflow for eight catchments and a combination of baseflow and stormflow for one catchment (CH4).

Presentation of lags.
We present mean lag times calculated from significant fits (P < 0.1) of the crosscorrelation and cumulative GAM functions. We did not calculate lag times for those catchments where intensive agriculture was < 10% of the catchment's area. These catchments are under MDC (viz. reference conditions) and are unlikely to exhibit significant nitrate-N leaching 70 . Land use in these catchments has not changed > 0.25% over the period of monitoring (Table 2 and Supplementary Fig. S2). We also filtered out rivers that were either impacted by hydroelectric schemes or had significant land works or forestry harvesting that caused uncharacteristically high sediment load 50 (Table 2).
To model mean lag times as a function of catchment attributes, we obtained data for 13 variables that described climate, hydrology, and land use at a catchment level (Table 1). These data were used in a best subsets regression to output an 'optimal' candidate based on maximising the adjusted R 2 , minimising the AIC, and exhibited a Mallows Cp value that closely matched the number of predictor variables to avoid overfitting.
In addition to calculating and predicting lag times we calculated the cumulative nitrate-N retained in a catchment over the period of record as the difference between annual nitrate-N leached relative to annual nitrate-N load in the river as: Nitrate-N retained is inclusive of loss processes such as denitrification in aquifers 71 or uptake by stream vegetation and benthic sediments 72 but excludes input processes such as atmospheric deposition. A similar approach has been used at the catchment level in other regions and globally 17,73 . We used the same input variables and analysis used to predict mean lag times to predict nitrate-N retained.

Data availability
R code for nitrate-N loads and lag-times and the datasets generated during the current study are available in the