Time lags in watershed-scale nutrient transport: an exploration of dominant controls

Daily discharge data was obtained from the Water Survey of Canada.[42] N concentration data was obtained from the Ontario Provincial Water Quality Monitoring Network (PWQMN; Ontario Ministry of Environment). PWQMN monitoring stations were chosen for the current analysis based on the following criteria: (a) location with the GRW; (b) proximity to a Ministry of Environment (Canada) flow-monitoring station with temporally corresponding discharge data; and (c) available sampling data over a period of at least 25 years. The longest record length for any of the stations was 51 years (1964–2015), and 87% of stations had data available over a period of at least 35 years.

We used the weighted regression on time, discharge, and season (WRTDS) modeling approach [43] to estimate flow-weighted concentrations (FWCs) from intermittently measured concentration and continuously measured flow data. Daily concentration values were calculated in WRTDS via the following equation:

$$\begin{equation} \begin{array}{*{20}l} {\ln (c)} \hfill & = \hfill & {\beta _0 + \beta _1 t + \beta _2 \ln (Q) + \beta _3 \sin (2\pi t)} \hfill \\ {} \hfill & + \hfill & {\beta _4 \cos (2\pi t) + \varepsilon } \hfill \\ \end{array} \end{equation} \tag{ 1 }$$

where, c is concentration [ML⁻³], β₀,β₁,β₂,β₃,β₄, are fitted regression coefficients, Q [L³T⁻¹] is daily streamflow, t [T] is time, and ε is an error term. A close 1:1 relationship was obtained between observed and predicted concentration values (mean error = −0.01 mg mg NO₃-N/L), with a mean percent bias for all stations of −0.01 (SI, part I, tables S1.1 and S1.2).

Seasonal and annual FWCs were calculated from measured daily discharge and estimated daily concentrations (equation 5.4) using the following equation:

$$\begin{equation} C_f = \frac{{\sum\nolimits_{i = 1}^n {Q_i C_i } }}{{\sum\nolimits_{i = 1}^n {Q_i } }} \end{equation} \tag{ 2 }$$

where the numerator is the annual or seasonal N load and the denominator is the total annual or seasonal discharge. Seasonal FWC values were calculated for the winter (December, January, February), spring (March, April, May), summer (June, July, August) and fall (September, October, November) seasons.

Trajectories of net anthropogenic N inputs (NANI) were developed for the period 1901–2011 for the GRW subwatersheds, as described in section 2.2. Across the GRW, annual NANI values nearly doubled, from 46 kg N ha⁻¹ yr⁻¹ to 87 kg N ha⁻¹ yr⁻¹, between 1901 and 1976, and then decreased to 72 kg N ha⁻¹ yr⁻¹ by 2011 (figure 2). NANI trajectories for the GRW subwatersheds show similar patterns, peaking between 1976 and 1980 and then subsequently decreasing (figure 2). The extent of that decrease has ranged from 10.5%–39.8%, with a median percent decrease of 16.5%. The decrease in N inputs can primarily be attributed to an approximately 10% decreases in atmospheric N deposition as well as increases in the N use efficiency associated with agricultural production [48].

Zoom In Zoom Out Reset image size

Flow-weighted nitrate concentrations show an increasing temporal trend in the 1964–1992 timeframe across all the 16 subwatersheds analyzed in this paper (SI, part 4, table S4). These increases were found to be significant at a 99% confidence level (p < 0.01) for all but three of the subwatersheds. For the period 1993–2011, however, concentration trends have varied across the GRW. Only 3 of the 16 subwatersheds show a decreasing trend for nitrate FWCs over this period, while 13 continue to show an increasing trend. As a whole, the data show that nitrate FWCs increasing before 1993 at a mean rate of mg l⁻¹ y⁻¹; after 1993, concentrations continue to increase in the mean, but at reduced rates (mg l⁻¹ y⁻¹) and with a greater distribution of behavior across the different subwatersheds.

Spatial analysis of the GRW subwatersheds for the period 2000–2010 shows a strong linear relationship between NANI (inputs) and mean flow-weighted nitrate concentrations (outputs) (R² = 0.64, p < 0.001). However, a look at relationships between inputs (NANI) and outputs (nitrate FWC values) over time for individual subwatersheds indicates a disruption in the linearity of the input/output relationship.

As seen in figure 3(a), annual NANI values for the Grand River at Brantford began to decrease in the late 1970s, but there has been no decrease in annual nitrate FWCs since that time, although the slope value for concentration vs time for the post-1993 period (0.022 ± 0.007 mg l⁻¹ y⁻¹, R² = 0.35) has decreased somewhat from that before 1993 (0.078 ± 0.002 mg l⁻¹ y⁻¹, R² = 0.98). For the Lower Canagagigue, however (figure 3(b)), the response to changes in NANI have been faster, with a peak for flow-weighted nitrate concentrations being seen in 1997 at 8.3 mg l⁻¹, and then with significant decreases in FWC values for the post-1993 period of 0.120 ± 0.017 mg l⁻¹ y⁻¹.

Zoom In Zoom Out Reset image size

The results described above are suggestive of a mismatch between annual NANI values and annual nitrate FWCs. For example, in figure 3(c) we see a positive linear relationship between NANI values and N outputs (annual FWC) for the slower-responding main stem of the Grand up until approximately 1978 (R² = 0.93). From 1979–1986, however, there is a period of nonlinearity during which NANI values remain relatively constant, while FWC values continue to increase. Beginning in 1987, a negative relationship develops between current-year NANI values and flow-weighted nitrate concentrations. This negative relationship between current inputs and outputs indicates a potentially lagged relationship and thus creates a visible hysteresis effect in the response curve of NANI versus nitrate concentration values. For the Lower Canagagigue (figure 3(d), the hysteresis loop is tighter than that for the Grand River at Brantford, and by 2010 shows a return to late 1970s-level nitrate FWCs. The tighter hysteresis loop for the Lower Canagagigue suggests that this smaller tributary responds more quickly to changes in N inputs than the main stem of the Grand, meaning that current nutrient dynamics in the watershed are less driven by N legacies and thus that time lags will be shorter.

The cross-correlation results provide us with estimates of time lags between changes in annual NANI values and subsequent changes in flow-weighted nitrate concentrations. Annual lag times for the 16 study watersheds were found to range from 12–34 years, with a median value of 25.5 years (SI, part 5, table S5).

Lag times were generally found to be the longest in the northern subwatersheds (30–33 years), decreasing in the middle GRW, and then increasing again in the lower GRW. The longer lag times in the northern watersheds correspond to areas of low permeability and high surface runoff, but also to areas with lower NANI values (figure 4), less cropland, higher fractions of forested area, and little tile drainage (figure 1). In contrast, the shorter lags in the middle GRW correspond to areas with more tile drainage as well as larger urban areas and their associated WWTPs. In the lower GRW, which is more dominated by groundwater flows, lag times again increase, with slow subsurface transport of nitrate via groundwater pathways likely contributing to the longer lags.

Zoom In Zoom Out Reset image size

Time lags in the GRW also show seasonal variation. As an example, we show in figure 5 the monthly nitrate FWC trajectories superimposed against NANI values for the Grand River at Glen Morris Bridge. Our cross-correlation analysis of annual NANI and FWC values indicates a lag time of 20 years for this site. The monthly trajectories shown in the figure, however, show a range of behaviors across the year. In the summer, concentration values are increasing consistently over the study period, with little or no response to the decreasing post-1976 net N inputs. In the fall and winter, however, the watershed appears much more responsive to decreases in N inputs, with distinct peaks and subsequent decreases in the N concentration trajectories. In fall and winter, for example, the slope values for the pre-1993 period (0.077 ± 0.002 mg l⁻¹ y⁻¹, R² = 0.98, fall; 0.105 ± 0.005 mg l⁻¹ y⁻¹, R² = 0.94, winter) are significantly different from those for the post-1993 period (0.002 ± 0.006 mg l⁻¹ y⁻¹, R² = 0.01, fall; 0.017 ± 0.014 mg l⁻¹ y⁻¹, R² = 0.08, winter). In summer, however, the slope values for the pre- and post-1993 periods are statistically indistinguishable (0.060 ± 0.002 mg l⁻¹ y⁻¹, R² = 0.98, pre-1993; 0.041 ± 0.002 mg l⁻¹ y⁻¹, R² = 0.95, post-1993). Across all the subwatersheds, seasonal cross-correlation analysis indicates that median N-related time lags are the longest in summer (median 26.0; iqr = 19.8–28.5), although variability among stations is significant. In the following section, we explore the dominant controls on annual as well as seasonal time lags to further elucidate the observed spatial and temporal variability in the identified lag times.

Zoom In Zoom Out Reset image size

The GRW is home to seven major reservoirs, suggesting that reservoir dynamics, which increase both water and nutrient residence times, may exert a major control on lag times. Our preliminary analysis suggests that the presence of a major reservoir above a sampling station can significantly increase seasonal variability in lag times, possibly due to reservoir operations, which store water during high flow times and release it during periods of low flow. Such an impact can be seen when comparing the Marsville and Shand Dam subwatersheds, both located on the main stem of the GRW (figure 4). Although the annual time lags are similar for these two stations (Marsville, 33 years; Shand Dam, 30 years), there is much more seasonal variability at Shand Dam (figures 6(a) and (b)). By most measures, these nested subwatersheds are quite similar, e.g. cropped area (33%), NANI (48 kg ha⁻¹), percent tile drainage (∼1%), and watershed slope (∼5°). Notably, however, the Shand Dam monitoring site is directly below Belwood Lake, a 12 km long reservoir with a storage capacity of 63.9 million m³ [52] while the Marsville station is unaffected by reservoir dynamics. Accordingly, as represented by the orange dashed lines in 6(a) and 6(b), the base flow index (BFI) shows little variation across the four seasons at Marsville, while BFI values at Shand Dam site show a strong peak in summer, and values for summer and fall are substantially higher than those for Marsville (>0.7, Shand Dam; ∼0.5, Marsville). During this period, river flows at the Shand Dam site are dominated by controlled releases from the upstream reservoir, thus artificially elevating summer BFI values. The influence of the reservoir on summer flows results in greater variation in the estimated time lags throughout the year (CV = 0.02, Marsville; CV = 0.14, Shand Dam), with the peak at Shand Dam occurring in summer.

Zoom In Zoom Out Reset image size

WWTPs, especially in areas with growing urban populations, can also have an important impact on summer concentration trajectories and on the calculated lag times. In the summer, wastewater N constitutes a larger portion of stream N than at other times of the year due to lower levels of discharge and higher uptake of non-point source nutrients via crop growth. The Speed River monitoring site, for example, is heavily impacted by a large upstream WWTP, which annually releases 2376 tons of nitrate into the river; summer flow-weighted nitrate concentrations at the Speed River station are 3.5 mg l⁻¹, with approximately 75% of that summer nitrate mass (2.7 mg l⁻¹) being associated with human waste (National Pollutant Release Inventory, Canada). The summer dominance of wastewater N inputs can directly impact lag times, as WWTP discharge into a stream is not impacted by flow-path variations and long-term retention with the landscape; in other words, wastewater N has zero lag. The relatively shorter summer time lags between decreases in net N inputs and changes in water quality for the Speed River (figure 6(c)), as well as those for other heavily WWTP-impacted watersheds such as the Nith (figure 6(d)), may reflect the more direct, and essentially lag-less, relationship between wastewater N and stream nitrate concentrations.