Concentration-discharge relationships derived from a larger regional dataset as a tool for watershed management.

Concentration-discharge (C-Q) relationships have been widely used to assess the hydrochemical processes that control solute fluxes from streams. Here, using a large regional dataset we assessed long-term C-Q relationships for total phosphorus (TP), soluble reactive phosphorus (SRP), total Kjeldahl nitrogen (TKN), and nitrate (NO3 ) for 63 streams in Ontario, Canada, to better understand seasonal regional behavior of nutrients. We used C-Q plots, Kruskal-Wallis tests, and breakpoint analysis to characterize overall regional nutrient C-Q relationships and assess seasonal effects, anthropogenic impacts, and differences between "rising" and "falling" hydrograph limbs to gain an understanding of the dominant processes controlling overall C-Q relationships. We found that all nutrient concentrations were higher on average in catchments with greater levels of anthropogenic disturbance (agricultural and urban land use). TP, SRP, and TKN showed similar C-Q dynamics, with nearly flat or gently sloping C-Q relationships up to a discharge threshold after which C-Q slopes substantially increased during the rising limb. These thresholds were seasonally variable, with summer and winter thresholds occurring at lower flows compared with autumn and greater variability during snowmelt. These patterns suggest that seasonal strategies to reduce high flows, such as creating riparian wetlands or reservoirs, in conjunction with reducing related nutrient transport during high flows would be the most effective way to mitigate elevated in-stream concentrations and event export. Elevated rising limb concentrations suggest that nutrients accumulate in upland parts of the catchment during drier periods and that these are released during rain events. NO3 C-Q patterns tended to be different from the other nutrients and were further complicated by anthropogenic land use, with greater reductions on the falling limb in more disturbed catchments during certain seasons. There were few significant NO3 hydrograph limb differences, indicating that there was likely to be no dominant hysteretic pattern across our study region due to variability in hysteresis from catchment to catchment. This suggests that this nutrient may be difficult to successfully manage at the regional scale.

Introduction

In riverine ecosystems, nutrient exports coupled with the catchment’s land use and land cover are used to assess the potential for eutrophication in downstream ecosystems (Bridgeman et al. 2012, Schindler et al. 2012). In order to implement mitigation techniques to prevent environmental degradation, watershed managers must take into account not just the source of nutrients but also the fact that the allochthonous contribution of these nutrients is inextricably tied to hydrology (e.g., Banner et al. 2009, Lapworth et al. 2013, Moatar et al. 2017). Concentration‐discharge (C‐Q) relationships therefore provide valuable information on processes controlling nutrient availability and connectivity in catchments (Creed et al. 2015, Ali et al. 2017, Moatar et al. 2017). Understanding these processes is especially important in light of the increasing frequency of extreme climatic events (Burn and Whitfield 2016) which will increase hydrological nutrient export (Bosch et al. 2014).

Solutes are typically delivered to stream systems through event‐driven (storm and snowmelt) transport from the terrestrial environment through subsurface and overland flows. This “new” water from storm and snowmelt events mixes with or displaces older soil water in which nutrients have accumulated in both dissolved and particulate form, delivering them to the stream through both types of flows (e.g., Buttle 1994). Overland flows also erode the soil surface, transporting sediments and associated nutrients. The resulting stream solute levels will either be chemostatic (concentrations are discharge‐invariant) or chemodynamic (concentrations increase or decrease in synchrony with discharge). Chemostatic behavior has been observed for nitrate (NO3) and phosphorus (P) in long‐term intensive agricultural catchments (Basu et al. 2011). Negative chemodynamic behavior is due to some degree of dilution from incoming flows (Ali et al. 2017, Moatar et al. 2017), while positive chemodynamic behavior can be the result of in‐stream removal and retention processes during low flows and terrestrial solute export during high flows from catchment connectivity (Ali et al. 2017, Moatar et al. 2017). Connectivity is an important concept in flooding; an increase in storm or melt water causes an increase in discharge through water and associated nutrients transported from source pools within the catchment, source pools that are effectively isolated from the stream during drier periods. Chemostasis or slight negative C‐Q relationships can also be attributed to weathering (e.g., Godsey et al. 2009). For all solutes the processes influencing C‐Q relationships may shift throughout the discharge range. Moatar et al. (2017), who conducted a regional study in France, reported that C‐Q slopes frequently differ on either side of a discharge threshold, while a prairie study reported models that incorporated a discharge threshold that tended to better describe C‐Q patterns (Ali et al. 2017). Several other studies have also reported discharge thresholds (e.g., Banner et al. 2009, Murphy et al. 2014). These thresholds indicate a change in the dominant process influencing nutrient concentrations and may be caused by a change from baseflow to stormflow.

Another important part of chemodynamic relationships is the differences that occur because of increasing (rising hydrograph limb) or decreasing (falling hydrograph limb) discharge. C‐Q relationship differences between hydrograph limbs are common and have been extensively reported in the literature (e.g., Cooke and Cooper 1988, Raymond and Saiers 2010, Siwek et al. 2013, Lloyd et al. 2016, Winnick et al. 2017). These studies and others have found a wide range of hysteretic patterns that can be broadly categorized as clockwise (rising limb concentrations are greater than those on the falling limb) and counterclockwise (greater falling limb concentrations). The former can be an indication of nutrient sources that are hydrologically connected and relatively close to the stream, while the opposite is suggested of the latter (Creed et al. 2015).

Most studies use high‐resolution data for a single or small number of catchments and focus on event‐specific hysteretic effects (e.g., House and Warwick 1998, Evans et al. 2003, Lloyd et al. 2016), the pattern of which depends on various catchment and event‐specific factors (e.g., Siwek et al. 2013, Feinson et al. 2016). Despite the variability of hysteretic patterns, an overall pattern may emerge when looking at lower‐resolution, longer‐term datasets that combine multiple sites using a big data approach (Moatar et al. 2017). Raymond and Saiers (2010), for example, examined dissolved organic carbon (DOC) C‐Q relationships separated by hydrograph limb from over 5,000 measurements of 30 small, forested sites taken over 31 years. By classifying discharge measurements into “bins” and averaging the discharge and concentrations within each bin, they determined that DOC concentrations were higher on the rising limb compared with the falling limb. Here, we used a similar method but expanded our analysis to include several nutrients and evaluated groupings of catchments of variable sizes and differing land uses.

Our objective was to characterize the overall C‐Q relationships, including threshold behavior, of total phosphorus (TP), soluble reactive phosphorus (SRP), total Kjeldahl nitrogen (TKN), and NO3 in 63 Ontario, Canada catchments in order to make inferences on the dominant sources and processes controlling nutrient concentrations. Our long‐term, low‐resolution dataset consisted of up to 14,450 measurements taken from periods of 10 to 50 years between 1964 and 2014. This type of data, in addition to being readily available, avoids much of the site and event‐specific hysteretic heterogeneity and complexity inherent with high‐resolution datasets. Analysis of this dataset will provide a general understanding of the processes influencing hydrology and solute biogeochemistry, which will be beneficial (Kirchner 2003). In this study we described a regional picture of C‐Q patterns and made inferences on what these patterns may mean for nutrient management. We divided the data by hydrological season and standardized all concentration and discharge data then, using a bin‐averaging technique, we assessed regional C‐Q behavior on the rising and falling limbs of the hydrograph. We predicted that C‐Q relationships of TP, SRP, TKN, and NO3 would vary depending on season, hydrograph limb, and land use. We expected that C‐Q relationships would generally be positive, but that chemostatic effects might be observed in catchments that were heavily influenced by anthropogenic land uses. We also expected to find similar seasonal patterns across nutrients.

Methods

Site selection

Water quality sites were selected from Ontario’s Provincial (Stream) Water Quality Monitoring Network (PWQMN; Ministry of the Environment and Climate Change 2015) based on proximity to a stream gauging station and the length of the data record. Stream gauging stations (Environment and Climate Change Canada 2014) were not usually co‐located with water quality sites, therefore we paired our sites with gauging stations only if they were on the same main channel, not separated by a dam, and if their catchments overlapped by at least 75%. We used the 75% overlap such that water quality sites and gauging stations were located closely enough that the effect of landscape features which may alter flow or nutrient patterns, such as riparian wetlands, would be minimized. The overlap on average was 96% (± 5% standard deviation). Catchments were delineated using the Ontario Flow Assessment Tools III (OFAT III; Ontario Ministry of Natural Resources and Forestry 2015) and dam locations were determined from the Ontario Dam Inventory (Ontario Ministry of Natural Resources and Forestry 2014). Additionally, only sites with a range of at least 20 years of discharge data overlapped with at least 10 years of water quality data were included in our analysis. The record length for water quality data was chosen to ensure a large number of samples could be included for each site. We considered using concentration means for each year and/or discharge exceedance level for each site; however, due to the low resolution of our dataset this meant that different sites contributed data unevenly. We therefore found it necessary to combine all data across all years in order to have enough data at the extreme ends of the discharge range to facilitate analysis. In total, there were ˜11,100–14,450 measurements depending on the nutrient (Appendix S1: Table S1). Even though some sites were sampled more frequently and for longer periods than others, we included them in the analysis, because they were often among only a handful of sites that contributed data at very‐low‐discharge and very‐high‐discharge bins. In addition, when excluding these sites entirely we found that the overall patterns of C‐Q relationships were similar to those found when retaining all sites. We therefore retained all data.

The catchments we studied varied in topography, land use, and hydroclimate (Table 1). It is well known that these factors influence stream nutrient levels (e.g., Li et al. 2006, McCullough et al. 2012, Becker et al. 2014), therefore a site‐specific standardization was applied to make values comparable across the various types of catchments (discussed in the Hydrology and Nutrients subsections). We did this with the assumption that, although individual C‐Q patterns were likely to vary somewhat between catchments and events (Appendix S1: Fig. S1), a dominant pattern would emerge that would allow us to make inferences on the types of processes influencing solute levels across the discharge range; information that would be applicable to regional water quality management decisions.

Table 1

Summary of catchment landscape variables for the Ontario catchments.

Landscape variable	Minimum	Maximum	Mean	Median
Longitude (DD)	−89.21	−74.49	−80.24	−79.83
Latitude (DD)	42.16	49.61	44.54	44.17
Catchment area (km2)	32	8,532	501	210
Mean elevation (m asl)	5	475	260	271
Mean catchment slope (%)	0.41	7.65	3.42	2.99
Main channel length (km)	15	380	72	54
Mean channel slope (%)	0.04	1.22	0.28	0.17
Catchment area affected by dam (%)	0	100	32	10
Catchment agriculture (%)	0	96	61	70
Catchment forest (%)	0	96	22	10
Catchment urban (%)	0	79	11	3
Catchment wetland (%)	0	13	2	1
Catchment open water (%)	0	13	1	0.1
Mean annual temperature (°C)	1.28	9.72	6.53	6.80
Mean annual precipitation (mm)	688	1,095	913	896
Mean annual discharge (m3/km2)	43	13,357	2,751	1,657

Hydrology

Environment and Climate Change Canada (2014) stream gauging station daily mean discharge data were used to calculate discharge exceedances for each site’s period of record using standard flow duration curves in the Streamflow Analysis and Assessment Software (SAAS; Metcalfe and Schmidt 2016). We manually defined seasons using each site’s median annual hydrograph for all years in which hydrological data were available during our period of record. While median hydrographs cannot capture interannual variability they are useful for determining when hydrological seasons typically occur at each site. We estimated the date ranges for each season rounded to the nearest half‐month for each site and defined them as snowmelt (large discharge peak occurring in the spring), summer (low‐flow period), autumn (slightly elevated flow period), and winter (reduced or stabilized flow); see Fig. 1 for a visual explanation of seasonal definitions. This method resulted in uneven season lengths (typically summer was the longest season, and autumn the shortest); however, we felt that we obtained a good representation of hydrological seasons. We noted that while the seasonal date ranges and lengths varied from catchment to catchment, sites located within the same secondary catchment (larger river systems within Ontario’s three primary watersheds) generally had similar hydrographical seasonal definitions.

Fig. 1

Example of hydrologic seasonal boundary definitions, with (a) winter as the low‐flow period during the cold season, (b) snowmelt as the point where discharge increases rapidly over time following winter, (c) summer as the point when discharge begins to drop less dramatically following snowmelt, and (d) the period of increased flows between summer and winter. Data are from gauging station 02BA003, which is mainly forested (95%) with an area of 4,222 km2.

Although discharge records were mostly complete, missing data were estimated using interpolation (for gap lengths of up to three missing days which were present in eight sites) or using the values derived from the linear equation of a highly correlated site (r > 95%) within the same secondary watershed (for up to 14 consecutive missing days that were present in 10 sites). SAAS was also used to define rising and falling hydrograph limbs with rising limb days designated as days in which discharge was increasing, including the day on which discharge peaked, and falling limb days defined as those in which discharge was decreasing, including the day with the lowest discharge before flow began increasing again (Fig. 2). Samples taken on days when estimated baseflow was equivalent to mean discharge were considered to be non‐events and, as there were not enough of these to facilitate analysis and they could not be considered rising or falling limb, they were removed from our dataset.

Fig. 2

Conceptual diagram of differentiation between rising and falling hydrograph limbs used in this study. Discharge peaks were considered to be part of the rising limb, while discharge lows were considered part of the falling limb. Non‐events, in which mean daily discharge was equivalent to baseflow discharge (not shown) were not included in analysis.

Similar to the method used by Banner et al. (2009) we converted the mean daily discharge associated with each water quality measurement to the respective site‐specific discharge exceedance percentage in order to standardize these measurements across all sites. That is, we calculated discharge exceedances for each site individually in order to have standardized discharge measurements that could be combined across all sites. Discharge exceedances refer to the percentage of time that streamflow is greater than or equal to a given discharge. For example, a discharge value associated with a 10% exceedance (Q10) indicates that that discharge was observed to be greater or equal to that value only 10% of the time.

Nutrient analysis

We obtained TP, SRP, total Kjeldahl nitrogen (organic nitrogen and ammonium; TKN), and nitrate + nitrite (referred to here as simply NO3 because NO2 was <15% of total nitrates) data for 63 sites (Fig. 3) from the PWQMN dataset (Ministry of the Environment and Climate Change 2015). We recognized that the larger, more northerly sites had the potential to be outliers in our dataset; however, our overall results did not change when they were omitted from analysis. PWQMN water quality samples were collected by Ontario’s Conservation Authorities under a standard protocol (Ontario Ministry of the Environment 1983) and analyzed at the Ontario Ministry of the Environment and Climate Change (MOECC) Laboratory Services Branch. Data were collected approximately monthly, although more or less frequent sampling occurred in some years at various sites (Appendix S1: Table S2). Some sites were also monitored for longer periods than others, thus the number of sample measurements varied for each site. All nutrients were measured using various colorimetry methods. Earlier analytical methods (1964 to approximately 1990) were frequently not documented for the PWQMN data; however, we assumed that they were similar to those where analytical methods were provided (˜1980–2014; MOECC analytical code E3367A and equivalents for TP and TKN, and E3364A and equivalents for SRP and NO3). It should be noted that both SRP and NO3 were first filtered until the 1980s, thereafter the supernatant of a settled sample was used. Through analysis, the MOECC determined that the change in method did not affect SRP. The nitrates/nitrites methodology involves the reduction of NO3 in the water to NO2 and assumes that no sediment was present in the sample extracted for measurements (Ontario Ministry of the Environment 1983). In 2013 and 2014, another method was used for some TP and TKN samples (E3516) which may underestimate concentrations in samples with elevated sediment levels. About 3% of the data were analyzed using this method, therefore we included these in the analysis as any possible underestimations were unlikely to affect our bin‐averaged results.

Fig. 3

Ontario, Canada catchments used in regional concentration‐discharge analysis.

TP, SRP, TKN, and NO3 concentrations were z‐scored by season within each site using the equation , where C is the z‐scored concentration, C is the measured concentration, and C mean and C stdev are the mean and standard deviation for the concentrations, respectively, measured at the given site in the given season. Z‐scoring in this case converts the concentration data to values that, instead of representing absolute concentration, indicate the number of standard deviations a particular concentration is away from the mean. Thus large z‐scored values are higher concentrations for that site and season while large negative values are lower compared with seasonal mean concentration levels.

Data treatment and analysis

Similar to the method used by Raymond and Saiers (2010), all data were divided into either the rising or falling hydrograph limb and bin‐averaged, which we did according to discharge exceedances (Appendix S1: Table S3). For example, all z‐scored TP measurements (across all sites) that occurred on the rising limb when discharge was at each site’s respective Q10 discharge level (10th discharge exceedance bin) were averaged together to return a single mean TP value. Bins were chosen such that high and low discharge exceedance bins were narrow to better represent high‐flow and low‐flow effects. Most water quality samples were taken in the middle of the discharge exceedance range, with more samples on the falling limb (Appendix S1: Fig. S2).

Our initial assessment of the binned z‐scored concentration‐discharge exceedance plots (visual evaluation and model fitting) showed nonlinear tendencies in C‐Q relationships, thus breakpoint analysis was run using the “segmented” package (Muggeo 2008) in R (R Core Team 2016), which estimates breakpoints using broken‐line regression models. Because our C‐Q relationships contained a maximum of 14 data points there were not enough data to facilitate linear regression or analysis of covariance (ANCOVA) tests to compare rising and falling limb slopes with any reasonable degree of power. Instead, we ran Kruskal‐Wallis (KW) tests on the complete, unaveraged binned data (at each bin level) to assess differences between seasons, as well as differences between rising and falling limb concentrations at an alpha level of 0.05 and 0.10 corrected using the Bonferroni method to reduce family‐wise type I error. We also briefly investigated land use to assess differences in C‐Q relationships at different levels of anthropogenic disturbance. As land use maps were not available for all years in our data range we calculated the percentage of each land use type for each site as the weighted average of the available years (for 1966 we used Ducks Unlimited Canada 2009; for the years 1990, 2000, 2010 we used the reference AAFC 2015; and for the years 2011–2014 we used the reference AAFC 2016). We defined anthropogenic disturbance as the sum of agriculture and urban land use, and roughly split our sites into high disturbance (where agriculture + urban > 60%, about two‐thirds of all sites) and low disturbance (≤60%, about one‐third of all sites). The 60% threshold was used because, at this percentage, the majority of the catchment land use is anthropogenic. While we recognize that there are differences in nutrient sources between urban (usually point‐source) and agricultural land (usually diffuse), both cause excessive nutrient contributions during flooding (St‐Hilaire et al. 2016); there is additionally a great deal of variability within these land use classes (for example, imperviousness and drainage systems in urban areas, tile drainage and crop type in agricultural areas) therefore we combined them for simplicity. We assessed the effect of disturbance on standardized nutrient levels at each bin level using the KW method described above for each season. Additionally, we calculated the raw concentration means for each site and used KW tests to determine whether these differed between the high‐ and low‐disturbance sites. We did this for all data across the entire discharge range as well as at only low flows (below Q90).

Results

The use of concentration‐discharge (C‐Q) relationships derived from longer‐term but relatively low‐frequency water quality sampling allowed us to assess regional patterns between nutrients and hydrology across Ontario, Canada. We expressed C‐Q relationships as bin‐averaged, seasonally z‐scored concentrations against discharge exceedance and, using this standardized dataset, we found overall positive C‐Q relationships on both the rising and falling limb of the hydrograph. There were, however, variations depending on the hydrological season. KW tests found many significant seasonal differences (Bonferroni adjusted P‐value <0.05) in concentrations for TP, SRP, TKN, and NO3 in multiple discharge exceedance bin levels (Appendix S1: Fig. S3) therefore we used seasonally separated data in further analyses.

We found similar seasonal C‐Q patterns for TP, SRP, and TKN with a few exceptions. C‐Q slopes for these nutrients tended to be nearly chemostatic or gently positively sloping below discharge thresholds around Q60, Q10, and Q40 for summer, autumn, and winter, respectively (Figs. 4, 5, 6; see Appendix S1: Table S4 for all threshold values). In terms of discharge flow values, this can be interpreted as breakpoints occurring at lower flows during the summer and winter compared with autumn, where the breakpoint occurred at high flows. We were unable to detect thresholds for snowmelt TP and TKN on the rising limb, but we did identify an SRP threshold around Q20. Falling limb thresholds for this season were variable across nutrients. Above these breakpoints, the C‐Q slopes showed a pronounced increase, particularly in autumn, and this was common across all three of these nutrients.

Fig. 4

Regional concentration‐discharge relationships for TP by hydrograph limb during each season. Concentration is represented by site‐specific z‐scored nutrient concentration and discharge is represented by site‐specific discharge exceedances. Each point is the mean of all data across all sites within the respective discharge bin. Z‐scores indicate how far above (positive values) or below (negative values) the concentration falls from the mean. Asterisks denote significant differences at a corrected alpha level of 0.05 (black) and 0.1 (gray) for individual discharge exceedance bins (Kruskal‐Wallis tests). Error bars represent standard errors. Points lacking error bars represent single measurements rather than averages.

Fig. 5

Regional concentration‐discharge relationships for SRP by hydrograph limb during each season. See Fig. 4 caption for details on axis units. Asterisks denote significant differences at a corrected alpha level of 0.05 (black) and 0.1 (gray) for individual discharge exceedance bins (Kruskal‐Wallis tests). Error bars represent standard errors. Points lacking error bars represent single measurements rather than averages.

Fig. 6

Regional concentration‐discharge relationships for TKN by hydrograph limb during each season. See Fig. 4 caption for details on the axis units. Asterisks denote significant differences at a corrected alpha level of 0.05 (black) and 0.1 (gray) for individual discharge exceedance bins (Kruskal‐Wallis tests). Error bars represent standard errors. Points lacking error bars represent single measurements rather than averages.

As extremely low (droughts) and high (large floods) flow events do not occur every year, some of the points at the extreme ends of the discharge exceedance ranges for all C‐Q relationships were composed of a single sample measurement (identified as a point lacking error bars; Figs. 4, 5, 6, 7, 8). We included these points to give full consideration to as much of the discharge range as possible in the interest of making a regional assessment. Retaining these points did not effectively influence our results, but we were cautious in the interpretation of these points when they did not follow the C‐Q pattern. For example, the summer rising limb TP concentration value for bin 13 (Q0.1) is quite a bit higher than the other points (Fig. 4), which may be uncharacteristic (an outlier). We therefore did not discuss thresholds that were driven by points derived from a single measurement.

Fig. 7

Regional NO3 C‐Q relationships for snowmelt (top row) and summer (bottom row) by hydrograph limb in catchments with low or high anthropogenic disturbance (agricultural plus urban land use). Concentration is represented by site‐specific z‐scored nutrient concentration and discharge is represented by site‐specific discharge exceedances, with the mean if all data across all sites taken within a discharge bin. Z‐scores indicate how far above (positive values) or below (negative values) the mean the concentration measurement is. Asterisks denote significant Kruskal‐Wallis tests at a corrected alpha level of 0.05 (black) and 0.1 (gray) for individual discharge exceedance bins. Error bars represent standard errors. Points lacking error bars represent single measurements rather than averages.

Fig. 8

Regional NO3 C‐Q relationships for autumn and winter by hydrograph limb. Error bars represent standard errors. Points lacking error bars represent single measurements rather than averages. Concentration is represented by site‐specific z‐scored nutrient concentration and discharge is represented by site‐specific discharge exceedances, with the mean if all data across all sites taken within a discharge bin. Z‐scores indicate how far above (positive values) or below (negative values) the mean is the concentration measurement.

Typically, rising limb C‐Q slopes were greater than those on the falling limb for TP, SRP, and TKN (Figs. 4, 5, 6). KW tests showed several differences between rising and falling limb concentration levels, largely during the snowmelt season; however, a lack of significant differences at individual discharge exceedance levels does not necessarily preclude differences between slopes. The presence of several significant hydrograph limb differences, and the fact that nearly all TP, SRP, and TKN C‐Q relationships show a similar pattern of pronounced above‐threshold rising limb concentration increases with discharge, suggests that above‐threshold rising limb levels are indeed greater than those on the falling limb. The snowmelt season is an obvious exception to this pattern, in which differences in hydrograph limbs are less distinct or more complicated.

Average nutrient concentrations were always significantly higher in catchments with >60% anthropogenic disturbance across the entire discharge range (Table 2). We did not find that any nutrient C‐Q relationships were affected by anthropogenic disturbance except NO3. While there were seasonal differences in NO3 C‐Q relationships, there did not seem to be any clear seasonal pattern such as those we found for TP, SRP, and TKN. The anthropogenic disturbance influence on NO3 C‐Q slopes only appeared to occur during snowmelt and summer, with KW tests showing significant differences (Bonferroni adjusted P < 0.05) at individual discharge bin levels on both the rising and falling limb, although the differences were more apparent on the falling limb (Fig. 7). Differences between rising and falling limb C‐Q slopes for these seasons were less apparent compared with the other nutrients. KW tests identified few significant differences between hydrograph limbs in these seasons; in low‐disturbance catchments, statistical differences were found at bins 7, 8, and 12 (Q60, Q40, and Q01, respectively) during snowmelt, but only at bin 10 (Q10) during the summer. In high‐disturbance catchments only summer showed significant KW tests for bins 7 and 8, although this latter example’s C‐Q pattern reflected those observed for other nutrients in this season.

Table 2

Nutrient concentration means for low and high‐disturbance sites, with P‐values from Kruskal‐Wallis tests.

Flow range	Disturbance level	TP (mg/L)	SRP (mg/L)	TKN (mg/L)	NO3 (mg/L)
Complete	Low	0.05 (± 0.04)***	0.01 (± 0.02)***	0.67 (± 0.87)***	0.47 (± 0.29)**
	High	0.18 (± 0.36)	0.14 (± 0.35)	1.00 (± 0.67)	2.38 (± 1.13)
Low flow (<Q90)	Low	0.04 (± 0.05)**	0.01 (± 0.03)*	0.71 (± 1.22)*	0.37 (± 0.30)***
	High	0.14 (± 0.39)	0.07 (± 0.38)	0.90 (± 0.68)	1.86 (± 2.21)

Standard deviations are in parentheses. Significance levels are as follows (for paired low and high disturbance level values): *P < 0.01, **P < 0.001, ***P < 0.0001.

We found neither anthropogenic effects nor hydrograph limb effects for NO3 C‐Q relationships in the autumn and winter. Autumn C‐Q slopes increased across the discharge range, while winter C‐Q patterns appeared to be chemostatic (Fig. 8).

Discussion

C‐Q relationships can provide information on the processes controlling stream nutrient concentrations and using a large, publicly available regional dataset (or “big data”) the dominant processes across a region can be identified, including thresholds indicating the point in the discharge range where these processes change. Using a large dataset allowed us to capture more hydrological points at the higher and lower ends of the discharge range that, at an individual catchment level, tend to be underrepresented when sampling frequency is low. From this we identified seasonal threshold behavior in TP, SRP, TKN, and NO3 C‐Q relationships which we use to suggest management strategies that are relevant to our study region.

TP, SRP, and TKN

Processes controlling above and below‐threshold nutrient concentrations

Effective management of stream ecosystems begins with understanding the processes that control stream nutrient concentrations, which includes the analysis of C‐Q relationship thresholds and hydrographic limb effects. The dominant regional behaviors of TP, SRP, and TKN were characterized by similar nonlinear seasonal C‐Q patterns and hydrographic limb effects. Below‐threshold C‐Q slopes tended to be flat or gently positively sloping in most seasons for these nutrients (Figs. 4, 5, 6). While C‐Q thresholds for low‐resolution data are rarely reported in the literature, there have been a few studies that provide threshold‐related results that suggest that flat or decreasing below‐threshold C‐Q slopes may be relatively common for these nutrients. Moatar et al. (2017) reported that more than 90% of their French sites exhibited this behavior for P below Q50, and other studies have assumed flat below‐threshold P C‐Q slopes as well (Banner et al. 2009, Ali et al. 2017). Following the conceptual framework presented by Moatar et al. (2017), we attributed this behavior to in‐stream biogeochemical processes such as decreases in biological assimilation or the release of sediment‐bound nutrients as below‐threshold flow increases. Greater variability of nutrient concentrations at low flows was also evident in our data (Appendix S1: Fig. S1) which further suggests in‐stream and transformation processes dominate rather than transport.

As our below‐threshold nutrient dynamics are not hydrologically driven (Godsey et al. 2009, Li et al. 2017), a management focus on projects such as stream channel restoration to improve benthic substrate or increased cross‐sectional area to lengthen travel time (Wilcock et al. 2002, Thompson et al. 2018) may increase nutrient retention thereby reducing ambient (low flow) in‐stream concentrations. The implementation of these types of strategies would be beneficial only if ambient concentrations are above levels of concern, which appears to be true for our sites (Table 2). For example, at low discharge (less than Q90) as well as across all flows the TP concentrations in the more disturbed sites are above the 75 µg/L eutrophic threshold suggested by Dodds et al. (1998). Eutrophic conditions suggest that low‐flow nutrient management strategies such as those mentioned above may be of some benefit to reducing ambient concentrations. It should be noted however, that even a successful reduction of baseflow nutrient concentrations would be likely to have little overall impact on annual nutrient loads. Banner et al. (2009), for example, modeled a reduction in baseflow TP concentrations and found that this reduced annual loads by less than 5%. Furthermore, the effectiveness of stream channel restoration is not certain; although studies report nutrient reductions at the reach scale (e.g., Thompson et al. 2018), others suggest that stream channel restorations which focus on in‐stream processes to reduce nutrients may not have a significant impact on nutrient concentrations (e.g., Filoso and Palmer 2011). We did not find below‐threshold dilution effects (negative C‐Q slopes) in our region; however, in these cases management targeted at below‐threshold discharge would be even less effective. In cases where a strongly increasing below‐threshold relationship exists, these low‐flow mitigation strategies may be of some benefit.

In contrast, high‐flow (above‐threshold) TP, SRP, and TKN C‐Q relationships tended to be chemodynamic on the rising limb, increasing substantially with flow (Figs. 4, 5, 6) suggesting a shift from in‐stream biogeochemical to hydrological controls on nutrient dynamics (Moatar et al. 2017). This indicates that nutrient export from terrestrial sources becomes appreciable at higher flows and the fact that this occurs primarily on the rising limb suggests that the sources of these nutrients connect relatively easily to the stream as event water increases (Creed et al. 2015). Event water delivered to streams is mainly composed of older water that has been stored in the catchment between events (Buttle 1994); newer event water displaces or mixes with this potentially nutrient‐laden stored water, delivering it to the stream in both subsurface and overland flows. A positive C‐Q slope implies that these nutrient pools are present at higher concentration in portions of a watershed that only become hydrologically connected under high‐flow conditions. This may be the result of accumulation in certain landscape features, such as wetlands, during periods of disconnectedness or of an increasing proportion of flow originating from near‐surface rather than deeper flow paths. Greater erosion from high‐intensity storms may also be transporting P through erosion, creating these positive TP C‐Q slopes. Additionally, “legacy” nutrients can be present in soils, particularly in agricultural catchments (Basu et al. 2011) which may also wash out during events causing nutrient pulses in the stream. We found that concentration peaked with discharge, with most seasonal TP, SRP, and TKN C‐Q relationships showing greater nutrient concentrations on the rising limb of the hydrograph compared with the falling limb (Figs. 4, 5, 6). This suggests that, at the event scale, clockwise hysteresis dominates. While individual studies differ as to which hysteretic pattern is most common, both clockwise and counterclockwise hysteretic patterns have been reported for these nutrients (e.g., House and Warwick 1998, Rodríguez‐Blanco et al. 2013, Lloyd et al. 2016). As event water tends to peak first followed by ground and soil water, the receding limb is likely to be made up mostly of the subsurface fractions which may have lower nutrient concentrations.

Implications for annual nutrient loads

We know that stream nutrient concentrations have a direct impact on annual nutrient loads. It has been well documented that the bulk of annual nutrient loads occur during high‐discharge events (e.g., Royer et al. 2006, Banner et al. 2009, Moatar et al. 2017) and the steeper the C‐Q slope, the larger is this annual proportion (Banner et al. 2009). We have shown evidence that a more pronounced increase in nutrient concentrations with discharge appears to occur after a threshold (Figs. 4, 5, 6), which could compound nutrient load estimates made assuming a more consistent linear increase. Therefore we are confident that our study catchments would generally benefit more from above‐threshold nutrient management that focuses on the reduction of nutrient inputs through hydrological pathways than the reduction of ambient nutrient concentrations. These might include measures that improve catchment hydrological connectivity, such as lowering and revegetating the floodplain (e.g., Baptist et al. 2004, Singh et al. 2018) and the reduction of nutrient loading through best management practices to diminish the size of upland source pools that could inhibit nutrient build‐up between events, reducing the likelihood of large nutrient pulses. This would not necessarily decrease the annual load but may spread that load more evenly across smaller events throughout the year instead of primarily during extreme events. We know that large nutrient loads, particularly P, can cause downstream algal blooms (e.g., McCullough et al. 2012), therefore these detrimental effects may be prevented by decreasing individual event loads. Other methods aimed at reducing peak stream flow through increased retainment and slowed flow velocities are likely to also be beneficial as these could give nutrients the opportunity to settle out and/or be processed before they reach the stream. Strategies might include creating or restoring riparian vegetation and wetlands or installing other types of reservoirs. Two‐stage agricultural drainage ditches have also been shown to reduce nutrient export (Christopher et al. 2017).

Seasonality of thresholds

Given the importance of high‐flow management to our catchments, it is important to consider the seasonality of high‐flow events; that is, at what level does the high‐flow threshold occur in each season? We found that the definition of high flow, in terms of discharge thresholds, did indeed vary seasonally. Banner et al. (2009) found seasonal median TP thresholds between Q32 and Q67 with similar patterns of pronounced positive C‐Q slopes above‐threshold, while Moatar et al. (2017) set Q50 as a discharge breakpoint and reported this pattern for 29% of their sites. In our study autumn breakpoints occurred at a greater discharge exceedance than the other seasons (Figs. 4, 5, 6) which indicates that nutrient contamination is not such a large concern in this season except at extremely high flows (above Q10). Additionally, the above‐threshold C‐Q slope is much steeper than in other seasons, therefore flows above Q10 are likely to have considerably high nutrient export. While the reason for this pattern is unclear, it does highlight the need for seasonally based nutrient management such as adjusting when fertilizer is applied to fields (Gildow et al. 2016). Nutrient loads or concentrations may be similarly reduced by adjusting the timing of dam reservoir or wastewater treatment plant release in urban areas.

Seasonality may also be an important factor in determining the form of the nutrients that are contributed. During snowmelt for example, TP and TKN rising limb C‐Q relationships do not show evidence of a discharge breakpoint, while the processes affecting SRP shifted at about Q20. All three of these nutrients are associated with overland flow (e.g., Cooke and Cooper 1988, House and Warwick 1998, Siwek et al. 2013, Lloyd et al. 2016); TP and TKN are composed of both dissolved and particulate fractions, while SRP is exclusively dissolved, therefore the difference between the threshold behavior (or lack thereof) of these nutrients may be due in part to nutrient form. Indeed, a Manitoba study found dissolved P to be the dominant form of P transported during snowmelt (Wilson et al. 2019), which the study speculated might be due to the lower erosive force of snowmelt compared with rainfall events. In other seasons it may be that particulate forms of TP and TKN originate in‐stream or near stream from eroded sources, while SRP is more influenced by whether overland flows are great enough to connect watershed sources with the stream. Regardless of the transport mechanism, these results suggest that while TP, SRP, and TKN may benefit from similar seasonal watershed management during most of the year, snowmelt mitigation strategies may require a nutrient‐specific approach.

NO3

A specific mitigation approach may also be required for NO3, as NO3 C‐Q patterns differ from the other nutrients, with fewer threshold effects and added complexity from anthropogenic disturbance. This difference is not surprising as, unlike the other nutrients, NO3 is largely associated with groundwater and subsurface flow pathways (e.g., Bowes et al. 2009, Lloyd et al. 2016). Autumn and winter NO3 C‐Q slopes were relatively linear, which suggests that NO3 supply is not exhausted during events. In addition, there were no statistically significant hydrograph limb differences (Fig. 7), which may be due to heterogeneity of hysteresis behavior across our catchments and events. A high‐resolution study of North Carolinian river reported that there were no dominant NO3 hysteretic patterns during most of the year (Baker and Showers 2019); autumn and winter hysteretic NO3 patterns were likely to be similarly variable in our catchments and prevented us from postulating on the dominant regional supply pathway of this nutrient. Other studies of high‐resolution data also report no clear seasonal pattern in NO3 hysteresis (Lloyd et al. 2016, Vaughan et al. 2017). Snowmelt and summer seasons, in contrast, showed a few hydrograph limb differences. These, however, tended to be found in the middle of the discharge range (Q60 to Q40) with no consistent pattern between the seasons and few similarities with the C‐Q patterns found for other nutrients. Therefore hysteresis heterogeneity is likely to extend across all seasons for NO3, although it is also possible that there are weak or absent hysteretic patterns that were not detectable from the regional, low‐frequency data we used here. If similar rising and falling limb concentrations are indeed representative of regional NO3 patterns, steps to reduce terrestrial source inputs may need to be taken, including projects such as creating riparian wetlands in areas of concern. However if there is no dominant C‐Q pattern across the province it would be difficult to form a regional management plan for this nutrient and a more catchment‐specific approach may be required.

Anthropogenic disturbance effects

Catchment‐specific levels of anthropogenic disturbance may also be important for NO3 management, however we found that the effect of urban and agricultural land, although present, was difficult to interpret. All nutrient concentrations were higher on average in more disturbed catchments (Table 2) which suggests that low‐flow management strategies would be beneficial in these areas. Of all the nutrients and seasons examined, only NO3 C‐Q relationships showed a land use effect, and only during snowmelt and summer. NO3 C‐Q slope differences between low and high anthropogenic disturbance on the rising limb were minimal (Fig. 7). When event water recedes (falling limb) however, there is a more pronounced reduction in concentration in high‐disturbance catchments as discharge decreases. This more rapid decrease in NO3 may indicate that this nutrient is primarily sourced from the terrestrial environment, with comparatively low baseflow concentrations compared with more natural catchments. This was not unexpected, as agriculture and urban land uses are both associated with elevated NO3 levels and nutrient dynamics that differ from more natural landscapes (e.g., Aguilera et al. 2012, Lapworth et al. 2013); however, we were expecting to also find some evidence of chemostatic C‐Q patterns in high‐disturbance catchments. Nutrient C‐Q relationships in intensive agricultural areas, for example, have shown chemostatic responses to flow across most of the discharge range due to “legacy” nutrients built up in the soils through years of fertilization (Basu et al. 2011). Tile drainage, which is common in Ontario agricultural areas, may contribute to this effect as event water completely bypasses riparian areas. Similarly, urban stormwater in highly impervious catchments also has the tendency to miss the riparian zone as it is carried through stormwater drainage (Hogan and Walbridge 2007) which, along with high groundwater and high stormwater NO3, may cause NO3 chemostasis in urban catchments (Long et al. 2014). Similar to NO3 concentration patterns observed in the present research, a study of intensive agricultural catchments in Manitoba, Canada did not find chemostasis in P concentrations and suggested that the presence of legacy nutrients alone was not sufficient to result in discharge‐invariant nutrient dynamics, particularly in catchments with seasonally variable hydrologic pathways (Ali et al. 2017). The lack of anthropogenic C‐Q responses may also be due to the fact that our catchments, although largely agricultural, were also mixed, containing a mosaic of other land covers (primarily wetland, forest, and urban, with smaller proportions of barren, shrubland, and open water).

We observed below‐threshold near‐chemostasis in many cases for TP, SRP, and TKN; however, the management strategies for nutrients that show this discharge‐invariance across the entire discharge range would require different approaches than those we discussed earlier. For NO3, for example, stream channel restoration aimed at improving denitrification rates would work poorly at high flows because absolute NO3 losses to denitrification would be low compared with the volume of discharge water (Royer et al. 2004, Smith et al. 2006). Chemostasis tends to occur in catchments that have high ambient concentrations (Moatar et al. 2017); catchment management in these cases should therefore focus on terrestrial‐based projects, such as working to reduce legacy nutrient stores and redirecting urban stormflow to riparian wetlands.

Conclusion

The identification of regional seasonal thresholds and slope patterns in C‐Q relationships may be an important first step for water quality managers as it gives an indication of broad, regional patterns which can help to focus mitigation efforts. By examining overall C‐Q patterns inferences can be made about the dominant processes influencing solute concentrations and at what point in the discharge range these tend to switch. We found that TP, SRP, and TKN C‐Q slopes generally increased substantially after a discharge exceedance threshold that varied seasonally suggesting that management decisions and efforts should primarily focus on strategies that reduce peak flows and/or increase nutrient retention within the catchment, such as having healthy riparian wetlands. Our results show that peak flows could be managed based on C‐Q relationship breakpoints which would be beneficial in reducing annual nutrient loads. To minimize large nutrient pulses during extreme events, strategies that increase hydrological connectivity during lower flows and/or reduce nutrient accumulation in the catchment would also be useful for preventing pronounced increases in C‐Q slopes. Catchment‐specific strategies are also needed to account for seasonal effects. These include adjusting the timing of agricultural fertilizer application and urban development activities. NO3 management may need to additionally consider anthropogenic impacts during snowmelt and summer. Considering that the connectedness of the catchments is also important, nutrient control measures may benefit from either reducing the concentrations in solute pools that are only connected to the stream during high‐flow events, or by improving the connection of these pools in the catchment to inhibit solute build‐up between events. In this study we used C‐Q relationships as a diagnostic tool to help identify the dominant processes controlling nutrients in streams, but our method can be easily adapted and applied to any region that already possesses long‐term multisite water quality and stream discharge datasets, giving watershed management agencies the ability to rapidly develop beneficial mitigation strategies. As intensive, high‐resolution monitoring is neither financially nor logistically feasible for these agencies to conduct, the use of regional C‐Q curves to inform management decisions represents a realistic and potentially powerful tool to improve our freshwater systems. While individual catchment processes will be variable across sites and events, identifying the dominant patterns will help watershed managers form regional policies to control nutrients.

Appendix S1

Click here for additional data file.