Nutrient inputs from subarctic rivers into Hudson Bay

Arctic regions and adjacent sub-Arctic areas have been altered rapidly by ongoing climate change, with noticeable increases in surface air temperatures, precipitation, and river discharge, as well as reductions in permafrost and the extent and thickness of sea ice and its associated snow cover (e.g., Chen et al., 2015; Dodd et al., 2015; Gao et al., 2015; Vonk et al., 2015; Déry et al., 2016). These changes are expected to have particularly severe impacts on the subarctic domain, where steep environmental gradients mark the transition between temperate and polar climates. Hudson Bay (HB), a wide and shallow inland sea that connects the Arctic and temperate regions of central Canada, can be regarded as an archetype and sentinel area for the study of climate change impacts on coastal waters in this transition zone (Macdonald and Kuzyk, 2011). By virtue of their geographical location, as well as considerable exposure to terrestrial inputs of freshwater and matter, sub-Arctic estuarine shelf systems such as HB are regarded as the most sensitive to climate change (Mackenzie et al., 2004; Liu et al., 2010).

Two of the largest rivers that discharge into HB, the Nelson River and the La Grande River, are regulated for hydropower production but also for flood protection, irrigation, and industrial and recreational purposes, and thus are considered to be moderately to strongly fragmented (Dynesius and Nilsson, 1994). The development of hydroelectric installations began in the headwaters of these two river systems in the early 20^th century, culminating in two river complexes which strongly alter flow regimes in the present day (Grill et al., 2015). The Nelson River is the largest and discharges to the southwest of HB, while the La Grande River discharges to the southeast in James Bay and became the second largest river (by discharge) draining into the bay after diversion from three contributing watersheds (Rupert, Eastmain, and Caniapiscau) for hydroelectric development. Previous studies of HB rivers focused mainly on how their discharge responds to climate change and anthropogenic disturbance (Déry et al., 2018). The results show that hydrological regimes have been altered in both the regulated and natural systems, with a projected increasing trend in bay-wide discharge under future climate scenarios (Déry et al., 2016; MacDonald et al., 2018). While these studies improved our understanding of river discharge into HB, the respective impacts of climate change and river regulation on the delivery of chemical constituents to the coastal zone and their consequences for marine biogeochemistry in HB remain to be addressed. A recent study proposed that, worldwide, damming leads to the preferential removal of P and Si over N in reservoirs, which, by increasing the N:P and N:Si ratios of waters delivered to the ocean, raises the potential for P limitation and possibly has a negative impact on diatom productivity in estuaries (Maavara et al., 2020).

Given their supporting role for biological production in coastal Arctic waters (McClelland et al., 2012), river nutrients, their seaward transport and their impact on shelf biogeochemistry were the focus of several investigations (Holmes et al., 2012; Tank et al., 2012; Le Fouest et al., 2013; Tremblay et al., 2014). These studies showed that the biogeochemistry of Arctic rivers is sensitive to temperature, precipitation, discharge, and permafrost thawing and that changes in these factors could alter the magnitude and variability of riverine nutrient deliveries as the climate continues to warm (Spencer et al., 2015; Vonk et al., 2015; Bowring et al., 2020; Juhls et al., 2020). For temperate and tropical rivers, the observed covariation between nutrient concentration and discharge is very useful for assessing the impact of changing river flow on nutrient transports under future climate scenarios (Sigleo and Frick, 2007; Masotti et al., 2018). Whether this covariation applies to HB or not remains to be established.

The aim of this study was to provide a comprehensive, contemporary portrait of inorganic nutrient concentrations in rivers around HB, to estimate the associated transport of these nutrients to the coastal zone, and to assess if rivers that are regulated for the production of electricity differ from unregulated ones in these aspects. To meet these goals, we combined discharge estimates generated by three global climate models with a database of new and historical nutrient data to estimate nutrient transports.

River sampling was carried out between 2016 and 2019, as part of the field campaigns of the Hudson Bay System Study (BaySys; Barber, 2014) and Coastal Oceanography of the eastern James Bay (COast-JB; Figure 1 and Table S1). Water samples were collected into acid-cleaned 15-mL polyethylene tubes after filtration through a GF/F filter inserted in a filter holder to remove large particles and either measured fresh onboard the CCGS Amundsen or immediately frozen upon collection and stored at −80°C until analysis in the laboratory at Université Laval and Université du Québec à Rimouski, Canada. Colorimetric determinations of nitrate plus nitrite (referred to as NO_x), phosphate (PO₄³⁻), and silicate (Si(OH)₄) were performed with an Autoanalyzer 3 (Bran and Luebbe) with routine methods adapted from Hansen and Koroleff (1999). Analytical detection limits were 0.03 μM, 0.05 μM, and 0.1 μM for NO_x, PO₄³⁻, and Si(OH)₄, respectively. Concentrations of ammonium (NH₄⁺) and total dissolved nitrogen (TDN) were measured during the BaySys expedition only. The former was measured using the sensitive method of Holmes et al. (1999) with a detection limit of 0.02 μM, and the latter was measured by high temperature combustion (HTC) using a total organic carbon analyzer (TOC-V_CPN, Shimadzu) with a chemiluminescent nitrogen detector (TNM-1 module) and a precision of 2%. Concentrations of dissolved organic nitrogen (DON) were calculated as the difference between TDN and the sum of NO_x and NH₄⁺. Concentrations below the detection limit were reported here as the designated value for a given nutrient.

Historical data sets for selected rivers were included to assess seasonal and regional patterns in nutrient concentrations and possible relationships with discharge from: the Coordinated Aquatic Monitoring Program (CAMP) 2008–2013; the environmental field study and monitoring for the Conawapa Generating Station (Conawapa GS) in the Nelson River 2004–2009; and previous publications (Hudon et al., 1996; Kuzyk et al., 2010; Tables S2 and S3). The CAMP project sampled widely in Manitoba, but here we worked only with data obtained at stations located in the lowest reaches of the Churchill (150 km from the coast), Hayes (80 km from the coast), and Nelson (100 km from the coast) rivers. For the Conawapa GS project, data originating from different sampling stations distributed between the mouth of the Nelson River and a distance of 100 km upstream were considered. Although the sampling frequency differed between years, the overall data sets contain data from all seasons. The CAMP and Conawapa GS data sets included the concentrations of NO_x, PO₄³⁻, Si(OH)₄, ammonia (NH₃), total nitrogen (TN), and total phosphorus (TP). Kuzyk et al. (2010) measured nutrient concentrations using similar analytical procedures, whereas Hudon et al. (1996) used different methods with detection limits of 0.07 μM for NO_x and 0.03 μM for PO₄³⁻, respectively (references therein). Because the lack of temporal and spatial overlap between these data sets precludes a rigorous comparative analysis, all data were taken at face value given the absence of apparent offsets in the ranges of nutrient concentrations and ratios reported by different investigators. More information on data collection and sample processing is addressed in Hudon et al. (1996) and Kuzyk et al. (2010).

Daily streamflow data for subarctic rivers from 2006 to 2015 were simulated by the Hydrologic Predictions for the Environment model (HYPE; Lindström et al., 2010). Hydrologic simulations were forced using daily precipitation and temperature outputs from the fifth Coupled Model Intercomparison Project (CMIP-5), bias-corrected for the pan-Arctic domain (Braun et al., 2021) to the 2^nd generation Hydrologic Global Forcing Dataset (HydroGFDv2; Berg et al., 2018). The atmospheric forcing datasets used to drive the hydrologic model were the Geophysical Fluid Dynamics Laboratory Climate model (GFDL-CM3), the Model for Interdisciplinary Research On Climate (MIROC5), and the Meteorological Research Institute Coupled atmosphere-ocean GCM (MRI-CGCM3). CMIP-5 models are widely used to reconstruct or predict climate variables (Taylor et al., 2012). HYPE simulations were validated against hydrologic gauge data analyses through a calibration process with added model processes designed to improve characterization of natural (Stadnyk et al., 2020) and anthropogenic (Tefs et al., 2021) water storage and routing. Further details on the selection criteria of the climate simulation and configuration of the HYPE model for the Hudson Bay Drainage Basin (HBDB) are provided in Stadnyk et al. (2019).

The study domain here was defined as the area bounded by Wager Bay in the northwest (65°N, 90°W), the top-left corner of Quebec province in the northeast (63°N, 79°W) and the southern tip of James Bay in the south (51°N). Figure 1 represents all of the subarctic rivers discharging into HB. These rivers were grouped into four different sectors (Northwest, Southwest, South, and East) to allow for a regional comparison across Canadian territories or provinces (Figure 1). This categorization was somewhat arbitrary, but it captures broad differences in nutrient concentrations and fluxes.

Given the seasonality of river discharge and primary production in estuarine and coastal ecosystems in general, we first divided the datasets into four distinct periods according to the subarctic nival regime in river flow (Déry et al., 2005; Woo et al., 2008): a winter period of low flow from November through April, the spring freshet from May to June, a period of intermediate flow during summer (July and August), and a period of secondary flow peaks during fall (September and October). Nutrient concentration data from the aforementioned projects and the literature were binned into these four “seasons” regardless of measurement year. NO_x and PO₄³⁻ concentrations below detection limits were replaced with 0.03 and 0.05 μM before the binning process (Tables S1 and S2). These bins were used to analyze inter-seasonal and regional variability in nutrient concentrations and ratios (hereafter referred to as observation-based data), noting that consistent multi-annual trends were not detected in the rivers for which several years of data were available in a given season (e.g., CAMP).

For each of the 158 rivers draining into HB, daily simulated discharge (in m³ s⁻¹) from three different model outputs (GFDL-CM3, MIROC5 and MRI-CGCM3) was integrated separately for each of the above seasons and then averaged over 10 years. A separate “climatology” of river flow rate for each season (total km³ for the season) was thus created by taking an ensemble mean of the three model outputs. Total annual discharge was calculated by the sum of all season values (in km³ yr⁻¹).

For the nine rivers in which sampling was conducted at least in two separate seasons (i.e., the Churchill, Nelson, Hayes, Severn, Winisk, Maquatua, Castor, La Grande, and Grande rivière de la Baleine rivers, which collectively account for 49% for total annual discharge into HB), all of the data available for a given season (sometimes represented by a single sampling event; Tables S1–S3) were averaged. Missing values for the other seasons were filled by linear interpolation (Table S4), as this approach is used widely to estimate missing data and increase resolution in coarse budgets (Le Fouest et al., 2013). To estimate seasonal nutrient fluxes, seasonal average concentrations were multiplied by the corresponding integrated discharge (10⁶–10⁸ g for each nutrient in a given season). Values from all seasons were summed to obtain annual nutrient transport estimates (10⁶–10⁸ g yr⁻¹ for a given nutrient).

For the other rivers in which seasonal data coverage was inadequate (one season only) or no data were available (149 rivers accounting for 51% of total annual discharge into HB), seasonal nutrient fluxes were estimated using two different approaches. In Method 1, the average value of all rivers for which samples existed in a given season was used as a proxy (Table S5) and multiplied by the corresponding modelled discharge for that season. A refinement was attempted for Method 2 by taking regional differences into account for the two regions where seasonally resolved data existed (the Churchill, Nelson, and Hayes rivers in the Southwest and the Maquatua, Castor, and La Grande rivers, and Grande rivière de la Baleine in the East; Tables S4 and S5). There, regionally averaged concentrations in each season (based on the aforementioned 7 rivers) were assigned to unsampled rivers (27 rivers in the Southwest, corresponding to an annual discharge of 19 km³ yr⁻¹ or 11% of the regional total; 50 rivers in the East, corresponding to an annual discharge of 89 km³ yr⁻¹ or 61% of the regional total) prior to the calculation of transport for each season. For both methods, annual transports were calculated by the sum of the seasonal totals. To obtain bay-wide transports of nutrients, either the estimates obtained by Method 1 (n = 149) were added to the fluxes for the nine representative rivers, or the estimates obtained by Method 2 (n = 77) were combined with those obtained by Method 1 for northwestern (n = 41) and southern (n = 31) rivers and those of nine representative rivers. Overall, the seasonal and annual fluxes provided by Methods 1 and 2 were comparable (Table S6), providing confidence that the simplifications inherent to the former allow for a valid coarse estimation.

Bay-wide estimates of nutrient transports were obtained by combining the fluxes estimated for the nine representative rivers and those obtained by Method 1 for the others. Given the inconsistency of measurement methods for some N forms (TN versus TDN or DON) and the fragmentary spatial and seasonal coverage for these forms and NH₄⁺, only NO_x was considered for establishing bay-wide transport estimates. However, the possible contribution of NH₄⁺ and DON to the seaward flux of available nitrogen is discussed based on the available data.

The uncertainty in discharge was determined based on one standard deviation calculated from the ensemble-modelled discharge data set. For nutrient transport, the uncertainties for seasonal nutrient fluxes in each river were propagated as follows:

where Q denotes the seasonal mean discharge; C, the seasonal mean nutrient concentration; and σ_Q, and σ_C, their associated errors, respectively. Given the limitations of the data sets, an estimation of error (or uncertainty) was only possible for the rivers in which at least three nutrient measurements existed for a given season (i.e., the Churchill, Nelson, and Hayes rivers). For a given river, the error was taken as one standard deviation from the average of all measurements available in a season. The average of relative errors (⁠$σcC$ = 0.5) for these rivers was then used for the other representative rivers where sampling occurred less than three times for a given season, as well as for all the other rivers in HB. Uncertainties for annual nutrient fluxes and regional or bay-wide totals were calculated by summing the errors attached to all relevant components.

Data were organized using the “tidyverse” package (Wickharm et al., 2019), and the Kruskal-Wallis test was performed for non-parametric comparison in R (R Core Team, 2020). We also used “ggplot2” package in R for data visualization (Wickharm, 2016). For consistency in the presentation of tables and figures, the names of rivers or groups of rivers were ordered and displayed anti-clockwise starting at the northwest end of Hudson Bay (i.e., the Wilson River).

Nutrient concentrations in rivers ranged from <0.03 (or below detection limit) to 12.9 μM for NO_x, from <0.05 to 0.97 μM for PO₄³⁻, and from 0.1 to 95.7 μM for Si(OH)₄ (Tables S1–S3). This variability was caused mainly by the rock formations, climates, and diversity of drainage basins surrounding HB and to a lesser extent by seasonality within rivers (Figure 2). These drivers are all subject to considerable differences across the basins that drain into HB. The HBDB includes the Canadian Shield, Interior Plains, and Hudson Bay Lowland (right-hand panel in Figure S1). Most HB rivers drain the Canadian Shield, which consists of poorly weatherable, intrusive igneous and metamorphic (crystalline) rocks of Precambrian origin, with some imbedded volcanic and sedimentary assemblages (Wheeler et al., 1996). Overlying soils are influenced by the geologic and glacial history and tend to be shallow. In the southwestern part of HBDB, the Interior Platform (the upstream part of the Churchill and Nelson river basins) consists of fine-grained and poorly consolidated sedimentary rocks (e.g., limestones, sandstones, and shales). Well-developed mineral soils (turbic crysols, luvisols, and brunisols) and peatland complexes composed of peat plateaus, flat bogs, and channer fens are characteristic of this region (Baldwin et al., 2020). The Hudson Bay Lowland in the southwestern and southern parts of Hudson Bay (the lowermost sections of the Nelson, Hayes, Severn, and Winisk river basins) consists of flat deposits of Paleozoic limestones and dolomites with an organic-rich soil above which exists a vast wetland/peatland complex.

Mean annual air temperatures range from 4°C in the Canadian Prairies (warm) to −10°C in Nunavik (cold), while mean annual precipitations range from 400 mm in northern Nunavik (dry) to 800 mm in the Boreal Forest region (wet; Brown et al., 2012; Stadnyk et al., 2019). Permafrost throughout the HBDB ranges from continuous in the north to fully absent in the south. The wide latitudinal range in climate and permafrost coverage spanned by the various HB river basins leads to a similarly wide range in vegetation among the various subbasins (left panel in Figure S1). Proceeding from north to south, we find the Arctic shrub tundra zone (Wilson, Ferguson, The-Anne, and Thlewiaza in northwest, Nastapoka, Innuksuac, and Povungituk in northeast, and Foucault and Deception in Hudson Strait), the subarctic tundra-woodland zone (Seal, Knife, and Churchill in the southwest, and Grand rivière de la Baleine and Petite rivière de la Baleine in the east), and the boreal forest and woodland zone (Nelson and Hayes in the southwest, Severn and Winisk in the south, and Eastmain, Vieux Comptoir, Maquatua, Castor, La Grande, Piagochioui, and Roggan in the east).

The impact of environmental drivers on the chemical composition of HB rivers is underscored by the spatial pattern in nutrient concentrations and molar ratios. During the BaySys and COast-JB expedition during spring and early summer (June to early July), for example, concentrations among rivers varied by two orders of magnitude for NO_x (min of <0.03 μM in the east, max of 3.2 μM in Hudson Strait) and Si(OH)₄ (min of 0.1 μM in the east, max of 39.1 μM in the south) and a factor of 6 for PO₄³⁻ (min of <0.05 μM in the northwest, max of 0.28 μM in the southwest). While PO₄³⁻ was generally depleted in rivers, the Nelson, Hayes, and Severn rivers showed >0.1 μM of PO₄³⁻. This finding is consistent with the fact that phosphorus levels in fresh waters are augmented in lowland areas with sedimentary rock deposits and organic-rich soils but are generally lowest in crystalline bedrock areas (Wetzel, 2001). Similarly, the southwestern and southern rivers showed elevated Si(OH)₄ concentrations compared to other rivers (28.5 μM and 35.4 μM versus 4.3 μM, 13.5 μM, and 17.8 μM on average for the northwestern HB, eastern HB, and Hudson Strait rivers). This contrast reflects the different rates of chemical weathering across river basins. Rosa et al. (2012) reported that rock denudation rates correlate positively with the percentage of volcanic and sedimentary rock cover (5.3–5.6 t km⁻² yr⁻¹ for the Interior Platform and Hudson Bay Lowland versus 1–4 t km⁻² yr⁻¹ for the Canadian Shield), which explains the high Si(OH)₄ concentrations in southwestern and southern rivers.

Molar ratios of the different nutrients reflected the spatial patterns discussed above (left panel in Figure 3). The rivers of northwestern HB and Hudson Strait had particularly high N:P (12.2 to 64.9) and low Si:N (0.7 to 29.9) ratios compared to those of the southwest and southern HB (N:P ranging from 1.2 to 6.4 and Si:N ranging from 56.2 to 276.8). Moreover, a spatial contrast in the relationship between NO_x and NH₄⁺ concentrations suggests that changes in vegetation type and abundance, as well as in the rate of biological activity, also affect nutrient concentrations. Northwestern rivers, Hudson Strait rivers, and the Povungnituk River in the northeast exhibited relatively high NO_x concentrations with wide-ranging NH₄⁺ concentrations, whereas the southwestern and southern rivers had reduced NO_x concentrations (<0.5 μM) with a relatively narrow range of NH₄⁺ concentrations (right panel in Figure 3). The northern part of river basins located above the treeline (above approximately 60°N) consists of sparse low-to-ground vegetation, and the terrestrial and aquatic ecosystems there were at a relatively early stage of seasonal development at the time of sampling due to the delay in the onset of suitable growing conditions (left panel in Figure S1; Baldwin et al., 2020). This situation results in NO_x concentrations not being depleted in contrast with the southwestern and southern rivers located in the boreal forest and woodland zones. The latter harbor relatively productive and diverse vegetation (e.g., higher plants, algae) that captures NO_x and NH₄⁺ through assimilation. These spatial patterns in NO_x and NH₄⁺ concentrations were inverse to those observed for DON (Figure 4). While a statistical analysis for the difference in DON concentrations between the grouped rivers was not feasible due to a small number of observations, DON concentrations were high in the southwestern and southern rivers (mean of 19.7 μM) where DON contributed nearly all of the TDN (mean of 98%). By contrast, the mean concentration of DON (7.0 μM) and its contribution to TDN (72.4%) were relatively low in the northwestern rivers, Hudson Strait rivers, and the Povungnituk River.

Nutrient loads in rivers and lakes are affected by chemical weathering, soil leaching, runoff from adjacent basins or deposition from the atmosphere (mainly nitrogen), mineral precipitation and biological activity (Frey and McClelland, 2009; Makarov, 2009). Nutrient concentrations showed seasonal patterns in the few rivers that were sampled on more than one occasion within a year (Figure 2). Concentrations of NO_x were highest in winter (seasonal mean of 5.0 μM; Kruskal-Wallis test, P < 0.05), when nitrogen remineralization by microbial activity exceeds nitrogen demands by autotrophs, decreased during spring through the onset of biological uptake within the rivers, and remained relatively low during summer and fall. Concentrations of Si(OH)₄ were lowest (mean of 17.8 μM) during the spring freshet and increased thereafter. The highest values (mean of 47.6 μM) occurred in the fall (Kruskal-Wallis test, P < 0.05), presumably due to input from active chemical weathering of silicate and aluminosilicate minerals driven by jointly increasing temperature and runoff (White and Buss, 2014), and held constant during winter (mean of 38.3 μM). By contrast, PO₄³⁻ concentrations did not show distinct seasonality and were generally low throughout the year, owing to the fact that P is generally the limiting element for biological productivity in freshwater, where it is usually fully depleted during the productive season (Schindler, 1978; Schindler et al., 2008). These nutrient-specific features are similar to those observed in the Arctic and Eurasian subarctic rivers (see Figure 1 in Le Fouest et al., 2013, and Figure 2 in Humborg et al., 2003) and suggest that physical and biochemical processes underlying nutrient variability are generalized across the circumpolar domain.

The mean molar C:N:P composition of 106:16:1 (or Redfield stoichiometry; Redfield et al., 1963) that characterizes the organic content of marine phytoplankton is used as a reference point to assess whether a particular ambient mixture of inorganic nutrients can satisfy algal demand, on average. Although the optimal stoichiometry of different algal species or groups varies widely (e.g., Geider and Laroche, 2002), ambient nutrients of N:P ratios >16 are generally considered to lead to P limitation and, conversely, those of N:P ratios <16 to lead to N limitation. This concept was expanded to encompass a critical Si:N threshold of 1 (Brzezinski, 1985), below which the diatoms that account for the majority of oceanic net primary production in temperate and polar waters (Sigman and Hain, 2012) may become Si-limited. Seasonal patterns were apparent in molar nutrient ratios (Figure 5), which are first calculated on the basis of NO_x alone to allow for comparison with other studies and the use of a greater number of rivers (i.e., other forms of nitrogen were measured in a subset of rivers). The N:P molar ratios in dissolved inorganic nutrients decreased from maximum values in winter (22.3 on average) to the lowest ones during fall (4.7 on average), whereas Si:N molar ratios showed an opposite trend with average values increasing from 7.7 during winter to 44.5 during fall. The results indicate that rivers are N-deficient (N:P < 16 and Si:N > 1) except during winter, when the transport of nutrients toward the coast then sets the stage for a sizable spring bloom in which the use of riverine PO₄³⁻ and Si(OH)₄ is maximized. As biological activity is reduced during winter, these nutrient inputs may disperse farther into HB than they otherwise do during the productive season. During the latter season (from spring to fall), however, the low N:P and high Si:N molar ratios constrain the use of PO₄³⁻ and Si(OH)₄ and allow for a differential accumulation of these nutrients in the bay.

NH₄⁺ and labile DON are known to be important nitrogen sources for primary production. We reassessed nutrient molar ratios including these nitrogen forms (considering that the sum of NO_x, NH₄⁺, and labile DON represents bioavailable nitrogen). Because previous studies of Arctic rivers have reported low and relatively invariant concentrations of NH₄⁺ throughout the year (Holmes et al., 2012; Wickland et al., 2012), with the exception of the Yukon and Ob’ rivers during the winter and spring freshet (ArcticGRO Water Quality Dataset: Holmes et al., 2022), our average springtime concentration of NH₄⁺ (0.6 µM, the Hudson Strait rivers excluded) can presumably be used for all seasons (Table S1). Given that DON concentrations are known to peak during the spring freshet and to generally decrease thereafter (Le Fouest et al., 2013), our springtime measurements of DON concentration (15 μM on average) should be considered as maximum values and possibly overestimates for other seasons. While the calculation of DON concentration is sensitive to the combined analytical error of three separate analyses (TDN, NO_x, and NH₄⁺), the generally low NO_x and NH₄⁺ concentrations relative to TDN measured during spring results in a low aggregate uncertainty (Lee and Westerhoff, 2005; Khosh et al., 2017).

To estimate the DON fraction that can be considered as labile within hours to a few weeks, we considered the relationship between salinity and DON concentration in the Nelson Estuary during the spring freshet. This approach does not formally assess the bioavailability of DON but assumes that local departures from the conservative mixing line between freshwater and marine end-members indicate net gains (points above the line) or losses (points below the line). The nearly conservative behavior of DON in this analysis showed that biological processes had a minor effect on concentrations during the sampling period. Negative departures ranging from 1 to 10% (mean and standard deviation of 6% ± 4%, n = 5) of the expected conservative value occurred throughout the salinity range, but a net release corresponding with a peak in NH₄⁺ was also apparent at a salinity of about 15 (Figure S2). The conservative estimate of 6% labile DON thus obtained is at the low range of values reported for other rivers across the Arctic using a variety of methods (4% to 38%; Tremblay et al., 2014; Kaiser et al., 2017; Sipler et al., 2017). Our estimate also contrasts with the high lability of dissolved organic carbon observed using the incubation method during late winter (mean of 30% lability in the Nelson River and estuary including the Hayes River; Kazmiruk et al., 2021). We speculate that DON was less available microbially due to a shift in dissolved organic matter composition at the time of sampling (late June). Another possible explanation is that DON inputs frequently occur through desorption from sediment in the estuary, and thus DON concentrations in the water column are not reduced. The latter is supported by the minor DON increase concomitant with the NH₄⁺ peak at intermediate salinities (Rysgaard et al., 1999) and implies some degree of uncoupling between net change in concentration of the reduced forms of N (i.e., NH₄⁺ and DON) and N cycling rates (e.g., uptake, regeneration, and ammonification) in the Nelson estuary (Lee et al., in preparation).

Using our estimates of NH₄⁺ and labile DON above, the revised N:P molar ratios (the sum of seasonal mean NO_x, NH₄⁺ (0.6 μM), and labile DON (0.9 μM) values divided by the seasonal mean PO₄³⁻ value using the binned data from Tables S1–S3) are 29.0, 25.3, 16.8, and 11.2 for winter, spring, summer, and fall, respectively. In contrast with the analysis based on NO_x only, these revised values indicate that river waters exhibit P shortage in winter and spring, have nearly balanced N:P molar ratios in summer, and experience N shortage in the fall. Revised Si:N molar ratios (5.9, 7.3, 7.4, and 18.5 for winter, spring, summer, and fall, respectively) remain well above the potential limitation threshold for diatoms.

The multi-annual data from the Churchill, Nelson, and Hayes rivers collected during the BaySys, CAMP, and Conawapa GS projects can be used to probe into possible temporal trends in nutrient concentrations. Because the sampling locations and dates of sampling varied across these projects, we focused on the winter period to minimize the possible masking influence of high environmental and biological variability from spring to fall (e.g., precipitation, atmospheric deposition, evaporation, weathering, snowmelt events, tributary mixing, and biological consumption). Nutrient concentrations did not differ significantly between sampling years in all rivers (Kruskal-Wallis test, P >0.05 for NO_x, PO₄³⁻, and Si(OH)₄; Figure 6) and no temporal trends were apparent over the periods 2009–2017, 2004–2017, and 2010–2017 in the Churchill, Nelson, and Hayes rivers, respectively. As these comparisons suggest the absence of climate-driven changes or cyclicity in the concentrations observed during winter, the fluctuations observed are therefore attributed to interannual variability. The available time series is relatively short, however, and would need to be extended into the future and better resolved seasonally to draw firm conclusions.

We did not find any clear evidence that nutrient concentrations differ between regulated (Nelson and La Grande) or partially diverted (Churchill and Grande rivière de la Baleine) rivers and the others (Kruskal-Wallis test, P > 0.05). While the Nelson River showed the highest NO_x and PO₄³⁻ concentrations of all rivers on some sampling occasions (Figure 2), the seasonal means were not significantly higher than those of nearby rivers (Kruskal-Wallis test, P > 0.05). The variability in NO_x and PO₄³⁻ for the Nelson River is linked to processes that occur far upstream in Lake Winnipeg and the Nelson watershed, as the CAMP and Conawapa GS data showed nearly spatial homogeneity in nutrient concentrations along the length of the river (not shown). There are two plausible explanations for the peaks in NO_x and PO₄³⁻ concentrations. Firstly, Lake Winnipeg is located in a lowland area among sedimentary rock deposits and boreal shield formations that hold numerous bogs and wetlands known for their abundance of organic phosphorus (Wetzel, 2001), showing high TP concentration (McCullough et al., 2012). Periodical release and mineralization of this phosphorus possibly modulates PO₄³⁻ concentration in the river system. Secondly, the Nelson watershed includes the Temperate Prairies ecozone where cultivated crops dominate land use (Déry et al., 2016). The leaching resulting from fertilizer inputs possibly contributes to the peaks in NO_x and PO₄³⁻ concentrations, which always co-occurred during the crop-growing season (Figure 2). This idea is consistent with the fact that seasonality in Si(OH)₄ concentrations was similar in all rivers of the southwest, including the Nelson River, whereas seasonally averaged Si:N and Si:P ratios for dissolved inorganic nutrients were generally lower in the Nelson (Si:N < 16 and Si:P < 80) than in the Churchill or Hayes rivers (Si:N > 25 and Si:P > 100) during the productive period (Figure 7). These patterns in nutrient concentrations and ratios were not observed in the regulated river of eastern HB, where no agriculture occurs in the watershed.

Given the diversity of processes influencing the Nelson River and the absence of pre- versus post-damming assessments, to discern a possible impact of regulation on in-river nutrient concentrations is not possible. The results also showed that Si:N ratios in winter were generally low and similar in all rivers, implying that differences in flow regime do not noticeably impact silicate loading and nitrogen cycling. While Si:N and Si:P ratios were varied across rivers during the productive period (spring to fall), the lack of data on the separation of regulation effects from biological uptake (high autotrophic demand for inorganic nutrients) precluded an analysis of contribution of different flow regimes to the variation in Si:N and Si:P ratios.

To provide context for the interpretation of riverine inputs in HB, we compared our observed nutrient data for the Churchill, Nelson, and Hayes rivers (where the seasonal resolution of sampling was relatively high) with those of other rivers at the pan-Arctic scale. Estimated mean annual concentrations of NO_x for the North American Arctic, Eurasian Arctic, and Eurasian subarctic rivers are 7.5 μM, 6.1 μM, and 4.7 μM, respectively (Table 1), which exceeds the HB average (2.9 μM) by two-fold approximately. Similarly, annual mean concentrations of Si(OH)₄ were nearly 2-fold higher in the Eurasian subarctic rivers (74.3 μM) and 6-fold higher in the North American Arctic and Eurasian Arctic rivers (190.3 μM and 211.6 μM), respectively, compared to HB (35.1 μM). By contrast, the annual average of PO₄³⁻ in HB (0.31 μM) was on par with the value for Eurasian Arctic rivers (0.31 μM) and twice as high as the values for the North American Arctic and Eurasian subarctic rivers (0.12 μM and 0.14 μM).

Distinct nutrient-specific patterns across the different watersheds of the pan-Arctic region can be explained by environmental and biological drivers similar to those discussed above. While average concentrations of NO_x during winter (defined as November to April here) and summer–fall (defined as July to October here) were respectively high and low across all pan-Arctic regions, rivers of the Eurasian and North American Arctic, and to a lesser extent the Eurasian subarctic, exhibited high peak values (8–25 μM) compared to the HB rivers (<7 μM). Winter plays a pivotal role for the seaward transfer of nutrients due to long periods of plant dormancy in the terrestrial system, a low nutrient demand by autotrophic organisms in rivers, and continued nutrient recycling by microbes at low temperatures (e.g., high ammonification and nitrification rates in terrestrial and freshwater systems; Cavaliere and Baulch, 2019; Sanders et al., 2021; Xu et al., 2021). Rivers of the Eurasian subarctic and the Eurasian and North American Arctic showed elevated mean concentrations of winter NH₄⁺ (ranging from 0.7 μM to 19.9 μM versus ≤0.5 μM for the HB rivers based on a single sampling event), suggesting relatively high ammonification rates in the rivers or soils there. Nitrification of this NH₄⁺ enriches the NO_x pool and may explain why the latter is relatively high there. These distinct winter patterns in NO_x and NH₄⁺ concentrations suggest that contrasted rates of net remineralization resulting from different distributions, community structures, and activities of microorganisms between watersheds account for inter-regional differences in river NO_x concentrations. While groundwater can be an important source of winter dissolved inorganic nitrogen (DIN) in some systems (Dornblaser and Striegl, 2007), the lack of data on this process in HB precludes an analysis of its contribution to regional differences.

Differences in annual averages of Si(OH)₄ concentration between rivers can be explained by chemical weathering. Eurasian and North American Arctic river basins typically have elevated chemical weathering rates compared to the Hudson Bay river basins, including the Canadian Shield and Interior Platform (Millot et al., 2003; West et al., 2005; Rosa et al., 2012; Kang et al., 2022). While Eurasian subarctic river basins consist predominantly of gneiss, granite, and acid volcanic rock that exhibit low weathering rates compared to volcanic and shield rocks (Gaillard et et al., 2003; Humborg et al., 2004), some rivers in this region nevertheless exhibit relatively high annual Si(OH)₄ concentrations. These high values are due in part to a greater proportion of coniferous forests (>30% of land cover) over broad-leaved deciduous forests (<15% of land cover), which enhances the weathering rate of resistant rocks (Humborg et al., 2004). This enhanced rate indicates that the type of vegetation within watersheds can be another controlling factor for the variation in Si(OH)₄ concentrations. Moreover, the climate is one of the important drivers to affect variation in Si(OH)₄ concentrations. By virtue of their position at the southern margin of the pan-Arctic region, HB watersheds experience a relatively warm climate and a protracted productive season that favor biological activity and Si(OH)₄ consumption in rivers, leading to low extremes in the annual average concentration for this nutrient.

For PO₄³⁻, contrasts across the high-Arctic region are likely caused by differences in the extent of permafrost-influenced or -free peatlands between the major Arctic river drainage basins and the Ob’ River drainage basin. Permafrost peatlands are distributed over large parts of the northern circumpolar region, holding greater amounts of organic soil mainly in the form of frozen peat (Joosten, 2018), which likely causes a decrease in hydraulic conductivity of upper peats that then reduces phosphorus transport. By contrast, upstream portions of the Ob’ River basin (below 60°N) lie within permafrost-free peatland and a warm climate zone, which have favorable conditions for hydraulic conductivity of nonfrozen peats as well as net mineralization by enhancing microbial activity. These localized processes in nonfrozen peatland areas provide a relatively large source of phosphorus (Frey et al., 2007) and drive the elevated annual averages of PO₄³⁻ in the Ob’ River. As a result, PO₄³⁻ concentrations tend to be high in rivers of the Eurasian Arctic. While river basins of the Eurasian subarctic and HB are composed of a large proportion of permafrost-free peatlands, the rivers in the former show low annual averages of PO₄³⁻ concentrations. Previous modeling studies found that anthropogenic disturbances play a vital role in the change of nutrient delivery along the water pathway (Sferratore et al., 2008), suggesting that trapping phosphorus into dams and reservoirs leads to variation in PO₄³⁻ concentrations. Furthermore, shifts in vegetation types and losses of vegetated soils resulting from inundation when reservoirs and dams are built are known to impact dissolved constituents in this region (Humborg et al., 2002; Smedberg et al., 2009). These impacts possibly explain the relatively low annual average of PO₄³⁻ concentration there. As discussed above, the results reported for Eurasian subarctic rivers contrast with our findings for the Nelson River, implying that localized upstream features such as large lakes (e.g., Lake Winnipeg) can substantially affect nutrient loading. However, the monitoring a possible shift in PO₄³⁻ concentration by cumulative alterations (e.g., flow regulation and damming) in relation to the local change in vegetation type in the Nelson watershed is needed.

Assessments of the influence of discharge on nutrient transport within a year, inter-annually, or in the context of climate change can take advantage of the correlation that is sometimes observed between discharge and nutrient concentrations. This relationship implies that nutrient transport does not simply vary in direct proportion to discharge (multiplied by a mean concentration) but is modified by the impact of streamflow on the balance of processes that add or remove nutrients in a river. For example, increased river flow can affect the concentrations of dissolved inorganic nutrients either negatively (dilution effect in high-latitude rivers; e.g., Sferratore et al., 2008; Holmes et al., 2012) or positively (flushing effect in agriculture-dominated watersheds; e.g., Outram et al., 2014). Here these effects can only be explored using nutrient data collected from the CAMP and Conawapa GS projects, for which comparisons between observed discharge at hydrometric gauge stations (Table S2, the water survey of Canada database; https://wateroffice.ec.gc.ca/) and discrete nutrient measurements can be established for the Churchill, Nelson, and Hayes rivers. Contrary to previous efforts of this type at mid-latitudes (e.g., Sigleo and Frick, 2007), we were unable to find any significant relationships using the data at hand. Other studies based on a high frequency of sampling also failed to detect relationships at high latitudes (e.g., the Mackenzie River; Tremblay et al., 2014) or in agriculture-dominated watersheds (Scholefield et al., 2005; Johnes, 2007; Wade et al., 2012). The reasons are unclear but we speculate that flow regulation or agricultural nutrient inputs in some of the watersheds mask the aforementioned effects of discharge in unperturbed systems.

Given the absence of discharge-concentration relationships, we elected to use a seasonal binning approach for matching nutrient concentration with discharge data, and we assigned the average concentrations of all HB rivers to those with no data. Despite its necessary simplicity, this approach reasonably captured seasonal patterns in both nutrient concentrations and discharge based on observations and thus allowed to assess the coarse temporal distribution of nutrient transports throughout the year. Because observational discharge records at gauged hydrological stations are not available downstream in most HB rivers (Figure 1), the bay-wide upscaling exercise here relied on the simulation of discharge using the HYPE model driven by three different outputs from the fifth Coupled Model Intercomparison Project (CMIP-5). The modeling approach allows for simulation of discharge dynamics at the bay-wide scale, showing improvements in comparison with observed discharge as well as comparable performance based on reanalysis of discharge datasets (Ridenour et al., 2019; Stadnyk et al., 2020; Tefs et al., 2021). Despite these substantial efforts to calibrate the models, the discharge estimates obtained still contain some uncertainty due to the scarcity of observational streamflow data available for ground-truthing. Thus, our approach leads to some underestimation or overestimation of the flux at different times of year (Figure S3 and Tables S7 and S8). It was, however, the most robust way to combine spotty nutrient concentration data with high-resolution model results (daily estimates for all rivers) of discharge at the bay scale.

In order to compare nutrient dynamics among rivers, we focused on nine rivers for which seasonal coverage was sufficient in nutrient data (interpolated values included). The selected rivers also represented significant contrasts with respect to watershed, size (annual discharge), the presence or absence of human disturbance, and the type of disturbance with respect to fragmentation and the regulation of streamflow. Based on the HYPE modeling results, the coefficient of variation (CV), which is the ratio of the standard deviation (σ) to the mean (μ), was calculated to assess annual and seasonal variability in discharge among models (n = 3). The seasonal variability (CV) in discharge was generally lower than 35% with exception of the Churchill (CV > 40% from GFDL-CM3 in all seasons and from MIROC5 in summer and fall), Maquatua (CV > 40% from GFDL-CM3 in all seasons and from MRI-CGCM3 in fall), and Castor (CV > 50% from GFDL-CM3 in all seasons and from MIROC5 as well as MRI-CGCM3 in fall) rivers (Table S9). The annual variability in discharge was generally less than 25% with exception of the Churchill (CV > 30% from GFDL-CM3 and MIROC5), Maquatua (54% from GFDL-CM3), and Castor (74% from GFDL-CM3) rivers. Overall inter-seasonal variability in discharge within rivers was higher than the variability of annual discharge across rivers in all model results, indicating the highly dynamic character of the HBDB hydrologic regime within a year.

Annual discharge and nutrient transports among the nine rivers varied by three and four orders of magnitude, respectively (Table 2). The Nelson River had the highest discharge (111 km³ yr⁻¹) and nutrient fluxes (7008 × 10⁶ g N yr⁻¹, 1961 × 10⁶ g P yr⁻¹, and 1371 ×10⁸ g Si yr⁻¹). Despite the fact that the La Grande River showed the second highest annual discharge (50 km³ yr⁻¹) and NO_x flux (2422 × 10⁶ g N yr⁻¹), it was outmatched by the Hayes and Severn rivers for annual PO₄³⁻ and Si(OH)₄ transports (175 × 10⁶ g P yr⁻¹ and 329 × 10⁸ g Si yr⁻¹). This result underscores that annual differences in discharge among rivers do not simply translate into proportional differences in annual nutrient flux, which is more variable due to seasonal patterns and watershed-based difference in nutrient concentrations. The seasonality in discharge also differed between rivers (Figure 8). Nearly half of the total annual discharge in the eastern rivers took place during spring (40% for the partially-diverted Grande rivière de la Baleine and 48% and 47% for Maquatua and Castor rivers, respectively), when sharp peaks occurred in the hydrograph, whereas in the southwestern and southern rivers (the partially-diverted Churchill River and the Hayes, Severn, and Winisk rivers) only a quarter of the annual discharge occurred during spring owing to broader peaks in the hydrograph starting in late winter and extending into early summer.

In regulated rivers, discharge was clearly skewed toward the winter, during which 50% and 65% of the annual water flow occurred in the Nelson and La Grande rivers, respectively. This shift, which results from peak demand in hydropower production during cold months, was most pronounced in the La Grande River system owing to the large storage capacity of its reservoirs versus the “run-of-river” regulation employed in the Nelson River (Déry et al., 2016). This seasonal shift in discharge had a significant impact on the proportion of annual NO_x, PO₄³⁻, and Si(OH)₄ fluxes occurring during winter, which respectively amounted to 76%, 56%, and 57% for the Nelson River and 82%, 71%, and 68% for the La Grande River. Seasonal fluctuations in PO₄³⁻ and Si(OH)₄ fluxes were generally proportional to those of discharge, with some minor or major deviations in spring and summer caused by biological uptake (lowers concentrations) or by catchment runoff and soil leaching processes (increases concentrations). By contrast, the seasonality of NO_x transports was generally decoupled from the seasonality of discharge, presumably owing to a relatively large impact of biological processes on this nutrient during the vegetative period (spring to fall). The decoupling in NO_x transport was most pronounced for the Churchill, Nelson, and Hayes rivers, where only 25% of the annual flux occurred during the vegetative period. Such decoupling was generally less pronounced in other rivers, although reduced NO_x transports were apparent during spring in the La Grande River and summer in the Maquatua River.

Our analysis of the regional partitioning of river discharge into HB showed that northwestern rivers accounted for less than 8% of the bay-wide annual total (Figure 9), with a similar contribution for the transport of each nutrient (range of 1% to 2%, discharge-normalized nutrient transports). Contributions of the southwestern, southern, and eastern rivers to overall discharge were much higher than in the northwest and comparable across the three sectors (range of 25% to 40%). Southwestern rivers made a disproportionate contribution to the transport of NO_x (46% to 58%) and PO₄³⁻ (41% to 74%), except during spring when their NO_x transport (29%) was similar to that of the eastern rivers (24%) and their PO₄³⁻ transport (41%) was similar to that of the southern rivers (38%). Southwestern and southern rivers (range 26% to 53%) made relatively high and similar contributions to total Si(OH)₄ transport, except during winter and spring when their contributions to total Si(OH)₄ showed a noticeable difference (50% and 27% in southwestern rivers versus 26% and 53% in southern ones), with eastern rivers accounting for 20%. The eastern sector was a minor source in the winter and spring when its contributions were only half those of the other sectors. Overall, on an annual basis, southwestern rivers collectively made the largest single contribution to nutrient deliveries, followed by modest contributions from the southern and eastern rivers, and minor contributions from the northwestern ones (Table 3). The winter season (defined here as a 6-month period from November through April), which was 3 times longer than each of the other 2-month periods (i.e., spring, summer, and fall), made the single largest contribution to annual bay-wide discharge (41%) and nutrient fluxes (74%, 50%, and 47% for NO_x, PO₄³⁻, and Si(OH)₄, respectively).

On average, nutrients in the ocean are consumed by phytoplankton according to a C:N:P:Si stoichiometry of 106:16:1:16, based on the work of Redfield et al. (1963) and Brzezinski (1985), recognizing that departures may occur locally or episodically depending on assemblage composition and growth conditions. Primary production is considered to be limited by the first nutrient to be fully depleted, following Liebig’s Law of the Minimum. The molar ratios of the different nutrients transported by rivers can therefore be used as a benchmark to assess which nutrient is likely to be used up first or to remain in excess in coastal areas. In this regard, the N:P molar ratio of riverine-transported nutrients varied between rivers and seasons (Figure 10). During winter, this N:P molar ratio was higher than the Redfield ratio of 16:1, except for southwestern rivers (12:1). From spring to fall, the N:P molar ratios were systematically below the Redfield ratio, indicating insufficient NO_x to support the full consumption of PO₄³⁻ in the river waters delivered to the coast. The Si:N molar ratio of nutrient transports systematically exceeded the average marine diatom consumption ratio (1:1) by a large margin, indicating that rivers provided much more Si(OH)₄ than diatoms would be able to consume given the amount of NO_x available.

So far, our analysis has not included the riverine NH₄⁺ and labile DON fluxes, which can contribute to primary production and lower the residual concentrations of riverine PO₄³⁻ and Si(OH)₄ in coastal waters. In the absence of sufficient data to calculate bay-wide, annual inputs of NH₄⁺ and labile DON to HB, their riverine fluxes can be approximated by multiplying overall NO_x transports by the NH₄⁺:NO_x and DON:NO_x ratios from the subset of samples that were collected from rivers during spring (see discussion above). The resulting estimates of NH₄⁺ transport thus obtained are 1.2 × 10⁸ mol N, 0.7 × 10⁸ mol N, 0.6 × 10⁸ mol N, and 0.7 × 10⁸ mol N for winter, spring, summer, and fall, respectively (by multiplying seasonal NO_x transport by 0.6 μM and dividing it by the seasonal mean NO_x concentration and the molar weight of nitrogen; Tables 3 and S5). The fluxes of labile DON would be on the order of 1.9 × 10⁸ mol N, 1.1 × 10⁸ mol N, 0.9 × 10⁸ mol N, and 1.0 × 10⁸ mol N for winter, spring, summer, and fall, respectively. The annual riverine input of bioavailable N (the sum of seasonal NO_x, NH₄⁺, and labile DON fluxes) was 2.21 × 10⁹ mol N yr⁻¹, which is 1.57 time higher than the estimate based on NO_x alone (1.41 × 10⁹ mol N yr⁻¹; Table 3). The latter is 46% lower than the previous annual estimate (2.62 × 10⁹ mol NO_x yr⁻¹) made by Kuzyk et al. (2010), a difference that we attribute to the greater seasonal and spatial resolution of the data available for the present study as well as the lower estimate of discharge used (e.g., the total annual discharge of 503 km³ yr⁻¹ for this study versus 713 km³ yr⁻¹ for Kuzyk et al., 2010).

Considering the indirect estimates of NH₄⁺ and DON above, the combined molar transports of bioavailable N, PO₄³⁻, and Si(OH)₄ by all rivers during winter and spring can provide an estimate of what is available for the spring bloom in the freshwater domain of estuaries (under tidal influence but with a salinity of 0). The combined molar transports amount to 16.4 × 10⁸ mol N, 0.8 × 10⁸ mol P, and 104.1 × 10⁸ mol Si, respectively. Assuming a C:N:P:Si stoichiometry of 106:16:1:16, full consumption of the winter/spring PO₄³⁻ supply by phytoplankton would leave behind residual pools of bioavailable N (3.6 × 10⁸ mol N) and Si(OH)₄ (91.3 × 10⁸ mol Si) that correspond to 22% and 87.7% of their riverine supplies, respectively. The situation is different for the summer and fall, when overall molar fluxes of dissolved nutrients are poor in bioavailable N (5.7 × 10⁸ mol N) relative to PO₄³⁻ and Si(OH)₄ (0.4 × 10⁸ mol P and 59.3 × 10⁸ mol Si, respectively). Following the bioavailable N depletion in this case, a PO₄³⁻ residual (0.04 × 10⁸ mol P or 11% of the river flux) would be left and Si(OH)₄ excess (53.6 × 10⁸ mol Si or 90% of the river flux) would remain high. This analysis shows that in the absence of other nutrient sources, the freshwater domain of estuaries would shift from a P-limited state during winter and spring to N-limited condition toward late summer and fall. In the salt transition zone, however, the admixture of marine waters, which are also unbalanced in terms of N:P ratio for dissolved inorganic nutrients, would allow the phytoplankton to tap into the residual riverine N during spring. In HB, Foxe Basin and Hudson Strait are the only possible sources of marine waters. These peripheral regions systematically present an excess of PO₄³⁻ at the surface and a systematically low N:P ratio of dissolved inorganic nutrients (<14) throughout the year (see Figure 8 in Tremblay et al., 2019, and Table 1 in Matthes et al., 2021). In this context, the surplus oceanic PO₄³⁻ would allow phytoplankton to consume the unused bioavailable riverine N when marine waters mix with freshwater in estuaries (either by horizontal mixing or vertical entrainment), potentially enhancing primary production during spring. This analysis changes our prior view of nutrient dynamics in HB (i.e., that rivers deliver an excess of PO₄³⁻ with respect to phytoplankton demand), which was based on data from late summer and fall only.

Given that nitrate is the limiting nutrient in marine HB waters (Ferland et al., 2011), the assumption that all of the bioavailable form of riverine N discharged into the bay is eventually used during the productive season is reasonable. The overall annual transport of 3.09 × 10¹⁰ g N or 2.21 × 10⁹ mol N can therefore be used to evaluate the amount of new production that rivers can support, which is equivalent to 1.8 × 10¹¹ g C yr⁻¹ (after multiplication of the molar N supply by the canonical Redfield ratio of 6.6 and the molar weight of carbon). This value can be compared with previously published estimates of bay-wide primary production. The recent work of Matthes et al. (2021) gives an annual value of 72 g C m⁻² based on a compilation of new and historical ship-based measurement of total primary production, whereas the review of Tremblay et al. (2019) gives a range of 10 to 80 g C m⁻² based on modeling and remote sensing approaches. Taking the most recent estimate of Matthes et al. (2021) and excluding the nearshore zone (i.e., considering only the portion of the bay that is deeper than 50 m, 562 × 10³ km²), annual bay-wide primary production would amount to 4.1 × 10¹³ g C yr⁻¹ in marine waters. To assess how much new production this represents, we used the relationship of Eppley and Peterson (1979), which gives an f-ratio (ratio of new versus total primary production) of about 0.18 for an annual level of total primary production of 72 g C m⁻². The corresponding bay-wide new production is on the order of 7.3 × 10¹² g C yr⁻¹, roughly 41 times larger than the new production that can be supported by river N inputs.

While rivers support a small fraction (2.4%) of bay-wide new production, their local impact on the coastal domain can be considerable. Given that N removal typically occurs within a few tens of kilometers of river or delta outlets despite reduced water transparency (Emmerton et al., 2008; Tremblay et al., 2014), annual new production (considering the riverine NO_x flux only in this case) in the Nelson estuary (estimated area of 1000 km² when extending 50 km offshore from the river mouth) potentially amounts to 39.8 g C m⁻² yr⁻¹ or 45.1 g C m⁻² yr⁻¹ when including the contribution of the nearby Hayes River.

This study represents the most comprehensive assessment of riverine NO_x, PO₄³⁻, and Si(OH)₄ fluxes into the HB and their biogeochemical implications for primary production at local and bay-wide scales. We report considerable spatial and seasonal variability in nutrient concentrations at the seaward end of rivers, which are attributed to environmental, geological, and biological drivers as well as to land use in the watersheds. This variability and the lack of pre-development benchmarks in the rivers harnessed for the production of electricity did not allow for the detection of possible impacts of flow regulation on nutrient concentrations. Our data should thus be considered as a contemporary reference point by which future changes relating to flow regulation and/or climate change can be monitored in major rivers (e.g., the Nelson and La Grande rivers) and sectors (e.g., the southwestern and eastern watersheds). One clear effect of flow regulation is to shift water transport and nutrient fluxes toward the winter period due to the relatively high demand for electricity at that time. This shift presumably allows river nutrients to spread over a larger portion of HB than otherwise occurs during other seasons when biological activity is relatively high in estuaries, thereby pre-conditioning the spring bloom over a wider coastal area than in “natural” rivers with relatively low annual discharge.

While river nutrients can make a large contribution to annual new primary production in shallow HB estuaries, the contemporary N supply by rivers can only support a small fraction (ca. 2.4%) of bay-wide new primary production. Despite the large number of rivers and the oligotrophic character of marine waters in HB (Tremblay et al., 2019), its very large size and the fact that NO_x concentrations are particularly low in the rivers of the region limit their overall impact on biological productivity. Based on these observations, we hypothesize that the largest effect of rivers on productivity in the HB system occurs through the modulation of vertical stratification by freshwater, which, contrary to nutrients at the surface, can accumulate far into the marine interior of the bay.

In terms of N:P:Si stoichiometry and the needs of primary producers, river waters contained a large excess of Si(OH)₄ with respect to marine diatom requirements throughout the year. Assessing the balance between N and P deliveries depends on whether one considers NO_x only or the sum of DIN and labile DON. The balance generally shifted from N-limiting to P-limiting conditions when NH₄⁺ and a small fraction of labile DON were included in the calculations. While the lability of DON was not directly assessed here, our indirect estimate derived from the nearly perfect linear relationship between this organic nutrient pool and salinity in the Nelson estuary can be considered conservative, as it falls at the low end of values reported for other systems. This comparison lends credence to the P-limiting scenario and suggests that full use of the bioavailable N provided by rivers must be supported by an admixture of excess PO₄³⁻ from nitrate-depleted Pacific-derived marine waters. Finally, our analysis implies that future studies of river nutrient fluxes into HB should attempt to better characterize organic nutrient pools and seek to monitor major regulated, fragmented, and unimpacted rivers at an intra-annual frequency sufficient to deconvolute the possible effects of regulation and climate change.

Nutrient datasets used in this research are included in this article (supplemental files). All simulation output used in this research is available on SMHI HYPEweb (https://hypeweb.smhi.se/explore-water/historical-data/arctic-long-term-means/) and made available through the BaySys data repository (http://canwin-datahub.ad.umanitoba.ca/).

The supplemental files for this article can be found as follows:

Tables S1–S9. Figures S1–S3. PDF

The authors wish to thank the officers and crew of the Canadian Research Icebreaker, CCGS Amundsen for their excellence in support of this study. They gratefully acknowledge Gabrièle Deslongchamps and Jonathan Gagnon for assistance with sampling and nutrient data processing. They would also like to thank the late Dr. David Barber, who was the driving force behind the BaySys project.

This work was conducted under the auspices of BaySys, a Collaborative Research and Development (CRD) network of the Natural Sciences and Engineering Council of Canada (NSERC) and Manitoba Hydro. BaySys data were completed with datasets obtained from the COast-JB project. Funding for this research was graciously provided by Manitoba Hydro, NSERC, Amundsen Science, Niskamoon Corporation, Québec-Océan (a strategic research cluster funded by Fonds de recherche du Québec), and the Canada Research Chairs program. This work contributes to the scientific programs of Québec-Océan, Institut Nordique du Québec, the ArcticNet Networks of Centers of Excellence, and the Arctic Science Partnership (ASP).

The authors have no competing interest to declare. JET is an Associate Editor at Elementa. He was not involved in the review process of this article.

Contributed to conception and design: JL, JET.

Contributed to acquisition of data: JL, AT, TS, VG.

Contributed to analysis and interpretation of data: JL, AT, TS, VG, MG, JET.

Drafted the article: JL, JET.

Revised the article: All authors.

Approved the submitted version for publication: All authors.

How to cite this article: Lee, J, Tefs, A, Galindo, V, Stadnyk, T, Gosselin, M, Tremblay, J-É. 2023. Nutrient inputs from subarctic rivers into Hudson Bay. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2021.00085

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Christine Michel, Department of Fisheries and Oceans, University of Manitoba, Winnipeg, Canada

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: The Hudson Bay System Study (BaySys)