Plume dispersion from the Nelson and Hayes rivers into Hudson Bay using satellite remote sensing of CDOM and suspended sediment

River water flowing into coastal waters often forms a “river plume,” a buoyant surface layer that overlies the more saline marine water. As the river water in the plume propagates offshore, it gradually becomes more saline by mixing with the ambient seawater (Hetland, 2005). This mixing creates a transitional salinity zone referred to as a region of freshwater influence (ROFI) (Horner-Devine et al., 2015). Several factors, which include river discharge, tidal amplitude and velocity, coastal bathymetry, coastline geometry, wind stress, ocean circulation, water mass stratification, presence of sea ice, and Coriolis force, determine the properties of the ROFI (e.g., Kuzyk et al., 2008; Pimenta et al., 2011; Osadchiev and Sedakov, 2019). Typically, the buoyant layer of river water arises from the density difference between fresh and salty ocean water (Garvine, 1995). The offshore propagation of this buoyant layer results primarily from the tidal and wind-forced drift of the surface plume and mixing with the underlying ocean water (Garvine, 1999). Ultimately, the balance between stratification and mixing drives river plume spread and dispersion in the open ocean. Due to higher river water inflow per unit area and ice meltwater, freshwater content exerts more control over upper water column stratification in the Arctic Ocean than in other oceans (Aagaard and Carmack, 1989; Carmack, 2007; Ehn et al., 2019).

Various tracers can be used to characterize the extent and shape (directionality) of a river plume. As a critical factor in density gradients, salinity is the most widely used physical tracer for river plume and ROFI characterization. However, plume definitions vary as to the salinity difference from offshore ocean water that delineates their extents, e.g., Mississippi River plume (salinity: 32–36), Congo River plume (<30–36), Ganges/Brahmaputra River plume (<26–34), Amazon River plume (30–36), Mekong River plume (32–33) and the Arctic Lena River plume (16–28) (Ahmed et al., 2020; Fournier and Lee, 2021; Zhuk and Kubryakov, 2021). Temperature may be used as a physical tracer of river plumes, but it also varies by the geographical and climatological setting of the rivers (e.g., Otero et al., 2009; Palacios et al., 2012; Hopkins et al., 2013). Other dissolved minor and trace constituents, such as barium (Ba; e.g., Guay and Falkner, 1997; Guay et al., 2009), phosphate (PO₄; e.g., Ikeda et al., 2018), nitrate (NO₃; e.g., Hsiao et al., 2014), aluminum (Al), and silica (Si; e.g., Middag et al., 2009; Bergeron and Tremblay, 2014) are also used to trace sources of riverine freshwater depending on their source concentration.

The rivers debouching into the Arctic Ocean contribute significantly to higher total alkalinity of the surface polar mixed layer (Anderson et al., 2004). Therefore, alkalinity serves as a fluvial water tracer (e.g., Lansard et al., 2012). The non-conservative nature of these chemical tracers, except for dissolved organic carbon (DOC; e.g., Guégen et al., 2011; Milleur et al., 2023), renders them less reliable compared to the more conservative proxies of salinity and temperature. Poleward depletion of the oxygen-18 (¹⁸O) isotope relative to oxygen-16 (¹⁶O) in precipitation has enabled the use of this isotope ratio ¹⁸O/¹⁶O (δ¹⁸O), a chemical tracer, and salinity to discern polar ROFI, which embeds river runoff (e.g., van Donk and Mathieu, 1969; Ostlund and Hut, 1984; Granskog et al., 2009; Granskog et al., 2011). Logistical inaccessibility of Arctic/sub-Arctic coastal waters, and expensive and time-consuming methodologies, further limit the widespread use of these chemical tracers in addressing questions related to spatial-temporal variability of river plumes.

The synoptic coverage of modern space-based optical sensors has proven to be an efficient way to characterize river runoff using optical proxies (e.g., Schroeder et al., 2012). Optical tracking of river runoff in Arctic waters effectively uses the concentration of colored dissolved organic matter (CDOM) due to its conservative nature (e.g., Matsuoka et al., 2012; Fichot et al., 2013). Similarly, the satellite-derived concentration of total suspended solids (TSS) has also been used to study the Mackenzie River plume in the Beaufort Sea (Doxaran et al., 2012), but with significant drawbacks, such as the inability of optical signatures to distinguish between river-derived sediments and resuspended bed sediments, and between vertical settling and horizontal dispersion of sediments.

River plume dynamics in Hudson Bay (HB) are poorly understood, yet essential to understanding the impact of climate change and river regulation on the HB marine system. Given the significance of cumulative fluvial discharges into HB (710 km³ yr⁻¹, including James Bay, 1984–2013; Déry et al., 2016; McCullough et al., 2019), knowledge of river plume dynamics is essential to understanding basin-scale oceanographic processes and associated ecological regimes. Because HB remains ice-covered for more than half the year from the winter to early summer periods of high runoff, previous studies have focused on the spread of river plumes under sea ice. Such plume dispersion can be more extensive than in open water due to reduced mixing by winds and waves (e.g., Messier et al., 1986; Ingram and Larouche, 1987; Li and Ingram, 2007; Peck et al., 2022). However, the use of optical satellite observations for river plume detection is not possible when sea ice is present in significant amounts. Therefore, to date, satellite remote sensing techniques have not been used to investigate river plumes in HB.

Here, we study the dispersion of freshwater from the Nelson and Hayes Rivers (NHR) in southwestern HB during the ice-free period using a satellite-based approach. The runoff from Nelson River (NR), augmented by Hayes River (HR), supplies more than half the freshwater loading from the southwestern watershed of HB (Déry et al., 2011) and is a significant source of terrestrial dissolved organic matter and particulates to the region (Guéguen et al., 2011; Guéguen et al., 2016; Capelle et al., 2020). The objectives of this study are three-fold:

• Evaluate the retrieval of CDOM and TSS from satellite imagery of the NHR estuary using available in-situ data;

• Develop a novel satellite-based method to evaluate the dispersal of CDOM and TSS-rich NHR plume into HB;

• Assess NHR plume extent and directionality into offshore HB using retrieved CDOM from satellite images.

The NR and HR flow northward through Manitoba (Canada) into a shared estuary in southwestern HB (Figure 1). The NHR estuary is a seasonally ice-covered, mesotidal estuary with a characteristic funnel shape (Wang et al., 2012). With a mean annual discharge of 109 km³ (in 1984–2013; Déry et al., 2016), the NR is one of the largest rivers draining into HB. It drains a watershed of 1.1 million km² across central and western Canada (inset in Figure 1). It transports significant amounts of suspended sediments and DOM from the Canadian Shield and the Canadian Interior Plains into HB (Rosenberg et al., 2005). Large lakes, such as Lake Winnipeg along the river course, and coastal mudflats act as both sinks and sources of sediments and terrestrial and algal CDOM (Goharrokhi et al., 2022). Since 1978, flow diverted from Churchill River has contributed 24 km³ yr⁻¹ to the flow in the lower NR (i.e., one quarter of the total discharge). Also, since 1978, the level of 23,750 km² of Lake Winnipeg has been regulated to provide seasonal storage (State of Lake Winnipeg, 2020), allowing for increased river discharge during the winter (Déry and Wood, 2005), reducing the spring and summer peak inflow into the lake (McCullough et al., 2019). Furthermore, over the post-diversion and regulation period, the discharge of the NR increased by about 42%, from 92 km³ yr⁻¹ in the decade 1984–1993 to 131 km³ yr⁻¹ during 2004–2013 (Déry et al., 2016), due to increased runoff driven by climate change in the upper watershed (Clair et al., 1998).

The HR, originating from Northern Manitoba, drains an entirely Precambrian Shield watershed into the same estuary as the NR (inset in Figure 1; Rosenberg et al., 2005). It is an unregulated river with an annual mean discharge of 19 km³ yr⁻¹ during 1984–2013 (Déry et al., 2016). The HR discharge has also increased over the last three decades, but by a smaller amount (28%) compared to NR (Déry et al., 2016). Over the same three decades, the two rivers accounted for 61% of the 210 km³ yr⁻¹ discharge into HB between the mouth of Churchill River and the entrance to James Bay (table 2 in McCullough et al., 2019). Since it drains through the northern wetlands and peatlands, HR has a higher CDOM and lower TSS than NR, which is readily visible in satellite imagery (Figure 2). Other significant rivers draining into the southwestern HB are the Churchill and Seal Rivers located north of NR. The combined mean annual discharge volume of these two rivers was determined to be 35% of the NHR discharge for the same 1984–2013 period (Déry et al., 2016), thus highlighting the significance of NHR discharge into southwestern HB. A geostrophic anticlockwise circulation regime prevails off the southwest coast, transporting a portion of terrestrial freshwater toward eastern HB (Prinsenberg, 1986; Ridenour et al., 2019).

Surface water samples for CDOM and TSS analysis were collected from 40 stations in the NHR estuary onboard a barge launched from the M/V Strait Signet from August 7 to 12, 2006 (except August 11, 2006; Figure 2; Table S1). Water samples filtered through a 0.7-μm Whatman 25-mm glass fiber filter were scanned for CDOM spectral absorption in the 275–800 nm wavelength range using an Ocean Optics S2000 spectrophotometer equipped with 10-cm long quartz cuvettes. Light transmittance (T) through the samples was converted to absorbance spectra (A) using Equation 1.

After baseline correction at 685 nm (Babin et al., 2003), the raw absorbance spectra were converted to spectral absorption coefficients using the following equation:

where a_CDOM is the CDOM spectral absorption coefficient (m⁻¹) at a wavelength λ (nm). A_CDOM is the CDOM absorbance (unitless fraction) and d is the path length (m) of the cuvette (0.1 m) used for the spectrophotometric measurements. Finally, the CDOM absorption spectra were modeled using Equation 3 to obtain the spectral absorption coefficient at 412 nm, a_CDOM (412), and the exponential slope, S, between 300 and 600 nm (Bricaud et al., 1981).

The CDOM absorption at the 412-nm reference wavelength was used in this study to measure CDOM when describing water properties.

TSS (mg L⁻¹) was determined by filtering a known volume of water sample through pre-weighed 0.7-μm Whatman 47-mm diameter glass fiber filters. After filtration, the filters were rinsed three times with Milli-Q water to ensure the removal of salt particles. The processed filter papers were oven-dried at 105°C for 24 hours. After drying, the filters were placed in a desiccator to cool to room temperature and weighed. This process was repeated three times to ensure a complete loss of moisture content. The difference between pre- and post-filtering of dry weights divided by the filtered water volume provided the TSS (Van der Linde, 1998).

These in-situ data were used to determine empirical relationships between CDOM and TSS and water surface remote sensing reflectance (Section 2.3) to enable satellite-based retrieval for NHR runoff characterization.

In-situ spectral remote sensing reflectance, R_rs(λ), the ratio of outgoing spectral radiance to incident spectral solar irradiance, was measured at 12 stations on August 8–9, 2006 (Table S1) using an Analytical Spectral Device FieldSpec 4. Surface water was sampled simultaneously for water quality analysis. Optical sampling was carried out at 1.4-nm bandwidth continuously through the wavelength range of 400–900 nm. The incident irradiance was measured using the reflectance of a calibrated Spectralon panel (placed horizontally) with known reflectance properties. Radiance was recorded at 2- or 4-second intervals using a 15° fore-optic positioned vertically downwards, about 2 m from the barge and 2 m above water, switching between over-water and over-panel (Spectralon) at five observation intervals. The R_rs spectrum calculated for each station is the median of at least 25 over-water observations from this series. These spectra were corrected for sun-glint following the similarity spectrum-based approach by Ruddick et al. (2006). They were the first and, so far, the only available in-situ remote sensing dataset collected for the NHR estuary region.

A blue (488 nm) to red (667 nm) R_rs band ratio (i.e., MODIS bands 10 and 13) was used to formulate an empirical relationship with in-situ CDOM absorption at 412 nm using a non-linear regression approach (Doxaran et al., 2005; Campanelli et al., 2017):

Following Nechad et al. (2010), non-linear empirical relationships were established for TSS retrieval using a single red band of R_rs 667 nm:

The fits of the data to these two relationships are shown in Figure 3a and b, respectively. These simple empirical relationships (Equations 4 and 5) were used to retrieve CDOM and TSS from geolocated Moderate Resolution Imaging Spectroradiometer (MODIS) images (on NASA’s Earth Observing System Aqua and Terra satellites). Level 1A radiance to R_rs conversion was carried out using the Management Unit of the North Sea Mathematical Models (MUMM) atmospheric correction module (Ruddick et al., 2000). This conversion was performed using the l2-gen processor of NASA’s SeaDAS image processing platform. The cloud-albedo threshold of 0.40 was used in the l2-gen processor instead of the default 0.20 to avoid flagging of highly reflective turbid waters of the NHR estuary as cloud pixels (Matsuoka et al., 2016). A match-up analysis, centered on station locations, was performed to ascertain the retrieval efficacy of the empirical relationships, resulting in 15 match-up points with a <2-hour difference between in-situ sampling and the satellite overpasses.

Thirty-three clear-sky and ice-free MODIS scenes obtained for August over the 2000–2010 period (Table S2) were used in this study for the NHR plume characterization. Clear-sky and ice-free scenes were found for the month of August of each year, apart from 2004 (not included among the 33 images), when mobile ice floes were present at the NR estuary during available cloud-free scenes. As Equations 4 and 5 were formulated based on August in-situ R_rs only, this study focuses exclusively on this month to limit uncertainty due to variation in the solar incidence angle, the presence of sea ice, and the paucity of in-situ observations.

Distributions of CDOM and TSS were derived from satellite imagery within an area extending 500 km offshore from the undiluted river water in the NR mouth (57.07°N, 92.51°W; Figure 1; dashed line in Figure 2). The salinity within this region ranges from 0 to approximately 29 based on field observations (Granskog et al., 2007; Granskog et al., 2009; Granskog et al., 2012; Guéguen et al., 2016). With the river mouth as the point of origin, the study domain was divided into circular sections 10 km wide (Figure 1). The initial section was set to 50 km from the river mouth to avoid the hydrodynamic complexity near the river mouth and a data gap in the available satellite images. For each section, quantile concentrations for probability ranging from 5% to 95% (Q_0.05–Q_0.95) of CDOM and TSS were retrieved (e.g., Figure 4).

Built on the work of Gangloff et al. (2017) on turbid plume thresholding, this study extends the use of Q_0.95 asymptotes (red lines in Figure 4) to partition the river and marine surface waters of southwestern HB using a quantile range. The Q_0.95 asymptote (hereafter “threshold”) represents the concentration of the most CDOM- or TSS-rich offshore waters of southwestern HB. Pixels with CDOM or TSS below this threshold are considered offshore marine waters with less than 5% possible river influence. Water with CDOM or TSS above the threshold is defined as belonging to the riverine coastal domain.

We used the threshold to estimate the dilution of NHR discharge as a function of distance from the river mouth. The quantile values equalling the threshold for each circular section were determined for CDOM and TSS. For example, in Figure 4c, the threshold line crosses the Q_0.50 values in the section located approximately 110 km from the river mouth. Thus, at 110 km from the NHR mouth, 50% of CDOM exceeds the threshold (i.e., $0.500.95$ or 0.53 marine water fraction). Therefore, we calculate a river-influenced water fraction of $1−0.500.95$ or 0.53 at 110 km. Similar calculations were performed for all semi-circular sections within the CDOM and TSS imagery to model the dispersion of river CDOM/TSS-influenced waters. These imageries corresponded to three daily discharge conditions: low discharge (LD; 1500–4200 m³ s⁻¹), mid discharge (MD; >4200–6500 m³ s⁻¹), high discharge (HD; >6500–9000 m³ s⁻¹), and tidal stages (flood and ebb). The NHR discharge data were obtained from Manitoba Hydro at the Long Spruce generating station. Tidal stages in MODIS images were determined from the Canadian Hydrographic Service’s modeled water level data. Based on a current speed of 3 m s⁻¹ for NR (Wang et al., 2012), the water discharged from Long Spruce station was estimated to take approximately 15 hours to flow the 162-km distance to the NHR estuary. Therefore a minimum of 2-day-averaged discharge of NHR was used for correlation studies involving the river-influenced water fraction (0.50, 0.25, and 0.00) and their corresponding $aCDOMQ0.95$ and TSS.

Finally, the NHR runoff directionality was derived in two steps. First, the bearing of each pixel (within the study domain) from the NR mouth was determined. The weighted mean bearing for each semi-circular section was calculated from these values after binning at 3° within each section. CDOM representing its riverine fraction, i.e., corrected for the marine background using Q_0.95 threshold (Equation 6), was used as the weight parameter. Negative values of $acdomriver(412)$ indicated pure marine pixels were precluded as weights for estimating the NHR plume direction.

Higher values of CDOM (a_CDOM (412), range of 1.07–7.69 m⁻¹) were measured in the estuarine samples near the HR mouth compared to the NR estuary (0.72–5.99 m⁻¹; Figure 5a). The HR drains through peat bogs of Northern Manitoba, which are a source of high CDOM (Rosenberg et al., 2005). At the confluence zone of the rivers (salinity of 19–20), a_CDOM (412) dropped to 1.04–1.5 m⁻¹, while a minimum of 0.72 m⁻¹ was observed beyond a salinity of 20 (coastal waters). These CDOM trends along the salinity gradient agreed with Guéguen et al. (2011). This study also reported that HR samples had 1.19 times more CDOM (at 355 nm) than the NR estuary. Coastal waters of HB have been reported to have a mean CDOM (at 355 nm) of 1.81 m⁻¹ at a salinity of 27, with values around 2 m⁻¹ for southern HB (Granskog et al., 2007), which are equivalent to 0.61 m⁻¹ and 0.68 m⁻¹ (calculated using Equation 3 and slope S_250–400 = 0.019 nm⁻¹) at 412 nm, respectively. These values are consistent with our minimum observed in the coastal waters of southwestern HB. Xi et al. (2013) reported HB mean summer and fall a_CDOM (412) to be around 0.24 m⁻¹ and 0.25 m⁻¹, respectively (calculated using S_350–500 = 0.019 nm⁻¹ and 0.0174 nm⁻¹). For comparison, the Mackenzie River, draining into the southern Beaufort Sea, had a_CDOM (412) values within the range of 0.03–1.89 m⁻¹ (calculated using S_400–800 = 0.02 nm⁻¹; Matsuoka et al., 2011). It contributed to 94% of the total non-water absorption at 412 nm (Bélanger et al., 2008). Such a range of CDOM absorption was also observed in European coastal waters (Babin et al., 2003). Therefore, coastal a_CDOM (412) (0.72–7.69 m⁻¹; S_300–600 = 0.0176 nm⁻¹) of southwestern HB was comparable to that of the Arctic and European coastal waters, highlighting the regional and global applicability of NHR runoff studies.

The mean (± standard deviation) TSS was higher for NR estuary samples (12.42 ± 13 mg L⁻¹, n = 26) than in the HR (6.67 ± 3.71 mg L⁻¹, n = 14), with a maximum of 67 mg L⁻¹ and 16 mg L⁻¹, respectively (Figure 5b). Around the river confluence zone and beyond, a mean TSS of 9 mg L⁻¹ and 5 mg L⁻¹ was observed for NR and HR, respectively. Compared to the HR, which has an entirely Shield watershed, the Nelson is a much larger river (about 6× greater discharge) with a watershed extending beyond the Shield into the sedimentary Interior Plains. Moreover, it remains impacted by bank erosion due to large impoundments and diversion projects supporting hydroelectric generation (Environment Canada and Department of Fisheries and Oceans, 1992; Rosenberg et al., 2005; Goharrokhi et al., 2022). Published TSS measurements in HB are scarce. Xi et al. (2014) reported that the highest inorganic particle load in HB was observed in coastal waters near the NR estuary. TSS maxima of 2.5 mg L⁻¹ and 1.8 mg L⁻¹ were observed in HB for September 2005 and July 2010, respectively, with their mean being 0.74 mg L⁻¹ and 0.55 mg L⁻¹ (Xi et al., 2013). However, this study excluded the estuarine TSS values of the NR. The TSS minimum (3 mg L⁻¹) measured in the coastal waters of southwestern HB closely matches the fall maximum reported by Xi et al. (2013). Seasonal variability of TSS inside the NR estuary peaks during the spring freshet, with turbidity measurements of 500 NTU (272 mg L⁻¹), decreasing in late summer to 200 NTU (88 mg L⁻¹; Wang et al., 2012). The turbidity equivalent TSS was calculated empirically using a linear slope of 0.61 mg L⁻¹ NTU⁻¹ and a bias of 35 mg L⁻¹ (personal communication from Greg McCullough).

Both CDOM and TSS were diluted at increasing salinity (Figure 5). While there were enough samples (n = 26) to suggest the power-law decay of NR CDOM (R² = 0.81), the data points from HR showed a linear decay (R² = 0.93; Figure 5a). (We use “decay” in its mathematical sense, and not to imply chemical decay.) Though NR had a lower CDOM relative to HR, it decayed relatively faster (0.72 m⁻¹ salinity⁻¹) than HR CDOM (0.39 m⁻¹ salinity⁻¹). This enhanced decay rate might be attributed to the higher hydrodynamic forcing of the larger NR. Guéguen et al. (2016) found that CDOM photodegradation decreases with increasing salinity. Therefore, in the context of this study, CDOM in the salinity range >20 will be least affected by photodegradation and thus can serve as a conservative tracer of river water. Such conservative mixing of CDOM in the coastal waters of HB was also reported by Granskog et al. (2007).

Similarly, for TSS, a higher decay rate was observed for NR (0.83 mg L⁻¹ salinity⁻¹; R² = 0.58) than HR (0.34 mg L⁻¹ salinity⁻¹; R² = 0.33; Figure 5b). The higher standard deviation of TSS samples suggests a higher variability in the particulate domain, thus non-conforming to a single source. Well known is that suspended solids are not generally conservative tracers of freshwater flowing into marine waters. They will tend to settle out as sediment transport capacity decreases, and, particularly in marine estuaries, the finer particles will tend to flocculate, forming larger particles marked by higher settling velocities (e.g., Kranck, 1984; Burban et al., 1989; Manning et al., 2010; Bainbridge et al., 2012; Many et al., 2016). Moreover, sediment dynamics in the Nelson and Hayes estuary are complicated by wave and tidal erosion of widespread mudflats, which provide an additional source brought into suspension by turbulence in the fluvial-marine mixing zone (inset in Figure 5b). Hence the non-conservative nature of TSS decay with increasing salinity in the NHR estuary (Figure 5b) limits the value of TSS as a river runoff tracer in our study area.

Because of reduced ice cover in HB, which led to higher evaporation of surface waters, an increase in cloud cover over southern HB was recorded for 1998–2018 (Fazel-Rastgar, 2020). The frequent presence of clouds poses a significant challenge to obtaining cloud-free optical satellite data and matchups with in-situ data for our study area in late summer. Nevertheless, of the 40 in-situ data points collected on August 7–12, 2006, 28 data points were available (excluding the 12 data points of August 8 and 9 used to formulate empirical relationships) for retrieval assessment of the empirical optical relationships (Equations 4 and 5). A single clear-sky MODIS satellite image was obtained for each day between August 7 and 10, 2006 (no sampling occurred on August 11, and no cloud-free images were obtained for August 12; Table S2). Among the 21 in-situ data points (excluding the 7 data points of August 12, 2006), stations D5, D6, and D7 of August 7, 2006, and B10, B11, and B12 of August 8, 2006, were precluded from the matchup analysis, as these points were sampled beyond the ±2-hour time frame of satellite data acquisition. The choice of this timeframe for matchup analysis was based on the study by Bailey and Werdell (2006) and accounted for the approximate 12-hour tidal period observed for NR. Therefore, 15 data points were effectively used for optical algorithm validation.

The satellite-retrieved a_CDOM (412) (⁠$aCDOMsat(412)$⁠; Equation 3) agreed well with the in-situ observations (R² = 0.74; slope = 0.60; Figure 6a), despite the interfering presence of high TSS. The root mean square error (RMSE) of 1.11 m⁻¹ can account for the intercept of 0.51 m⁻¹. The linear statistical relationship for TSS (Equation 4) was weaker than for a_CDOM (412) with R² of 0.69 and slope of 0.23, possibly due to its higher spatial-temporal variability caused by estuarine TSS dynamics (settling, resuspension, lateral flow; Figure 6b). The high y-intercept of 7.13 mg L⁻¹ can account for a higher RMSE of 12.34 mg L⁻¹, which may be due to the low variability within the TSS training data set (not subject to spatial-temporal misfits between optical and water quality data). Also, some combination of spatial resolution bias of satellite sensors, uncertainty in atmospheric correction approaches (Dorji and Fearns, 2017), and the TSS retrieval algorithm itself (Balasubramanian et al., 2020) might cause retrieval errors.

The river CDOM-influenced water fraction (hereafter RCWF), estimated from $aCDOMsat(412)$⁠, based on the quantile statistics (Section 2.4) obtained for each circular section (Figure 1), decreased exponentially with increasing distance from the NR mouth (Figure 7a). The exponential decay constants (b_CDOM; km⁻¹) were calculated using the following equation:

where d is the distance (km) from the NR mouth. The modeled b_CDOM in Equation 7 correlated negatively with daily discharge ranging between 1500 and 9000 m³ s⁻¹ (Figure 8a), simply showing that the NHR plume expanded offshore with increasing discharge. Therefore, instances with higher discharge were marked by extended offshore limits relative to the lower discharges shown by the black contours in Figure 9a, c, and e. For all 33 images, the mean extent of 0.50 and 0.25 RCWF (i.e., at 50% and 75% river CDOM dilution) were found to be at 88 km and 125 km, respectively (Table 1).

Moreover, the offshore depletion of CDOM covaried positively with NHR discharge (Figure 8b), mainly observed for the 50% and 75% attenuation of the river CDOM values (e.g., Figure 4a–c). This trend showed that NHR discharge contributes significantly to CDOM in southwestern HB. However, only 13% of the variability of a_CDOM (412) could be explained by b_CDOM variability (Figure 8c). As b_CDOM and a_CDOM (412) were functions of NHR discharge, their discharge-based normalizations were warranted. The variability of discharge-normalized CDOM (⁠$aCDOM′(412)=aCDOM(412)Q$⁠; s m⁻⁴) in ROFI accounted for 67% variance (p < 0.001) of discharge-normalized b_CDOM (⁠$bCDOM′=bCDOMQ$⁠; s km⁻¹ m⁻³; Figure 8d). Furthermore, together with NHR discharge (Q), $aCDOM′(412)$ variability explained about 79% of $bCDOM′$variance (p < 0.001; Equation 8).

where α and β are regression coefficients of Q and $aCDOM′(412),$ respectively, and ε is the constant. In a simplified way, $aCDOM′(412)$ is an inverse measure of CDOM load (discharge volume × CDOM). Therefore, the offshore decrease of CDOM is influenced by both CDOM load and hydrodynamic forcing of NHR discharge. Thus, from Equation 8, one can say that high NHR discharges will substantially load the coastal waters of southwestern HB with an enhanced volume of CDOM. This river water gets dispersed relatively farther into the bay, which was observed from the reduced decay rate of CDOM (b_CDOM) and characterized by a high Q_0.95 CDOM threshold (offshore CDOM maximum; Figure 9e).

Similarly, modeled b_CDOM for flood and ebb tidal flows correlated negatively with Q (Figure 10a and b) and vice versa for CDOM (Figure 10c and d). An enhanced rate of b_CDOM was obtained for flood tide (FT) relative to ebb tidal (ET) flow, resulting in reduced extents of river CDOM-influenced waters in coastal HB (Table 1; Figure 9g). For 100% dilution, ET dilution rate and CDOM correlated strongly with NHR discharge (Figure 10a and c) and vice versa for FT flows (Figure 10b and d). The tidal stage of the images, represented in terms of hours remaining until the next high tide, was not a statistically significant parameter in describing the variability of $bCDOM′$⁠, possibly an artifact of spring and neap water level variation influencing CDOM dilution.

The turbid water fraction (hereafter TWF) in southwestern HB, estimated from TSS^sat, also decayed exponentially with increasing distance from NR mouth (Figure 7b).

However, unlike b_CDOM, b_TSS was poorly and non-significantly correlated with Q, except for 0.25 TWF, which was statistically significant (Figure 11a). This trend was reflected in the invariable TWF extents estimated for the different discharge conditions (Figure 9b, d, and f; Table 1). Similarly, TSS was also poorly and non-significantly correlated with the NHR discharge (Figure 11d). Therefore, the traced TSS was not exclusively of riverine origin. Significant suspension of sediment is known to occur within the NHR estuary due to shallow bathymetry and the presence of mudflats (Figure 5b; Wang et al., 2012). Also, finer particles get flocculated to larger particles as the river water enters the marine domain with relatively higher salinity, leading towards an enhanced settling rate, as observed in Figure 5b. These factors contributed to the poor and non-significant correlation with river discharge, which was also observed for ebb (Figure 11b and e) and flood tides (Figure 11c and f). Reduced turbid plume extents and lower TSS were estimated for flood tides, unlike CDOM (Figure 9h; Table 1).

Once the NHR plume extents have been calculated, providing a spatial estimate of ROFI, determination of their flow direction seeks to understand river water transport pathways in southwestern HB. The NHR river waters propagated offshore with a mean azimuth of 54° from the river mouth (Figure 12a), that is, closely following the southern shore of HB. However, an overall range of 3–85° was observed within $aCDOMriver(412)$ azimuths for the semi-circular sections (Figure 12b). This range suggests that NHR waters often flowed straight into the bay and also propagated southeast towards James Bay. Thus, 0.25 and 0.50 RCWF, marked by higher a_CDOM(412), aligned along the southern coast of HB, while the NHR plume threshold was relatively offshore (Figure 12c). Therefore, the role of winds was envisaged to explain this variability as a function of atmospheric vorticity (Dmitrenko et al., 2020).

The averages of the 33 NHR plume azimuth (Az) ranged from 36° to 67°. Northerly (184°) and anticyclonic (clockwise; relative vorticity (⁠$rv¯$⁠) = $−1.8e−6$ s⁻¹) winds, obtained from ERA5, prevailed in southwest HB. The zonal and meridional winds were predominantly easterly (⁠$u¯$ = −4.1 m s⁻¹) and southerly (⁠$v¯$ = 0.34 m s⁻¹) with a mean wind speed of 5 m s⁻¹. Three classes of rv were obtained for 33 images: strong anticyclonic wind (⁠$rv¯$ = $−5.3e−6$ s⁻¹), weak anticyclonic wind (⁠$rv¯$ = $−1.3e−6$ s^–1, and cyclonic wind (⁠$rv¯$ = $8.2e−7$ s⁻¹). Strong and weak anticyclonic winds were classified based on mean negative rv.

Table 2 shows that the most angular $Az¯$ was associated with cyclonic winds (Figure 12a). The cyclonic (anticlockwise) wind rotation generates a cyclonic oceanic gyre, which forces offshore waters towards land through Ekman pumping (Dmitrenko et al., 2021), forcing the NHR water eastward along the southern shore of HB. A strong significant (p < 0.001) relationship was obtained between Az and intensified clockwise wind rotation (Figure 13a). Such wind reversal reverses the flow direction of the oceanic gyre, carrying onshore waters towards central HB (Ridenour et al., 2019). The offshore displacement of NHR water lowers the plume azimuth (Table 2; Figure 12a). With the reduction in the intensity of anticyclonic wind rotation, the plume azimuth became relatively shoreward, though NHR water still made its way into central HB (Table 2; Figure 12a). Under the prevailing anticyclonic wind regime, this effect was observed from the mean plume azimuth of all 33 images (Figure 12a and b).

The shoreward plume deflection during cyclonic winds was associated with strengthened southerly meridional winds (Figure 13b), while stronger easterly zonal winds accompanied by higher wind speeds resulted in the offshore displacement of the NHR waters (Figure 13c and d; Table 2). Such exchange of river water with central HB was reported by St-Laurent et al. (2011) and Granskog et al. (2009). Also, wind-induced plume transport occurs beyond Cape Tatnam (approximately 110 km from NR mouth) with a riverine regime preceding it (Figure 12a), suggesting the role of coastal geometry in NHR plume propagation.

Dispersion of the Nelson and Hayes river plume from the Nelson River estuary into southwestern Hudson Bay was investigated using satellite-derived colored dissolved organic matter and total suspended solids. The quantile concentration of CDOM and TSS (Q_0.05–Q_0.95) declined exponentially with distance from the river mouth. Based on paired in-situ salinity and a_CDOM(412) observations, $aCDOMsat$ (412) and its dilution over the distance from the NR mouth correlated well to NHR discharge and river water fraction. So, within the NHR plume, the defined river-CDOM influenced water fraction decayed exponentially with distance. Higher river discharge was associated with higher $aCDOMsat$(412) in the inner plume region, an extended offshore boundary, and a higher offshore CDOM maximum. The river water was diluted by about 75% at 125 km, ranging ±25 km due to discharge variability and ±35 km during the greatest flood and ebb tides. Based on pixel bearings of maximum $aCDOMsat$(412) in each quantile, wind vorticity contributes significantly to the NHR plume direction in HB. Clockwise winds carry the NHR plume towards central HB, while anticlockwise winds force alongshore eastward transport. Such geostrophic river water transport occurs elsewhere in the Arctic (Proshutinsky et al., 2009; Giles et al., 2012; Morison et al., 2012; Carmack et al., 2015).

Unlike CDOM, TSS did not correlate significantly with NHR discharge. Hence the offshore TSS maximum corresponds to the turbid plume boundary distinct from the NHR CDOM and river water plume boundary. The turbid plume carries sediments sourced from nearshore mudflats and/or bank erosion, in addition to sediments from the NHR watershed, and is reduced by settling. That is, TSS in southwestern HB is not a conservative tracer of river water and, in any case, is not solely sourced from the NHR. Hence, TSS is disregarded as an optical tracer of NHR water dispersion. Nevertheless, because both CDOM and TSS share similar absorption characteristics across the visible spectrum, TSS distribution needs to be understood so that the retrieval of CDOM from satellite images can account for its signature.

Due to time constraints and logistical difficulties, field measurements of Arctic estuaries are limited to a few transects throughout the freshwater-marine salinity gradient. These measurements are influenced by river discharge, tidal, wind, and ocean background current variability (Horner-Devine et al., 2015). We have demonstrated a novel satellite-based quantile regression approach for analyzing complex two-dimensional river water dispersion into the ocean at an estuary scale, which addresses such issues. Arctic and sub-Arctic river runoff is expected to increase due to increasing precipitation driven by climate change (Gagnon and Gough, 2002; Peterson et al., 2002; Serreze et al., 2006; Haine et al., 2015; Carmack et al., 2016; McCrystall et al., 2021). Our plume mapping approach can compensate for the dearth of in-situ data and advance our understanding of the response of river plume dynamics to increasing freshwater discharge into the Arctic, sub-Arctic, North Atlantic, and European coastal waters. Further studies on CDOM persistence and residence time are necessary to support the integration of the time dimension into river plume tracking by CDOM dilution. Furthermore, due to the poor quality of brightness temperature of seas close to landmasses, satellite-based sea surface salinity products like ESA-Soil Moisture and Ocean Salinity (SMOS) are not well suited for estuaries and nearshore waters (Oliva et al., 2012; Martín-Neira et al., 2016). Thus, a salinity data gap in the land-ocean continuum prevents continuous and synoptic river plume tracking from space. By adopting this technique based on the CDOM and salinity relationship, salinity can potentially be mapped for nearshore waters, thereby addressing estuarine and coastal salinization and desalinization, informing coastal management. Space-based mapping of river water mixing with seawater aims to deliver on Ocean Decade Challenge 1 and 2 (https://oceandecade.org/challenges/), and Sustainable Development Goal 14 of the United Nations.

Primary remote sensing dataset: Level 1 MODIS Terra and Aqua were obtained from NASA EODIS LAADS DAAC.

Primary wind dataset: ERA5 reanalyzed hourly wind data were acquired from Copernicus ECMWF.

Secondary dataset: In-situ CDOM, TSS, salinity, remote sensing data, have been provided as supplemental material Excel files listed under Supplemental material.

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

Table S1. List of in-situ sampling stations, time and date of sampling, tidal flow stage, and optical application for algorithm development and matchup. Docx

Table S2. List of satellite data used for Nelson/Hayes River plume dispersion study, and discharge type, tidal flow, and wind regime. Docx

Empirical_Optical_Relationship_Data.xlsx

Rrs Data_Compiled.xlsx

In-situ_Validation_All.xlsx

This study is a contribution to the Arctic Science Partnership (ASP). The authors thank Manitoba Hydro for the 2006 Nelson estuary field campaign and Mats Granskog for collection and processing of CDOM samples. They thank the continued support of various members of the Centre for Earth Observation Science, University of Manitoba, where this research was conducted. We also thank anonymous reviewers, associate editor, and editor-in-chief whose comments contributed to changes which improved the quality of this manuscript. As a special mention the authors thank the late David Barber for his inspirational support for this research.

This study is a contribution to the Natural Sciences and Engineering Council of Canada (NSERC) Collaborative Research and Development project: BaySys (CRDPJ 470028-14).

The authors declare that they have no conflict of interest.

Contributed to conception of the article: AB.

Contributed to the design of the article: AB, JE, SB, GM.

Contributed to acquisition of the field data: GM.

Contributed to processing and analysis of the data: GM, AB, AM.

Contributed to the interpretation of the data: All authors.

Drafted the article: AB.

Contributed to revision of the article: AB, GM, JE, SB, DD, AM, KS.

Approved submission of the manuscript: All authors (except the late David Barber).

How to cite this article: Basu, A, McCullough, G, Bélanger, S, Mukhopadhyay, A, Doxaran, D, Sydor, K, Barber, D, Ehn, J. 2023. Plume dispersion from the Nelson and Hayes rivers into Hudson Bay using satellite remote sensing of CDOM and suspended sediment. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2022.00076

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

Associate Editor: Laurenz Thomsen, Department of Marine Sciences, University of Gothenburg, Gothenburg, Sweden

Knowledge Domain: Ocean Science

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