Modelling historical variability of phosphorus and organic carbon fluxes to the Mackenzie River, Canada

The transport of organic carbon and nutrients from watersheds to receiving water bodies is known to be a hydrologically driven process. Previous studies of the spatial and temporal patterns of fluxes indicate that hydrologic connectivities through the surface and subsurface runoff are the key controllers of mobilization and transport processes, including in the Arctic and subarctic (Stieglitz et al. 2003; Laudon et al. 2011). Specifically, continuous/discontinuous permafrost underlay most of the Arctic and subarctic regions, and landscape controls, such as wetlands, glaciers and their connectivities, strongly influence various geo-chemical fluxes (Frey et al. 2007; Laudon et al. 2011). Hydrochemical connectivity also varies on an annual cycle, with runoff from seasonal snowmelt, frozen soil thaw and summer rainfall creating flow pathways. Hence, higher spring/summer flows dominate the annual cycles of the dissolved and particulate nutrients, as well as organic matter fluxes across river basins in the Arctic/subarctic (Cai et al. 2008; Holmes et al. 2011). However, the spatio-temporal variability of the hydrochemical connectivity leads to highly complex and non-linear relationships between geochemical fluxes and discharge (e.g., Holmes et al. 2011) that could vary between the rising and falling limbs of a hydrograph (e.g., Gareis & Lesack 2017).

The Arctic freshwater system is changing rapidly in response to a warming climate (Bring et al. 2017), with significant implications for a range of hydro-ecological regimes (Prowse et al. 2015; Wrona et al. 2016). Previous studies indicate increasing annual streamflow across the Eurasian and North American river basins due to the enhanced poleward moisture transport (Zhang et al. 2013) and increasing dissolved organic matter (DOC) fluxes from the Mackenzie River (Tank et al. 2016). A key reason behind the changing geochemical fluxes is thawing permafrost, which is exposing substantial quantities of stored organic carbon and macronutrients such as nitrogen and phosphorus, thus making them available for mobilization and transport to aquatic systems (Frey et al. 2007; Schuur et al. 2015; Kendrick et al. 2018). Further, a warmer climate enhances microbial processing and plant uptake of nutrients, leading to changes in nutrient and DOC cycling. For instance, nitrogen mineralization in excess of plant uptake could lead to an increase in supply to surface water (Jiang et al. 2016; Kendrick et al. 2018). The leaching from minerals could increase phosphorus concentration in stream water (Kendrick et al. 2018), while cycling of phosphorus in the headwaters could potentially lead to reductions (Jiang et al. 2016). In the case of DOC, while the thawing permafrost makes it available for transport, a warmer climate could accelerate the microbial breakdown of organic carbon and cause the atmospheric release of carbon dioxide and methane (Schuur et al. 2015), and potentially lead to decrease in DOC fluxes (Kendrick et al. 2018). These changes could affect the riverine transport of geochemical constituents to downstream receiving bodies such as estuaries and the adjacent ocean environments. In particular, Arctic estuaries and near shore coastal environments are vulnerable ecosystems and enhanced nutrient and organic matter transport will have implications on the biogeochemistry and productivity of these systems (Holmes et al. 2011; Casper et al. 2015; Wrona et al. 2016).

In this context, numerical modelling capability is important for establishing baseline conditions and projecting future climate change effects on biogeochemical fluxes. Deterministic modelling of water constituent fluxes for the high-latitude regions is constrained by limited data and lack of understanding of key hydro-geochemical processes. An alternative approach is statistical modelling based on the relationship between causative and resultant variables. To this end, a widely used approach is the regression based Load Estimator (LOADEST) model (Runkel et al. 2004), including applications in modelling historical changes in nutrients and carbon fluxes in the high-latitude regions (e.g., McClelland et al. 2007; Holmes et al. 2011; Tank et al. 2016). For example, McClelland et al. (2007) used LOADEST to estimate nitrate and dissolved organic carbon (DOC) export from the upper Kuparuk River, Alaska, and found that annual nitrate export increased by ∼5 fold and annual DOC export decreased by about one half from 1991 to 2001. Tank et al. (2016) used the LOADEST model to estimate changes in DOC fluxes from the Mackenzie River and its tributaries and found increasing trends.

It is also important to note that the regression-based estimates of water constituent fluxes and trends have inherent uncertainties. In particular, major uncertainties arise from the fact that the statistical models cannot fully capture the relationships between different variables. Additionally, a typical monitoring strategy collects only a limited number of water constituent samples in a given year and, hence, estimates such as a change in concentration can be highly uncertain (Hirsch et al. 2015). Further, sparse and often irregular sampling makes the statistical analysis (such as trend analysis) of water constituent data problematic. This is especially relevant given the limited water constituent observations in most high-latitude regions such as the Arctic Archipelago and eastern shores of Hudson Bay (Li Yung Lung et al. 2018). Therefore, the analysis of uncertainties should be integrated into regression-based estimates of water constituent fluxes.

Further, by evaluating the relationships between variables used in the regression models, avenues for model improvements could be explored. In this respect, concentration–discharge (C–Q) relationships could provide important insights on the presence of solute mass storage within the catchment, and landscape and hydrologic controls that regulate transport to the stream (Goodridge & Melack 2012; Moatar et al. 2017). Further, statistical predictors may be inferable from the nature of relationships (Musolff et al. 2015), thus facilitating the development of an appropriate statistical modeling framework such as LOADEST for estimating loads.

Given these challenges, the main objective of this study is to provide an improved modelling framework for estimating the water constituent fluxes together with uncertainties. This objective was accomplished via analyses of water quantity–constituent relationships and statistical modelling, using the total phosphorus (TP) and DOC dataset for the Hay, Liard and Peel Rivers in subarctic western Canada. We employed C–Q plots to analyse the seasonal, inter-basin and inter-constituent relationships. Based on plausible physical relationships, we explored improvements in the regression estimates of fluxes by (i) using the logarithmic LOADEST and the linear regression models, (ii) considering air-temperature as an alternate input variable and (iii) employing quantile mapping for bias-correction. Further, we analysed uncertainties in load estimation by using a bootstrapping method. Using the improved framework, we evaluated the annual and seasonal variability, and trends in fluxes.

This study was carried out using the datasets from the Hay, Liard and Peel Rivers, which are major western tributaries of the Mackenzie River (Figure 1). These river basins have low population densities and are mostly in pristine states, but oil and gas exploration and development, and forest harvesting activities are present in parts of the Hay and Liard basins (Mackenzie River Basin Board 2003). As a result of these anthropogenically driven activities, including climate-driven changes in the hydrologic regime, there is increasing recognition of ecological implications such as contaminant and nutrient enrichment and transport in these river systems (Wrona et al. 2016).

The three basins are located in the Boreal and Taiga ecoregions in the subarctic, with Liard the largest and Hay the smallest of the three basins (Table 1). Based on the permafrost extent map from Heginbottom et al. (1995), continuous (92%) to extensive discontinuous (8%) permafrost underlay the Peel basin, while permafrost is less prevalent in the Liard (continuous 1%, extensive discontinuous 34%, sporadic discontinuous 57%, isolated patches 8%) and Hay (sporadic discontinuous 75%, isolated patches 17%) basins. The Liard and Peel basins are physiographically heterogeneous with more than 50% of the basin area located in the Cordillera Mountains and the remainder in the interior plains. The mountains in these basins receive proportionally more precipitation and produce more runoff than the downstream plains (Slaymaker & Kovanen 2017). The entire Hay basin lies in the interior plains. Vegetation cover is diverse in the three basins. Based on the classification by Latifovic et al. (2010), needleleaf forest (27%), mixed forest (17%), shurbland (17%) and mixed forest (16%) are the main land-cover types in the Hay River basin. In the Liard basin, needleleaf forest (51%) is the dominant land-cover class with mixed forest (16%), shrubland (10%) and barren land (9%) as other main land-cover types. The Peel basin has grassland (35%), shrubland (24%), barren land (20%) and mixed forest (7%) as the main land-cover types.

The climate system in the region reflects the complex interactions between the large-scale forcings and the physical environment, with topography exerting a major influence on atmospheric circulation, especially over the Cordillera Mountains causing strong precipitation gradient (Szeto et al. 2008). The basins encompass a wide variety of climatic conditions, including the cold temperate, mountain and subarctic zones and characterized by short cool summers and long cold winters (St Jacques & Sauchyn 2009). The runoff regime of all three basins is subarctic nival, with the bulk of runoff produced by snowmelt during late spring/early summer (Woo & Thorne 2003). Overall, these three basins cover about 22% of the drainage area of the Mackenzie River while contributing about 33% of discharge (based on the Water Survey of Canada Hydrometric Station data, described later).

Recent studies suggest hydro-climatic regimes are changing in these basins. For example, St Jacques and Sauchyn (2009) found a statistically significant increase (at 10% significant level; p < 0.1) in winter baseflow but non-significant annual change for the Hay (1964–2007), Liard (1973–2007) and Peel (1975–2007) basins. Rood et al. (2016) used extended records to 2014 and found non-significant changes in the mean annual discharge of the Peel and Hay Rivers and significant increase (p < 0.05) in the mean annual discharge for an upstream Liard station (at Fort Liard). In the case of water constituent loadings, Tank et al. (2016) found significant increasing trends (p < 0.05) in the LOADEST simulated annual DOC, alkalinity and sulphate fluxes from the Liard basin and non-significant trends in these variables for the Hay basin. Analyses of trends and seasonal variability of TP fluxes for the Liard and Hay basins, as well DOC and TP trends for the Peel basin remain to be completed.

We used water constituent data near the outlets to evaluate the aggregated DOC and TP fluxes from the Hay, Liard and Peel basins (Figure 1). The datasets were obtained from the Water Quality Monitoring and Surveillance branch of Environment and Climate Change Canada (https://www.canada.ca/en/environment-climate-change/services/freshwater-quality-monitoring/online-data.html). We considered TP and DOC because the focus of this study is nutrient and carbon fluxes, and long-term records for only these two constituents are available. The periods of records used correspond to the length of all available records, which for TP (DOC) consisted of the periods 1988–2012 (1988–2012), 1972–2012 (1979–2012) and 1979–2012 (1979–2012), respectively, for the Hay, Liard and Peel Rivers (see Supplementary Materials, Table S1 for details, available with the online version of this paper). Sampling frequencies of these datasets are about six points in a year and, hence, there are limitations in evaluating historical variability and trends using these datasets.

Streamflow data for the three basins nearest to the freshwater quality stations (Figure 1) were obtained from the Water Survey of Canada online archive (http://wateroffice.ec.gc.ca/) (Table 1). Historical annual and seasonal streamflow trends were obtained by using the non-parametric Mann–Kendall test (Kendall 1955) with trend-free pre-whitening employed to remove the effects of serial correlation (Yue et al. 2002). The runoff trends matching the periods of available water constituent records are mostly non-significant at the 5% significance level (p < 0.05), except for winter flow increases for the Liard and winter and autumn increases for the Peel (Table 1).

Additionally, basin average daily air-temperatures were extracted from the gridded 1/16° spatial resolution Pacific Northwest North-American meteorological (PNWNAmet) dataset (Werner et al. 2019). Based on the datasets, the mean annual temperatures of the Hay, Liard and Peel basins are −0.8 °C (1988–2012), −2.2 °C (1979–2012) and −5.9 °C (1979–2012), respectively (see Table S2, available online). Mean annual temperatures for all three basins are increasing, which is statistically significant (p < 0.05) in the case of Peel and Liard basins.

Before developing a statistical water constituent model, it is important to understand the relationships between driving and resultant variables. To this end, we first analysed correlations of TP and DOC concentrations with river discharge (Q) and mean daily air-temperature (T_mean). We considered the correlation with T_mean as it is a useful indicator of seasonal and interannual variability. Secondly, we analysed the basin-scale hydrologic and biogeochemical relationships by using the C–Q plots.

However, instead of pre-selecting the most appropriate regression model, we applied both linear and logarithmic regression models in all cases, and then selected the best performing model. Specifically, a single MLR model and nine LOADEST logarithmic regression models were fitted using the entire dataset, with the adjusted maximum-likelihood estimation method used for estimating the regression parameters. The best performing model out of the 10 regression models was selected using the Akaike information criterion (AIC) (Akaike 1974), which is the standard approach in LOADEST (see Runkel et al. 2004), and previous studies (e.g., Tank et al. 2016; Li Yung Lung et al. 2018) also used the same approach. AIC provides a robust method against overfitting by trading goodness of fit with model complexity by including a penalty term against increasing numbers of parameters (Akaike 1974). Thus, the criteria lead to the selection of a regression model with parsimonious parameterization, and in contrast to calibration–validation approach allows the model selection using the entire dataset. Model selection using the AIC was repeated separately for all variables (TP and DOC) and basins (Hay, Liard and Peel).

Second enhancement in the load estimation model consisted of using an alternate input variable. The ‘Default Method’ of LOADEST uses streamflow and decimal time (Table 2). While time provides a representation of the seasonality variability, it cannot represent the interannual variability as dates change uniformly for all years. Although T_mean does not directly influence the constituent concentration and fluxes, it can represent both seasonal and interannual variability. Additionally, a previous study found improvement in the model performance when they used T_mean as an explanatory variable in the fuzzy-rule based model for estimating nitrate fluxes (Shrestha et al. 2007). Thus, given the prevalence of the relationship of concentration with air-temperature (discussed later), we considered T_mean as a replacement for decimal time (herein referred to as Method 1).

Third enhancement in this application pertains to the way the model handles bias. LOADEST has been known to produce biased load estimates (Runkel 2013) and uses a constant correction factor (delta factor) for adjusting the biased load (Default Method). However, a closer examination of the LOADEST model biases revealed that the model usually underpredicts the higher fluxes and overpredicts the lower fluxes (not shown). A constant delta-factor (either addition or subtraction) cannot correct both underprediction and overprediction bias at the same time. On the other hand, quantile mapping, which is a widely used method for statistical downscaling of climate models (e.g., Maurer & Hidalgo 2008), performs bias-correction by transforming the cumulative distribution function of the simulations to match that of observations. Hence, quantile mapping can handle both underprediction and underprediction type biases together and employed in this study. We applied quantile mapping to both Default Method and Method 1 by replacing the delta factor approach in the LOADEST and MLR models, which are herein referred to Methods 2 and 3, respectively (Table 2). Quantile mapping was performed by using the R ‘qmap’ package (Gudmundsson 2016).

We compared the statistical performance of the fitted models using the Nash–Sutcliffe coefficient of efficiency (NSE) (Nash & Sutcliffe 1970) and residual bias, and selected the best performing models from the combinations of four methods. These metrics were calculated on non-logarithmic loads. The NSE can range from −∞ to 1, with higher values (closer to 1) indicating better agreement. The residual bias refers to the cumulative % error after bias-correction using either delta factor or quantile mapping method. The lower value of residual bias (closer to 0) indicates better agreement.

As previously stated, regression-based load estimates are affected by the uncertainty that arises from the fact that the model is not able to fully explain the variability of observations. Further, uncertainties arise from errors in the observed data. Hence, a resampling method was implemented for the selected best performing method, and uncertainties in the MLR and LOADEST fluxes were estimated. We used the non-parametric residual bootstrapping method (Kreiss 1997), which generates ‘new’ samples of ‘observed’ data from the original samples. This consisted of (i) fitting a set of regression parameters and quantile mapping to the observed samples, and then calculating the residuals as fractions ((observed−simulated)/observed), (ii) resampling the residuals with replacements, then multiplying the output from step (i) with 1-replacements to generate ‘new’ samples and (iii) fitting the regression parameters to the ‘new’ samples, and then performing quantile mapping. In this case, we used 10,000 bootstrap iterations to generate a time series of fluxes.

We used the fitted model and bootstrapping results to evaluate the annual and seasonal variability and trends in fluxes for the periods of records under consideration. In this case, also, trend analysis was conducted by employing the non-parametric Mann–Kendall test (Kendall 1955) with the trend-free pre-whitening (Yue et al. 2002) employed to remove the effects of serial correlation. Trends were considered statistically significant at p < 0.05.

The TP and DOC concentrations show statistically significant (p < 0.05) correlations with Q and T_mean providing evidence of the underlying relations (Table 3). The correlations for TP versus Q are especially strong; thus, discharge alone can explain about 50% or more variance in TP. An exception is the Hay River DOC concentration, which is uncorrelated with either Q or T_mean, suggesting changes in these variables do not directly affect the concentration.

Further, the log–log plots showed strong relationships (Figure 2). Specifically, they exhibit increasing concentration with increasing discharge (chemodynamic relationship) that is in agreement with a strong positive correlation. Additionally, the relationships for the winter low-flow (air-temperature <0 °C) and summer high-flow (June–August) periods showed clustering of low concentrations with low flows and high concentrations with high flows. Such a positive or ‘up’ chemodyamic relationship has been described as ‘transport limited’ (Musolff et al. 2015; Moatar et al. 2017), i.e., enhanced hydrologic connectivity established during high flows reconnect and convey the constituents from the areas where a given element is more abundant. The slope of the C–Q relationship is steeper for TP compared with DOC, suggesting that the connectivity is especially important in the transport of phosphorus. The exception to the chemodynamic behaviour is the Hay River DOC–Q plot, for which there is no change in concentration with increasing discharge (chemostatic). Such chemostatic behaviour has been attributed to the homogeneous and uniform distribution of elements in the catchment, meaning that changes in hydrologic connectivity and flow path do not affect the solute concentration (Duncan et al. 2017; Moatar et al. 2017). Further, Creed et al. (2015) found chemostatic DOC–Q behaviour across rivers in the USA, and attributed to a shift from headwater-dominated transport processes to downstream-dominated instream and near-stream biogeochemical processes. The relationships also exhibit considerable spreads indicating complex interactions with flows such as the concentration effect of decreasing flow and dilution effect of increasing flow (Shrestha et al. 2013).

It is interesting to note that while TP concentrations for three basins are in similar ranges, the DOC concentration decreases progressively from the Hay to the Liard to the Peel basin (Figure 2). The proportion of the vegetation covered area also decreases progressively across the three basins, and is likely responsible for the differences in DOC across the basins. The TP concentrations are low in all three basins, and their similar ranges in all three suggest vegetation cover is not an important factor in the generation of TP. Hence, it is important to understand different forms and the dominant transport mechanism of these constituents. Based on the observed time series of TP and particulate P, TP mostly occurs in particulate form in all three basins (average Particulate P to TP ratios for the Peel, Liard and Hay Rivers are 0.76, 0.74 and 0.60, respectively). Runoff from snowmelt and surface runoff are known to be the main drivers of particulate TP transport (Cai et al. 2008; Heckrath et al. 2008). By contrast, the subsurface flow has been found to be the main driver of DOC transport and, hence, factors such as frozen ground and permafrost thaw could play important roles in the DOC cycling (Striegl et al. 2005; Walvoord & Striegl 2007).

The first step in regression modelling of TP and DOC loads is selecting the model with the best overall performance from an MLR model and nine LOADEST models. Based on the AIC criteria, the MLR models were selected for estimating DOC fluxes from the Hay River for all four methods summarized in Table 2. For all other combinations and rivers, different forms of LOADEST logarithmic models were selected. This reinforces the fact that the use of logarithmic LOADEST equation is only suitable when there is a logarithmic relationship between variables. Hence, there is a need to carefully examine the relationship between variables before applying regression models for estimating fluxes.

Table 4 summarizes the statistical performance of the best performing models for the four methods. The results for the Default Method and Method 1 are generally similar. Only a small improvement in the NSE was obtained when the decimal time was replaced by T_mean (Default Method vs. Method 1). However, with the application of delta method for bias-correction, large residual biases remained in both the methods. Replacing the delta bias-correction with quantile mapping (Methods 2 and 3) led to considerable improvement, especially in the case of residual % bias. This also led to improvements in the NSE values. Thus, quantile mapping provides an effective approach for improving the overall statistical performance of the regression model. Further, there are small improvements in the model performance when comparing Methods 2 and 3. This is likely due to the fact that T_mean (used in Method 3) can represent interannual variability in contrast to decimal time (used in Method 2), which varies uniformly for all years. The representation of interannual variability could be useful for continuous daily simulations; thus, it was considered more appropriate to use T_mean as an input variable.

Further, Figures 3 and 4 show comparisons of the observed and simulated seasonal and interannual variability of the TP and DOC fluxes using Method 3, with the results compared for the days when observed data are available. Corresponding comparisons for Method 1 are shown in Figures S1 and S2 (available with the online version of this paper). The results clearly reveal improvements with the quantile mapping approach, especially for the low flux values (both monthly and annual results), for which Method 1 has a tendency to overpredict (Figures S1 and S2) and Method 3 provides a better representation (Figures 3 and 4). As explained previously, the delta method adjusts biases using a constant value for the entire dataset, while quantile mapping, by definition, transforms the cumulative distribution function of simulations to match that of observations. Overall, the results in Figures 3 and 4 indicate a good ability of the model using quantile mapping to replicate the monthly dynamics, i.e., lower fluxes from November to April and higher fluxes from May to October, as well as interannual variability.

More importantly, given that one of the objectives of this study is to analyse historical trends, we considered the ability of the models to simulate trends. Specifically, the blue and red lines compare trends for the annual medians of the observed and simulated data points, respectively (lower panels Figures 3 and 4; Figures S1 and S2). The results show that despite some differences in the magnitude of trends, the directions of trends (mostly increasing) are correctly replicated by both Methods 1 and 3. Hence, this gives us some confidence in using the entire set of model results for the evaluation of the historical trends. Overall, given the improved statistical performance and a better representation of variability, Method 3 clearly provides an improved framework for evaluating historical baseline and trends and, thus, selected as the preferred model. Further, discussions are based on this improved model (Method 3), i.e., Q and T_mean as input variables and quantile mapping bias-correction (see Table S3 for details, available online).

It is important to note here that the sole purpose of these trend lines plots (Figures 3 and 4; Figures S1 and S2) is to assess the model's ability to replicate trends. As the number of observations is too few (typically about six samples in a year), these lines are totally dependent on whether the samples hit low, medium or high flux values and, thus, do not represent historical trends. The discussion below will shed light on the lack of representativeness of these trend lines.

After selecting Method 3 for further analysis and performing bootstrapping on this combination, we estimated the annual and seasonal variability and trends in TP and DOC fluxes, as well as uncertainties (Figures 5 and 6). All flux values are normalized to kg/day, and the ranges reflect uncertainties due to regression-based load estimates. The results show large uncertainty ranges, especially for the higher flux seasons of spring and summer compared to the lower flux seasons of autumn and winter (the fluxes are plotted on log-axis). However, the uncertainty ranges follow the variability of fluxes obtained from the selected models (red dots), thus, indicating general robustness of the models to represent variability.

Considering the seasonal variability of TP fluxes, all three basins exhibit lower fluxes in winter (DJF) and higher fluxes in spring (MAM) and summer (JJA) (Figure 5). The lower winter fluxes correspond to the lower winter runoff values for the three basins (Table 1), but the spring and summer responses vary considerably. For example, spring to summer runoff ratios are 0.7, 0.3 and 0.5 for the Hay, Liard and Peel basins, respectively, while the respective flux ratios for TP are 1.4, 0.2 and 0.9. Given that TP mostly occurs in a particulate form in these basins, the higher TP flux ratios compared to runoff ratios imply a larger contribution of surface runoff-driven TP transport during the rising limb of the hydrograph. This also implies that TP peaks precede discharge peaks, as has been observed in other cold-regions basins (McNamara et al. 2008; Shrestha et al. 2012). The Liard is a much larger basin than the Hay and Peel, with longer travel time for runoff to reach the basin outlet. Although, groundwater input is the dominant and most persistent streamflow source in the Liard basin, isotope analysis by St. Amour et al. (2005) found some contribution of surface runoff during the post-freshet period. Additionally, parts of the lower Liard valley is wetland dominated, where runoff is slowly conveyed to the basin outlet through a series of minerotrophic channel fens (Connon et al. 2014). Thus, the runoff from nutrient-rich lowland areas requires a longer time to reach the basin outlet. These factors likely contributed to the higher fluxes during summer months and, hence, higher summer to spring flux ratios. In the case of autumn, the fluxes are lower as runoff decreases (Table 1) in all three basins.

In the case of long-term trends, the results revealed small non-significant annual increases in TP for the Hay and Liard basins and a non-significant annual decline for the Peel basin (lower panel, Figure 5). While the increasing annual trends match the direction of indicative trend lines for the Hay and Liard basins (Figure 3), the directions are opposite for the Peel basin. This discrepancy stems from the limited number of sampling data available for plotting trend lines in Figure 3, thus, reinforcing the fact that trend analysis using limited water constituent monitoring data is problematic. Seasonally, the results revealed increasing trends for the winter months (Figure 5), with statistically significant increases for the Liard and Peel basins, which correspond to the significant streamflow trends (Table 1). However, given the smaller fluxes in winter, the contribution of this winter increase on annual fluxes is small. The trends for spring and summer months are small and non-significant except for the Peel basin, which showed a significant decline in the summer months. Autumn changes are affected by warmer temperature-driven increases in rainfall fractions and, thus, runoff increases (Table 1) with increasing trends for all three basins (significant only for the Liard basin). Overall, there is a match in the direction of TP flux trends and runoff trends, except for the Hay basin spring results. The results also show large uncertainties for some of the trends, which is related to the errors in load estimates. However, the uncertainty ranges for all statistically significant trends are in the same direction as the selected model, indicating the robustness of the trends.

The DOC flux results are largely similar to the TP fluxes, i.e., lower winter fluxes, higher spring and summer fluxes, and intermediate autumn fluxes (Figure 6). A notable difference compared with TP is a close agreement of the spring to summer DOC flux ratios with runoff ratios. Specifically, the runoff ratios are 0.7, 0.3 and 0.5 for the Hay, Liard and Peel basins, respectively, while the respective ratios for the DOC fluxes are 0.7, 0.4 and 0.7. This suggests a dominance of groundwater flow in all three basins, as has been found to be the case for the Liard basin (St Amour et al. 2005), which is the primary transport mechanism for DOC. Given that the bootstrap uncertainty ranges for the DOC fluxes are smaller than those for the TP fluxes, there is higher confidence in the regression model estimates of the DOC fluxes. In the case of trends, the DOC fluxes show non-significant increasing trends for the Hay and Liard basins and a small non-significant decrease for the Peel basin (Figure 6). Again, the Peel DOC trend differs from the indicative trend of the median fluxes (Figure 4), further illustrating the inadequacy of the limited sampling data to represent an annual trend. The seasonal trends mostly follow the streamflow trends (Table 1), i.e., a larger increasing trend in winter, and smaller increasing or decreasing trends for spring, summer and autumn. In this case, only the increasing trends for the Peel winter and autumn fluxes are statistically significant. The uncertainty ranges of the DOC trends are primarily smaller than the TP trends (except for the Liard winter trends), thus, implying higher confidence in the regression-based DOC trend estimates.

Comparing the results of this study with previous studies, there is a general agreement in the direction of trends, while magnitude and statistical significance of trends differ for some basins. For instance, the non-significant annual trend in DOC fluxes for the Hay basin (this study) agrees with that in Tank et al. (2016). In the case of the Liard basin annual DOC flux trend, while this study found a non-significant increase, Tank et al. (2016) found a significantly increasing trend (p < 0.05). The differences probably arise from the use of different bias-correction methods and inputs in the two studies. Further, across the region, previous studies indicate significantly increasing (p < 0.05) DOC fluxes in the Mackenzie River for the years 1972–2012 (Tank et al. 2016). The trends in the Yukon DOC fluxes indicate decreases for the years 1976–2003 (Striegl et al. 2005) and no trend for 2001–2014 (Toohey et al. 2016). Thus, there is no consensus in the direction of DOC changes. This may be due to the fact that the increased availability of organic matter from thawing frozen soil and permafrost could be offset by the increase in CO₂ sequestration or atmospheric release (Tank et al. 2012). In the case of phosphorus, there are no previous estimates of trends in the studied basins, while increasing the direction of the trend for the Yukon River (Toohey et al. 2016) match with the Hay and Liard basins in this study.

Overall, the variability of the magnitudes and trends of fluxes largely follow those of basin runoff. The impacts of a warming climate could be seen in the significant increases in winter streamflow and water constituent fluxes. While the historical annual and seasonal trends are mostly small and non-significant, it is interesting that both the TP and DOC trends are amplified compared with runoff trends. Specifically, compared with the annual runoff trends, which vary between −0.9% and +3.2% change per decade, the TP and DOC trends vary between −11% to +5% and −3% to +10% per decade, respectively, with the direction of the TP and DOC trends matching the direction of runoff trends. Similar results are also obtained for the seasonal trend values, i.e., the TP and DOC trend values are some order of magnitude greater than the runoff trends. This suggests the possibility of greater changes in the water constituent fluxes if significant changes in the streamflow were to occur in the future climate. This is especially relevant given the projected enhanced moisture transport and amplified hydrologic response in the Arctic/subarctic. For instance, for the Liard basin under future climatic conditions, Thorne & Woo (2011) and Shrestha et al. (2019) projected the occurrence of an earlier spring freshet, increased autumn to spring discharge, and an overall increase in annual discharge. Given the amplified trends in historical water constituent fluxes, such future changes in discharge could bring more pronounced change in the water constituent fluxes.

This study provides an assessment of historical variability and trends of TP and DOC fluxes from key tributary basins to the Mackenzie River using an improved statistical modelling framework. We evaluated the relationships between TP and DOC with river discharges for the Hay, Liard and Peel River basins. The analysis revealed primarily chemodynamic relationships, i.e., increasing concentration with increasing discharge except the chemostatic DOC–Q relationship for the Hay basin. The chemodynamic relationship suggests that the fluxes are transport limited and the enhanced hydrologic connectivity during high flows reconnects and conveys water constituents from areas where they are most abundant. The chemostatic DOC–Q relationship for the Hay basin suggests homogeneous and uniform distribution of elements in the catchment, and/or dominance of downstream riverine and near-stream processes.

Based on these insights, we implemented enhancements in the regression-based load estimation first by using the multiple linear regression (MLR) in addition to the USGS Load Estimation (LOADEST) models. Secondly, we replaced decimal time with air-temperature for a better representation of the interannual variability. Thirdly, we replaced the delta method for bias-correction with a quantile mapping method. Model selections based on the LOADEST logarithmic models and MLR models led to the logarithmic models when the relationships between constituent concentration and discharge are chemodynamic, and the MLR model when the relationship is chemostatic. Therefore, there is a need to carefully examine the relationships between variables before applying regression models for estimating fluxes. Further enhancements in this study, especially the use of quantile mapping led to improvements in the model performances and provided an improved framework for evaluating historical variability and trends. Although the improvements using air-temperature as an input variable are small, we decided to proceed with it because of its ability to represent interannual variability. Additionally, there is a need to evaluate uncertainties in load estimates such as the bootstrapping method implemented in this study.

The results revealed that variability of TP and DOC fluxes primarily follows those of basin runoff, reinforcing the prevailing understanding that water constituents are hydrologically driven. Additionally, as indicated by the differences in the TP and DOC responses, their forms (dissolved or particulate), flow paths and travel times all affect the transport processes. The trends of the TP and DOC fluxes are primarily small and non-significant, and generally follow the direction of runoff trends, but the magnitudes of the trends are amplified. The largest significant increases in all three basins occur in the winter months. Annually, the results revealed small non-significant increases in TP and DOC for the Hay and Liard basins, and a non-significant decline for the Peel basin. The annual trends for the Peel basin are in an opposite direction to the indicative trends derived from the available observed data. This highlights the inadequacy of using limited water constituent monitoring data, which are about six points in a year, directly for trend analysis. The sparse data also lead to limitations in modelling water constituent fluxes. Hence, for improved quantification of variability and trends, the higher temporal resolution of water constituent sampling should be prioritized. This will also facilitate an improved modelling of the historical water constituent fluxes as well as a projection of the future climate conditions.

Overall, the primarily small and non-significant trends in TP and DOC fluxes for the historical period indicate limited effects of a warming climate for the period under consideration. However, the trends in constituent fluxes are amplified compared with discharge suggesting the possibility of more pronounced changes if larger changes in streamflow were to occur. This leads us to a subsequent question about future changes to water constituent fluxes in these basins. While future hydrologic projections suggest an increased runoff and earlier timing of snowmelt-driven runoff as a result of enhanced moisture transport in the high-latitude regions (Thorne & Woo 2011; Bring et al. 2017), the implications of such hydrologic changes on water constituent fluxes remain to be addressed. Of particular importance is that thawing permafrost and enhanced microbial activity in a warmer climate will make trapped nutrients and organic matter more hydrologically available, which could have significant implications for water constituent fluxes. The evaluations of such future changes should be the focus for future research, and the improved insights about historical variability and trends obtained from this study would provide a useful baseline for such assessments.

We acknowledge the Water Quality Monitoring and Surveillance branch of Environment and Climate Change Canada for providing water constituent data used in this study. We thank two anonymous reviewers for their comments, which led to an improved version of the manuscript.