Quantification of Discharge-Specific Effects on Dissolved Organic Matter Export From Major Arctic Rivers From 1982 Through 2019

.

Supporting Information: Supporting Information may be found in the online version of this article.et al., 2006;Peterson et al., 2002).Furthermore, the Arctic Ocean is shallower and receives a much greater proportional volume of river water than other major ocean basins (McClelland et al., 2012) Thus, quantifying and understanding the connections between river hydrology, ocean biogeochemistry, and the coastal carbon cycle has implications throughout the Arctic.The widespread environmental changes are likely to alter multiple aspects of the Arctic carbon cycle, having implications for the region and global carbon budgets and carbon-climate feedbacks.
Arctic rivers export organic matter that includes a plethora of terrestrially sourced material, including an annual average flux of 34-38 Tg C of dissolved organic carbon (DOC) (Holmes et al., 2012).Approximately half of the annual total DOC export (18.1 Tg C) originates from the six largest rivers (listed in descending order by mean annual Q): the Yenisey, Lena, Ob', Mackenzie, Yukon and Kolyma (Holmes et al., 2012).Although the DOC flux only equates to ∼10% of the satellite estimated annual Arctic marine net primary production (Lewis et al., 2020), there are undoubtedly direct and indirect influences of riverine DOC export on the net carbon balance of the Arctic marine environment (Drake et al., 2018;Hessen et al., 2010;Terhaar et al., 2021).Q and DOC concentration are positively correlated in the six largest Arctic rivers (Holmes et al., 2012), and thus increasing annual Q is expected to be accompanied by increasing DOC export (the product of DOC concentration and Q).DOC-Q relationships may be affected by long-term permafrost thaw, vegetation shifts, and many other factors driven by climate change, but positive relationships between DOC concentration and water discharge found in the Arctic, and a majority of large rivers throughout the world (Raymond & Spencer, 2015), suggest that discharge will continue to serve as a key predictor of DOC concentrations in Arctic rivers in the future.
Accompanying a potential change in total DOC export (and concentration) with increased warming and river flow is the inherent molecular-level complexity of the total DOM pool that is fundamentally related to the terrestrial source and in-water processing from soils to sea (Amon et al., 2012).Chromophoric dissolved organic matter (CDOM), or the portion of DOM that absorbs UV-visible light, has been shown to be strongly linearly correlated with DOC concentration in Arctic rivers (Mann et al., 2016;Novak et al., 2022;Stedmon et al., 2011).CDOM spectra can offer insight into the bulk composition of an individual riverine sample such as molecular weight and aromaticity (Helms et al., 2008;Weishaar et al., 2003), as well as biomarkers such as the lignin content as an unambiguous tracer of vascular plant inputs (e.g., Mann et al., 2016;Spencer et al., 2008).In addition, the absorbance of light by CDOM can have widespread carbon cycle impacts because of the potential for photochemical reactions in Arctic rivers and the receiving coastal ocean (Grunert et al., 2021;Osburn et al., 2009;Ward et al., 2017), and by decreasing the amount of light available for phytoplankton growth (Clark et al., 2022).DOC and CDOM can potentially be transported far afield in the ocean, and the total annual export from the six largest rivers is ∼6 times greater (∼15 Tg C) than the annual particulate organic carbon (POC) flux (McClelland et al., 2016).CDOM also contributes to the heat budget of the Arctic with greater CDOM resulting in more rapid short-wave radiation absorption which increases sea surface temperatures (Granskog et al., 2015;Hill, 2008;Pegau, 2002).As the numerous complex and interrelated processes that govern the flux of terrestrial material into riverine systems continue to be affected by climate change impacts, the DOC-CDOM relationship may shift.Potential changes in DOC-CDOM relationships can complicate efforts to use remote sensing platforms that can estimate CDOM to predict the DOC concentration in Arctic inland and coastal waters (e.g., Griffin et al., 2018;Matsuoka et al., 2022;Novak et al., 2022).The relationship between DOC and CDOM optical properties will be influenced by source and biogeochemical processing but may also change in time due to shifting hydrology and ecosystem dynamics.Assessing how baseline changes in discharge (i.e., seasonality and volume) may affect DOC and CDOM in the absence of other confounding factors, such as shifting land and vegetation processes, is important for predicting future changes in Arctic river concentration and export with continued warming.
Here, a record of 329 observations of DOC concentration and CDOM absorbance collected from 2009 to 2019 by the Arctic Great Rivers Observatory (ArcticGRO) (Holmes et al., 2021a(Holmes et al., , 2021b)), coupled with river discharge records from the furthest downstream locations (Shiklomanov et al., 2021) were used to estimate DOC export for the six largest Arctic rivers, the "Big 6" (Figure 1).A method that utilized spectral information to derive chromophoric dissolved organic carbon (CDOC) mass and concentration was implemented to determine chromophoric versus total DOC concentration and export over 2009-2019 and hindcast backwards to 1982 utilizing the river discharge information and the predicted relationship between Q and concentration.Trends were calculated for each river and the sum total to assess potential discharge-linked changes in DOC and CDOC concentrations, optical properties, and export.DOC-Q relationships for the 2009-2019 period were used for hindcasting.Holding the DOC-Q relationship constant allowed us to assess the potential effects of changing hydrological seasonality

Study Sites and Methods
The six largest Arctic rivers (Figure 1) drain an area of 11.3 million km 2 and had a combined median annual freshwater flow of 2145 ± 120 km 3 yr −1 from 1982 to 2019, which is 41% of the total mean freshwater export to the Arctic region estimated for all combined watersheds with an area of 22.1 million km 2 (Feng et al., 2021).Between 1982 and 2019, 25 years had complete records for all six rivers; therefore, aggregated analysis was based on these 25 years, while analysis on individual rivers used years where all 365 days were measured.This approach prevents potential artifacts from filling in missing data via interpolation or extrapolation.DOC, CDOM and many other chemical properties were measured across all seasons under the ArcticGRO project, and the data and collection methodology are publicly available (Holmes et al., 2021a(Holmes et al., , 2021b)).
Concurrent measurements of DOC concentration and CDOM absorbance collected from 2009 to 2019 were used to estimate the CDOC concentration as a fraction of the total DOC pool, with the non-chromophoric DOC (NCDOC) as the remaining portion (Figure S1 in Supporting Information S1).An automated two-step process was used to filter out anomalous CDOM absorbance spectra, abs(λ), before converting to Naperian CDOM absorption (a(λ); m −1 ), which is detailed in the Supporting Information S1.Paired measurements of a(λ) and DOC concentration were combined to estimate CDOC specific absorption spectra (a*CDOC(λ); m 2 g C −1 ) for each season and each river.Seasonal observations were sorted into winter (November-April), spring (May-June) and summer (July-October) following Holmes et al. (2012).For each seasonal set (denoted by subscript s), ordinary least squares (OLS) linear regressions of DOC concentration as a function of absorption at 300 nm (a 300 ; m −1 ) were estimated using the MATLAB function fitlm.m (Equation 1) where the y-intercept (  NCDOC ) represented the DOC concentration where a 300 was 0.0 and operationally defined the mean NCDOC concentration, and α is the slope of the DOC:a 300 relationship (g C m −2 ).The individual initial CDOC (iCDOC s ) concentration was then estimated as the difference between each DOC measurement grouped by season (DOC s ) and the corresponding  NCDOC (Equation 2).For each river, this provided a range of individual estimated iCDOCs and three seasonally constant initial  NCDOC (Table S1 in Supporting Information S1).

DOC𝑠𝑠 = NCDOC𝑠𝑠 + 𝛼𝛼(𝑎𝑎300𝑠𝑠)
(1) The iCDOC s estimates and coincident measured a(λ) from ArcticGRO were then used to derive seasonally averaged   * CDOC() (m 2 gC −1 ), where the mean seasonal   * CDOC() was the average of each derived estimate grouped by season (Equation 3).Rather than directly using iCDOC as the operational estimate of CDOC for each seasonally grouped set of DOC and CDOM measurements in the river flux modeling, the spectral information contained in a(λ) was utilized to better constrain the contribution of different seasonal CDOC s (winter, spring, and summer) to the total CDOC concentration for each pair of a(λ) and DOC measurements.An iterative linear least square fitting routine utilizing the full a(λ) (270-700 nm) was employed that predicted the total known a(λ) as a function of the product of the average seasonal   * CDOC() estimated from the OLS regressions described above and the unknown seasonal CDOC s concentration (Equation 4) (Clark et al., 2019).This was accomplished using the MATLAB curve fitting function lsqnonneg.mand a three-term model to estimate individual concentrations of CDOC (g C m −3 ).Each CDOC estimate from the curve fitting routine was a combination of winter, spring, and summer CDOC operationally defined by the a*CDOC(λ) s .One-and two-component models were also tested, but the three-component model had the best performance as a predictor of the observed a(λ) (see Section 3.1).Each set of estimated CDOC concentrations at each time point was then summed as the total CDOC for each time point.Finally, individual NCDOC were recalculated as the difference between measured DOC and the predicted CDOC concentration from Equation 4 to yield a time varying NCDOC for each river.The measured DOC and estimated CDOC and NCDOC values were then assembled into input files for use in the United States Geological Survey (USGS) Load Estimator (LOADEST) river mass flux prediction software (Runkel et al., 2004).These concentration estimates can also be used for other types of mechanistic models that require concentration as a state variable to drive numerical equations of transport and biogeochemical reactions.
A 100-member ensemble modeling approach for each river was used to assess the variance in predicted loads with different sets of inputs of concentration and Q and to offer 100 separate estimates of predictions for each river.
For each model run, 80% of each set of measurements was randomly sub-sampled to generate the LOADEST estimate.The 100 sets of predictions were then used to provide a range of estimated mass and concentration for each river for each year.This methodology is similar to the bootstrap method, and the median values were similar to those produced using all measurements as the input for one LOADEST simulation.This method provides a range of potential values for each variable and year to show the variance of estimates through time.All measurements were taken after 2008 and therefore estimates of past concentration assume that the relationship between concentration and Q has not changed substantially from 1982 to 2019; any past estimated loads and emergent trends are therefore due to variation in the hydrographs (i.e., amount of water and seasonality of its flux).This method inherently excludes changes in concentration and Q due to environmental processes such as changing source concentration and composition or inputs of DOC with increased warming and vegetation greening.We leveraged this modeling constraint to isolate and evaluate expected changes in DOC fluxes that can be directly attributed to changes in quantities and seasonal distributions of Q.Our approach has the added benefit of being relatively simple to implement whereas more complex models that attempt to represent changing environmental drivers (that are likely unconstrained in the distant past) can rapidly become complex.Mann-Kendall trend analysis using the MATLAB file exchange package ktaub (Burkey, 2023) and Theil-Sen's slope regression were performed to estimate trends and statistical significance for each river and in total and Mann-Kendall p-values are reported for significance levels (see Supporting Information S1 for details on trend analysis and the specific configuration of the LOADEST model).

Inter-River Variation in DOC and CDOM, 2009-2019
The greatest average annual DOC watershed yield was 2.5 g C m −2 yr −1 in the Lena River Basin followed by the Yenisey and Yukon Rivers (Figure 1).There is substantial variation in measured DOC concentration and a(λ), but all rivers exhibited characteristic seasonality with the greatest DOC concentration and a(λ) in the spring and summer when flows were highest.Annually averaged a 300 was greatest in the Lena and the Ob' and lowest in the Mackenzie, while the DOC to a 300 ratio was highest in the Kolyma and Lena Rivers, followed by the Yukon (Table S1 in Supporting Information S1).There was a significant (p < 0.05) positive correlation between DOC concentration and a 300 for all rivers in all seasons (Figures 2a-2f), but the DOC: a 300 ratio (slope of each linear regression) and the non-chromophoric DOC (y-intercept) varied seasonally and between rivers (Table S1 in Supporting Information S1).On average, measured DOC and a 300 were 155 ± 44% and 218 ± 69% greater in the spring relative to the summer.Spring DOC and a 300 were 203 ± 53% and 429 ± 162% greater than winter DOC, indicating that the total DOC pool is strongly enriched in CDOM during the freshet season months when runoff is greatest (Figure S2 in Supporting Information S1).The greatest percent difference between mean seasonal DOC and a 300 was in the Kolyma (205% and 296% in DOC and a 300 between the spring and summer) and the Yukon (575% and 255% in DOC and a 300 between spring and winter).This is consistent with substantial evidence of predominantly surficial flow during spring freshet that transports DOC enriched in fresh aromatic rich material that is highly light absorbing (Holmes et al., 2012;Mann et al., 2016;Raymond et al., 2007;Spencer et al., 2009).
iCDOC ranged from 1.3 to 13.0 g C m −3 across the rivers, and iCDOC was greatest in the spring for each river but the Ob', which had a mean iCDOC of 8.8 g C m −3 in the summer (Table S1 in Supporting Information S1).When all rivers were combined, the winter, spring, and summer coefficient of determination (R 2 ) for DOC:a 300 was 0.87, 0.87 and 0.88, the slope was 0.18, 0.19, and 0.19 (g C m −2 ), and the  NCDOC was 2.1, 2.0, and 2.0 (g C m −3 ).Previous analysis of linear relationships of DOC and CDOM absorption at 350 nm estimated a greater  NCDOC (2.41 ± 0.23 g C m −3 ) and a steeper linear slope (0.34 ± 0.01) for all rivers (Mann et al., 2016).The lower NCDOC presented here follows logically with the use of the shorter cutoff wavelength (300 vs. 350 nm), which was chosen because it reduced the uncertainty in the linear least squares fitting routine by up to 42%.The large proportion of the total DOC that is CDOC is consistent with molecular level and biomarker evidence (Behnke et al., 2021;Spencer et al., 2008) and the strong relationship between DOC and a 300 is indicative of similar environmental control of DOC and CDOM source.
The mean seasonal mass-specific CDOC absorption,   * CDOC() , is a relative measure of the UV-Visible light absorption capacity of the chromophoric portion of the total DOC pool (Figures 2g-2l).A greater   * CDOC() implies there are relatively greater proportions of efficiently absorbing compounds and/or molecular structures, whereas a lower   * CDOC() indicates less efficient absorption by the CDOC pool.a*CDOC at 300 nm, a*CDOC(300), was greatest in the spring (6.15 ± 3.5 m 2 g C −1 ) on average and had the largest variance among rivers.If the Mackenzie is omitted due to the large variance and relatively poor DOC:a 300 linear relationship (Table S1 in Supporting Information S1), the mean a*CDOC(300) was 5.09 ± 1.52, 5.53 ± 1.73, and 5.41 ± 1.37 (m 2 g C −1 ) in the winter, spring, and summer.Emerging from the absorption spectrum analysis is the seasonal variance in S 275-295 , which is greatest in the winter and lowest in the spring for all rivers (Figures 2g-2l).S 275-295 is strongly negatively correlated with carbon normalized lignin yields in the six rivers (Mann et al., 2016) and lower S 275-295 corresponds to higher molecular weight DOC (Helms et al., 2008).The seasonal variance of each river's spectra and concentration qualitatively corresponds to the hydrophobic organic acid (HPOA; a measure of aromatic compounds) fraction and the specific UV absorption measured at 254 nm (SUVA 254 ) in 2009-2010 (Mann et al., 2016).The shallower S 275-295 in spring and the enrichment of CDOC is evidence of efficiently absorbing material derived from vascular plants.This suggests that, a priori, CDOC is qualitatively related to biomarkers such as lignin phenols and molecular level markers of terrestrial material such as condensed aromatics and polyphenolics (Behnke et al., 2021).
Although some river spectra are similar in magnitude and shape across seasons (e.g., Kolyma and Ob'), seasonally separating each set of spectra and DOC measurements to derive a*CDOC(λ) allowed for the most accurate mathematical reconstruction of a(λ) and derivation of the final CDOC concentration.Model residual analysis indicated that across the spectra and for each river, the three-component seasonal model was much better at recreating the full observed a(λ) spectra (Figure 3).Compared to a two-component model, the sum of the square of the residuals of estimated a 355 nm was reduced by 76-95%.Therefore, the representation of the CDOC pool as a combination of "winter base flow," "spring runoff" and "summer intermediate" was appropriate.CDOC made Figure 3. Spectral residual % error for the predicted aCDOM(λ) using Equation 4 from the main text. 1 CDOC refers to a single a*CDOC spectra used to predict the observed aCDOM(λ), 2 CDOC refers to the use of 2 a*CDOC spectra (winter and combined spring/summer) and 3 CDOC refers to the use of all 3 seasonal a*CDOC spectra to predict the observed aCDOM(λ).SSR is the sum of the square of the residuals at 355 nm (the lower the number, the better).
up most of the DOC pool in all the rivers except the Mackenzie, with median CDOC concentration comprising 74%, 73%, 43%, 69%, 56%, and 68% for the Kolyma, Lena, Mackenzie, Ob', Yenisey and Yukon Rivers (Figure 2, third row).The fraction of the total DOC that is chromophoric varies seasonally as well (Figure 2, third row inset) with all rivers but the Mackenzie exhibited a peak of CDOC around the spring freshet followed by variations through the summer and into the fall.The Lena and Yenisey had the least amount of seasonal variation, and Lena CDOC made up a larger portion of the total DOC across all months.There were prominent secondary peaks in the Lena and Yukon CDOC, which is due to the lower a*CDOC(λ) for the summer and fall for these rivers.This indicates that there is an elevated proportion of the total DOC that is relatively less enriched in highly absorbing CDOM.Overall, warmer months exhibit a greater total fraction of CDOC relative to NCDOC, corresponding to the seasonal variation of specific UV absorption at 254 nm, HPOA's, and other terrestrial molecules that are characterized as light absorbing (Behnke et al., 2021;Mann et al., 2016).
The Mackenzie is the only river where NCDOC is the majority of the total for more than half of the year (Figure 2, lower panel) and is the river with the largest relative open water area (7%) (Stedmon et al., 2011) and the second largest lake to watershed area ratio (Table 1) (Lehner et al., 2022;Linke et al., 2019).Over half of the watershed resides above the Great Slave Lake, which has a residence time of ∼16 years (Gibson et al., 2006).The longer residence time allows for more processing of CDOM and DOC during transit (Liu et al., 2022), allochthonous inputs from the lentic systems, and the potential for more photochemical degradation of CDOC into NCDOC (Osburn et al., 2009;Ward et al., 2017).The Mackenzie had the lowest lignin phenol concentration and carbon normalized lignin yields in 2009-2010 (Mann et al., 2016), but the   * CDOC() calculated here suggest a highly absorbing, small concentration of CDOC in the spring.The linear relationship of DOC and a 300 was the weakest in the Mackenzie, so further exploration of optical and molecular signatures of the CDOM pool is needed at the ArcticGRO station and along the river.

Effect of Increasing Annual Discharge
Annual river discharge (Q) increased significantly for the Big 6 from 1982 to 2019 by 5.46 km 3 yr −1 (0.25% yr −1 ; p = 0.01) equating to a median increase of 9.2% in Q (Figure 4).Annual Q increased for the Kolyma (0.70 km 3 yr −1 , 0.72% yr −1 ; p = 0.09), Lena (3.10 km 3 yr −1 , 0.54% yr −1 ; p = 0.02), and Ob' (1.36 km 3 yr −1 , 0.34%  (Lehner et al., 2022;Linke et al., 2019) yr −1 ; p = 0.06) Rivers, and while the Mackenzie (0.65 km 3 yr −1 ; p = 0.2), Yenisey (0.13 km 3 yr −1 ; p = 0.82), and Yukon (0.27 km 3 yr −1 ; p = 0.4) Rivers also had a positive trend in Q, the Theil-Sen's slope was statistically inconclusive (Figure 4).This estimate is slightly less than the estimate of 5.8 km 3 yr −1 (p = 0.01) for the Arctic Basin from 1975 to 2015 (Durocher et al., 2019), very similar to the Eurasian River estimate of 5.4 km 3 yr −1 estimated from 1965 to 2000 (McClelland et al., 2006), and half of the recent pan-Arctic estimate from combined model-satellite data products of 11.3 km 3 yr −1 from 1984 to 2018 (Feng et al., 2021).Feng et al. ( 2021) used a combined satellite-hydrological model methodology to estimate flow from small streams to the largest rivers and included 22.1 million km 2 of watershed (vs.11.3 million km 2 for the Big 6 watersheds), indicating that small stream flow is a substantial contributor to pan-Arctic river flow trends.The largest river basins in the system, the Lena and Ob', form 81.2% of the total increase of 5.46 km 3 yr −1 .Kolyma and Lena have the smallest lake volumes per watershed area (0.105 and 0.096 m 3 m −2 ) and historically have the coldest annual average air temperatures (−13.1° and −9.92°) (Table 1) (Lehner et al., 2022;Linke et al., 2019).All of the river watersheds have a significant surface warming trend across all seasons from 1950 to 2012, with the greatest rate of warming occurring in the winter, averaging 0.43 K per decade (GISTEMP, 2023;Lenssen et al., 2019).Kolyma and Lena have the greatest rate of May-June warming at a rate of 0.38 and 0.36 K per decade, which is 48% and 39% faster than the mean of 0.26 K per decade of the other rivers (Figure S3 in Supporting Information S1) (GISTEMP, 2023; Lenssen et al., 2019).
For each LOADEST ensemble, all models exhibited a mean percent error of <2% for the DOC predictions and estimations of DOC and CDOM absorption showed little bias (Figure S4 in Supporting Information S1).Trend analysis was conducted on the median annual values of each ensemble distribution for each river, and in total by summing the daily mass estimates to get annual mass flux (export) or taking each set of daily mass estimates and dividing by the daily river discharge to estimate the concentration.Our total annual river DOC mass flux had a median annual value of 18.4 Tg C yr −1 , which comprised 13.2 Tg C of CDOC (Table 2).Holmes et al. (2012) estimated an annual DOC flux of 18.1 Tg C yr −1 which shows that our model methodology, on aggregate, is consistent when the entire model period is analyzed and we have extended the estimated time period and included new observations.Our hindcast annual total DOC flux estimates showed a significant increasing trend from 1982 to 2019 with substantial interannual variability of DOC export (Figure 5a).The total annual DOC and CDOC mass flux increased by 49.2 and 38.3 Gg C yr −1 (0.27% yr −1 and 0.29% yr −1 ) (Table 2) and the maximum DOC and CDOC mass exports were 21.8 Tg C and 16.3 Tg C in 2002.The Lena and Ob' Rivers each had a statistically significant increase in DOC mass flux, with an annual trend of 39.1 and 20.4 Gg C yr −1 (Figures 5c and 5e).The annual percentage increase in CDOC was marginally greater than that of DOC in the Lena River (0.68% vs. 0.66%) but substantially greater in the Ob' River (0.57% vs. 0.49%), indicating overall enrichment in CDOC through time.This equates to an estimated increase in a 3.8% greater proportion of CDOC export from 1982 to 2019 for the Ob' River.All other Rivers but the Yenisey had a positive increase in DOC and CDOC, which contributed to the total Big 6 export trend.(Hollander et al., 2013).

Changing Discharge Enriches Chromophoric Dissolved Organic Carbon
Hindcast estimates from the Kolyma, Lena, Ob', Yenisey, and Yukon Rivers each show a pattern of enrichment of annual DOC export in CDOC as annual discharge increases, indicated by the darker colored CDOC:DOC ratio at higher annual discharge values (Figure 5).The Mackenzie has an opposite pattern where in lower flow years CDOC is more enriched relative to DOC.The enrichment in CDOC is less clear in the Lena, but there is a consistent pattern of CDOC enrichment during high river flow years that manifests in the Big 6 on aggregate (Figure 5a).There was also a drop in total river flow and hindcast CDOC:DOC in the Yenisey over the last decade (Figure 5f).This suggests that in general, more CDOC is exported during high river discharge years and especially during peak years (e.g., 2002, 2014, and 2015 in the Big 6).Therefore, as annual river discharge has also increased in the Big 6 (McClelland et al., 2006; this study), we should expect a priori relatively more CDOC export as time progresses in the absence of other drivers of change.Indeed, we see a small difference between the Big 6 annual CDOC mass flux trend (2.90% per decade) and DOC mass flux trend (2.67% per decade), but at the individual river level, it is more difficult to detect decoupling between CDOC and DOC mass flux trends through time.
Hindcast annual mean CDOC to DOC concentration ratio (CDOC:DOC) increased linearly as a function of annual water yield (discharge divided by watershed area) for each river across the time series (Figure 6).The Mackenzie is the exception where CDOC:DOC and annual water yield had a negative correlation indicative of CDOC dilution at high annual flow rather than enrichment.This may be a statistical artifact due to the poor CDOM:DOC linear relationship (Table S1 in Supporting Information S1) or sampling difficulty in capturing the freshet which could bias the model prediction.As discharge increased in each river, the fraction of the total DOC that was CDOC increased, and this relationship holds across the time series and during the observational period of 2009-2019.For the Big 6, there is a significant positive trend in hindcast DOC (Figure 7a) and CDOC (Figure 7b) concentration from 1982 to 2019, and CDOC concentration increased at a faster rate relative to DOC.Although the CDOC:DOC mass flux ratio did not significantly increase over the full modeling period (Figure 7c), there is a positive trend in the modeled CDOC:DOC concentration ratio since 1982 (Figure 7d).DOC concentration increased at a rate of 10.8 mg C m −3 yr −1 and CDOC increased at a rate of 9.6 mg C m −3 yr −1 which equates to relative rates of 0.18% and 0.25% yr −1 .This emergent enrichment of the hindcast DOC with CDOC is driven by statistically significant increases in the Kolyma, Lena, and Ob' DOC and CDOC concentration, with CDOC increasing at a more rapid rate in all three rivers.Total lignin phenols generally increased with higher runoff for the six rivers and spectral characteristics indicative of higher total absorption and aromatic molecular composition, such as S 275-295 and the slope ratio (S R ), positively correlate with Q and more surficial flow (Mann et al., 2016;Spencer et al., 2009).Molecular level composition is also highly seasonal (Behnke et al., 2021) with higher molecular weight terrestrial compounds indicative of vascular plant inputs more prevalent during warmer months and overprinted on a more degraded basal DOC pool that is ubiquitous across the big Arctic rivers.The positive correlation between the CDOC:DOC ratio and high river flow and the significant positive trend in the CDOC:DOC ratio over time, on aggregate, indicate that independent of non-discharge factors that may have altered the CDOC:DOC-Q relationship, changes in discharge would have pushed Arctic River DOC toward CDOC enrichment over the last four decades.
While these results address how changes in the quantity and seasonal distribution of river discharge translate into changes in mass fluxes and chromophoric content of DOC, it is helpful to consider how other factors may and concentration (mg C m −3 yr −1 ) and the CDOC:DOC ratio of each quantity.The p-values were calculated using a z-test that assumes the calculated trends come from a normal distribution with mean of 0 (no slope).(Hollander et al., 2013).The plots without regression lack a significant (p < 0.1) trend over the period analyzed.
also influence DOC export via changes to the discharge-concentration relationship that are not considered here.Earlier thawing of the active layer and into permafrost layers would allow for more subsurface flow earlier in the year not only potentially releasing carbon stores in the upper soil profile (e.g., Wild et al., 2019) but also altering the timing of when water passes through the landscape into the river channels.Release of POC and subsequent hydrolysis during riverine transport may be an additional in-water source of DOM that manifests downstream (Wild et al., 2019), but the rates of particulate degradation in the water column and potential release of DOM and DOM enriched in chromophoric components are poorly characterized.The varied Arctic landscape and the state factor of the landforms are directly related to the impact of permafrost thaw on the release of DOC and  other materials from permafrost into riverine systems (Tank et al., 2020).Thaw of the Yedoma that is widespread in Siberia is linked to the release of DOC, CDOM, and POC into stream networks (Tank et al., 2020), but much of the released DOM is highly biodegradable (Spencer et al., 2015).Indeed, the molecular signature of permafrost-derived DOC in the Kolyma River is rapidly lost downstream of inflow locations (Spencer et al., 2015).Alternatively, an increase in median river depth and more rapid transport once material reaches surface water and river networks because of increased water discharge volume may shield DOC and CDOM-rich water from degradation by limiting exposure to warm summer temperatures and sunlight, both important degradation factors (Grunert et al., 2021;Ward et al., 2017).
Shifting plant communities and increasing plant greening (e.g., Myers-Smith et al., 2020) can alter the hydrology of the system by changing evapotranspiration, and plants can act as a source of new DOC both directly through leaching and indirectly through senescence and detrital breakdown.There is evidence of increased plant greening throughout the Arctic (Berner et al., 2020) and increased productivity in delta regions such as the Yukon-Kuskokwim (Frost et al., 2021), which could export more DOC and potentially CDOM as more atmospheric carbon is fixed in the plant community.The enhanced greening and a shift from short tundra plant species to tall woody shrubs could also contribute an altered source of vascular plant derived DOC as these plants colonize and move into historically tundra regions (Mekonnen et al., 2021).There are many direct and indirect processes that a changing plant community has on the carbon flux from terrestrial systems into rivers and beyond and it requires a better mechanistic understanding and representation to characterize specific drivers of changes in DOC and CDOM in the past and into the future.Especially challenging is the collinearity of many of these processes, which necessarily requires higher-order statistical techniques and mechanistic models to tease apart.There are few significant trends in DOC and CDOM since the ArcticGRO sampling began in 2003, even though there has been noticeable greening as determined by satellite observations (Berner et al., 2020).Since 2003, only the Ob' River has a significant (p < 0.05) positive trend in Q (4.7 km 3 yr −1 ).This emphasizes that the long-term change in river flow volume and seasonality is driving the hindcast trends in DOC.Nonetheless, the manifestation of these large-scale ecosystem changes will likely impact the net organic carbon flux from terrestrial ecosystems into Arctic rivers.These ecosystem-driven changes in biogeochemistry would act in addition to the hydrologically driven trends as the timing of discharge has shifted and the annual volume of water has increased in Arctic Rivers.

Seasonal Effects
Hindcast Big 6 winter DOC and CDOC mass flux has significantly increased since 1982 by 12.3 and 9.0 Gg C yr −1 at a decadal rate of 6.4% and 7.8%, respectively (Figure 8a).All rivers except the Ob' and Yenisey had a significant increase in winter DOC and CDOC flux and the rate of increase in winter CDOC flux was greater than that of total DOC and Q (Figure 8a).Hindcast trends were greatest for the Kolyma with a decadal rate of increase of 21.2% and 22.5% for DOC and CDOC flux, followed by the Lena at a rate of 13.4% and 13.5%.The combined effects of increasing wintertime discharge and increasing hindcast concentration led to an overall greater trend in DOC and CDOC mass flux in the winter.The spring had a muted response with low statistical significance and a large decrease in spring Yenisey CDOC and DOC export nullifying the modest increases in springtime flux from the other rivers (Figure 8b).The summed Big 6 export in the summer increased at a rate of 30.9 Gg C yr −1 for DOC and 20.14 Gg C yr −1 for CDOC (Figure 8c).In November-April, the consistent inter-river patterns in direction and relative magnitude of DOC and CDOC flux are strong evidence that increases in wintertime flow were a strong driver of hindcast increases in annual DOC export.The Kolyma, Lena, and Mackenzie Rivers have seen a significant increase of 43.5%, 32.0%, and 12.8% in median winter DOC watershed yield (g C m −2 yr −1 ) from 1982 to 2000 to 2000-2019.The difference from 1982 to 2000 to 2000-2019 equates to a median DOC yield increase of 16.3% for the Big 6 when aggregated.There are similar changes in CDOC yield and CDOC concentration, indicating that shifts in discharge seasonality are a large driver of the hindcast annual increase in DOC and CDOC flux and concentration.
Measured winter river discharge has increased significantly between 1982 and 2019 for all rivers except the Yenisey, but the rate of increase is substantially less than that modeled for CDOC and DOC export over the same period (Figure 8a).Historically, increases in winter discharge have been linked to river flow regulation by hydroelectric dams (McClelland et al., 2004), but major dam construction and reservoir infilling on the rivers analyzed here was largely completed before 1982.In addition, the Yukon River has the third greatest rate of increase in winter discharge and the Yukon has virtually no dams that regulate flow (Table 1).This indicates that the presence of dams likely did not impact the seasonal trend analysis.The total winter discharge increased by 1.24 km 3 yr −1 at a rate of 3.4% per decade, while modeled DOC and CDOC flux increased by 6.4% and 7.8%.The faster relative rate of increase in DOC relative to river discharge is because of the positive correlation between discharge and DOC concentration, and the faster rate of increase in CDOC is due to the further enrichment of DOC with CDOC as flow increases (Figure 6).Winter has often been under sampled due to logistical challenges in these remote locations, and historically, characterizing the bulk of the flux that occurs in spring and summer was prioritized.The total annual median winter flow (369 km 3 ) is 17.2% of the total annual flow (2145 km 3 ) but accounted for a relatively outsized increase in the total annual DOC and CDOC export because of the large proportional changes in hydrology.The winter flow rates are becoming generally more elevated; therefore, there are more days where CDOC enriched DOC is being transported by the river water.A priority moving forward should be to concentrate research efforts in the shoulder months when increased river flow and associated increases in carbon export may occur due to shifting hydrology.

Conclusions
This analysis demonstrates how changes in the magnitude and seasonal distribution of river discharge might have influenced DOC and CDOC export from the six largest Arctic rivers.Model results show that DOC and CDOC export, as well as the proportion of CDOC to total DOC, can be driven by changing discharge.The trends reported herein do not account for other factors such as permafrost thaw and changes in vegetation that may also influence DOC mass fluxes and composition but do show the strong effects that hydrologic changes alone can have.Hindcast trends over the study period were most evident in the Kolyma and Lena Rivers, which led to an overall trend in the six river ensemble.The modeled increase in DOC was generally enriched in CDOC due to a higher discharge volume in winter, especially during shoulder months.This coincided with higher relative amounts of CDOC because the greater relative mean discharge is correlated with a higher CDOC:DOC ratio.There were no significant trends in discharge or associated DOC and CDOC fluxes during the observational period from 2009 to 2019; only when hindcast values driven by longer-term trends in river discharge were analyzed did long-term trends in DOC and CDOC flux emerge.Moving forward, efforts should be expanded to further include winter months in field campaigns, especially the shoulder months that are seeing the greatest proportional increase in river flow.The source and chemical nature and reactivity of the organic carbon in these systems need to be further characterized so that realistic mechanistic models can continue to be developed.If accurate, mechanistic models can explicitly link changing landscape characteristics and hydrology, all overlaid with rapid regional Arctic climate change, to altered organic carbon processes.These models must incorporate relationships between ecosystem processes, landscape state factors and varying landforms, and carbon export into rivers in order to accurately represent climate change driven shifts in the carbon cycle across the region.

•
Modeled dissolved organic carbon export was 18.4 Tg C yr −1 (median) from 1982 to 2019 for the six largest Arctic rivers • Proportional contributions of chromophoric to total dissolved organic carbon (CDOC & DOC) are positively correlated with water discharge • Increasing discharge and shifting seasonality, independent of other factors, would have increased CDOC and DOC export from 1982 to 2019

Figure 1 .
Figure 1.A map of the Arctic with each of the Great River watersheds shaded by the estimated mean dissolved organic carbon yield (g C m −2 yr −1 ) from 1982 to 2019.

Figure 2 .
Figure2.(a-f) Dissolved organic carbon (DOC) as a linear function of chromophoric dissolved organic matter (CDOM) absorption at 300 nm (a300) for each season (first row) where the color of each diamond is the runoff (daily discharge per watershed area), (g-l) the mean chromophoric DOC specific absorption spectra (m 2 gC −1 ) for each season with the spectral absorption slope (μm −1 ) from 275 to 295 nm in the legend for winter (W), spring (S), summer (SS), and total (T), and (third row) the total DOC, chromophoric DOC and non-chromophoric DOC concentration for each river.The inset plots in the third row show the fractional climatological contribution of chromophoric DOC and non-chromophoric DOC smoothed with a 30-day moving average filter.

Figure 4 .
Figure 4. Time series of the annual river discharge from the sum of the Big 6 rivers and each of the individual rivers from 1982 to 2019.The Big 6 total is from 1982 to 2017 because the Mackenzie River had incomplete gauge records in 2018-2019.Solid lines through each plot represent the Theil-Sens slope regression and the dashed lines represent the confidence interval of each regression(Hollander et al., 2013).
The trend is the change in dissolved and chromophoric dissolved organic carbon (DOC & CDOC) mass flux (Gg C yr −1 )

Figure 5 .
Figure 5. Modeled total dissolved organic carbon (DOC) mass flux for the Big 6 rivers summed over each year (a), and each individual river (b-g).The shaded region is the standard deviation of the 100-member ensemble model predictions for each year, and the color shading of each mark is the chromophoric dissolved organic carbon (CDOC) to DOC mass flux ratio.Note that each panel color scale is different.Solid lines through (a), (c), and (e) are the Theil-Sen's slope regression line fit and the dashed lines represent the confidence interval of each regression(Hollander et al., 2013).The plots without regression lack a significant (p < 0.1) trend over the period analyzed.

Figure 6 .
Figure6.Chromophoric dissolved organic carbon (CDOC) to dissolved organic carbon (DOC) concentration ratio for each of the 6 rivers as a function of watershed runoff (median annual discharge per watershed area) and the mean across rivers (Big 6).The solid black lines represent the ordinary least squares linear fit, R 2 is the coefficient of covariance and α is the slope of each regression.Every regression was significantly correlated with a p < 0.0001.

Figure 7 .
Figure 7. Dissolved organic carbon (DOC) and chromophoric dissolved organic carbon (CDOC) concentration over time (a and b), the CDOC:DOC mass flux ratio (c), and the CDOC:DOC concentration ratio (d).Solid lines through each time series represent the Theil-Sen's slope regression and the dashed lines represent the confidence interval of each regression (Hollander et al., 2013).(c) does not have a significant trend (p < 0.1) over time.The shaded region is the model ensemble prediction mean ±1 standard deviation.

Table 1
Watershed, Land cover, and River Characteristics From the HydroATLAS Global River Database for the Six Largest Arctic Rivers From the HydroATLAS Database