Flux and age of dissolved organic carbon exported to the Arctic Ocean: A carbon isotopic study of the five largest arctic rivers

2.1. Field Sampling

[6] The sampling sites for each river were Dudinka, Zhigansk, Salekhard, Tsiigehthchic and Pilot station for the Yenisey, Lena, Ob', Mackenzie and Yukon, respectively. All discharge measurements were from gauged sites. For the Yenisey and Lena proximate gauged sites (Igarka and Kyusyur, respectively) were used for discharge. In general the sampling protocols used in the PARTNERS Project were modeled on the USGS equal discharge increment (EDI) sampling protocol (http://pubs.usgs.gov/twri/). During the open water period (no ice cover) this calls for the collection of five depth-integrated samples across the river channel on each sampling date. The five depth-integrated samples were then mixed together to yield a single composite sample intended to be representative of the biogeochemical flux of the whole river on that date. During the ice-covered period, near-surface water was collected from a single hole in the ice, located midchannel. Procedures for sample processing were identical to those used during periods without ice cover.

[7] During the open water period, water samples from all five rivers were collected using 60-kg depth-integrating sampler (model US D-96) fitted with a Teflon nozzle and Teflon bags. For all rivers but the Yukon, the five depth-integrated samples were collected at equal intervals across the river channel, with the transit rate of the sampler through the water column being constant during each sampling date (Equal Width Increment method according to USGS terminology). Thus, at deep locations or in fast current, a full 3-L sample would be collected, whereas at shallower locations or in slow currents, the sample volume would be considerably less. On the Yukon River, the channel was divided into five sections which each section accounting for 20% of the total discharge (Equal Discharge Increment method according to USGS terminology). In that case, the sampler transit rate was adjusted such that a full 3-L sample was collected at each location.

[8] The elaborate depth-integrated sampling strategy employed during the open water period was designed particularly for accurately quantifying suspended sediment flux. This was necessary because of the uneven distribution of sediments through the water column. In contrast, we found no evidence of vertical heterogeneity in concentrations of dissolved constituents. Our most complete comparison of dissolved constituents in depth-integrated and near-surface samples was in summer 2003 when the sampling protocols were being established. On each river, depth-integrated and near-surface samples were collected at five locations across the river channel. Because H₂¹⁸O is a powerful tracer of water origin and should indicate if water column mixing was incomplete, we measured H₂¹⁸O on each of the samples. No significant differences were found between near-surface and depth-integrated samples. Moreover, on the Yenisey River where equipment problems meant that most samples were not depth integrated (see footnotes in Table 2), similar temporal heterogeneity and summer/winter contrasts are apparent in DOC concentrations and age. Therefore, while near-surface samples must be interpreted with caution for particulate constituents, we are confident that dissolved constituent concentrations are accurately reflected in our small number of winter near-surface samples and in the Yenisey where the depth-integrated samples were not always obtained.

[9] The five individual depth-integrated samples were then mixed in a 14-L churn, resulting in the single composite sample. Processing for DOC samples consisted of filtration through a precombusted (550°C for 3 hours) 47-mm-diameter quartz filter (Whatman QMA; 1 micron) into an acid-leached 500-mL polycarbonate bottle. Samples were frozen on site and later shipped frozen to Woods Hole, Massachusetts. DOC samples were then sent to Yale University for analysis. The US D-96 samplers used for the PARTNERS project where equipped with Teflon nozzles and Teflon sample bags to minimize contamination during the collection process. The churns for compositing the individual depth-integrated samples were also made of Teflon.

2.2. DOC and DOC Isotopes

[10] Dissolved organic carbon concentration was determined using high-temperature oxidation and UV oxidation methods. High-temperature oxidation was performed on a Shimadzu TOC analyzer. UV-oxidation was performed using a UV system that was directly connected to vacuum lines where the resulting CO₂ was cryogenically purified for isotopic analysis using methods described elsewhere [Raymond and Bauer, 2001a; Raymond et al., 2004]. Briefly, 120 mL of river water were placed into a clean quartz tube. The samples were acidified with 0.2 mL of ultra-high-purity 40% phosphoric acid. The samples were then sparged with UHP nitrogen to remove any inorganic carbon. Pure UHP O₂ was then sparged through the system to provide an oxidant for the UV oxidation of DOC. The sample was then oxidized with UV. The resulting CO₂ was transferred to the vacuum line and cryogenically purified. The purified CO₂ was analyzed for carbon isotopes at the National Ocean Sciences AMS at the Woods Hole Oceanographic Institution, or the University of Arizona's AMS facility.

[11] The agreement between the two oxidation methods was good (data not shown, r² = 0.95 slope = 1.03). In this manuscript the UV-oxidation results are reported for each river (Table 2). Over the course of this study, five duplicates (Table 2) were run for carbon concentration and isotopes. For these five samples the standard deviation was 0.27 mg L⁻¹ for DOC and 9‰ for Δ¹⁴C. The large standard deviation for Δ¹⁴C was mainly due to a Lena sample from 8/24, where the duplicate samples were +14 and +64‰. For the other four samples the average standard deviation was 3‰.

[12] Carbon export from the three stations were estimated using USGS LOADEST (Load Estimator) program (R. L. Runkel et al., Load estimator (LOADEST): A FORTRAN program for estimating constituent loads in streams and rivers, in US Geological Survey Techniques and Methods: Techniques and Methods Book 4, chap. A5, 2004, available at http://pubs.usgs.gov/tm/2005/tm4A5). LOADEST uses daily element concentration and flow data to establish relationships between flow and concentration and applies these relationships to the complete daily discharge record. DOC concentrations from Table 2 and daily discharge from the USGS for the Yukon, Arctic RIMS (rims.unh.edu) for the Eurasian rivers and the Water Survey of Canada's Hydrometric Database for the Mackenzie were used. For the Yukon River, additional DOC data from 2004 and 2005 that were collected by the USGS and available online (http://ak.water.usgs.gov/yukon) were utilized for the flux estimates.

3.1. Carbon Export

[13] A large export of DOC during ice-out was recently demonstrated for the Yukon and Kolyma rivers [Finlay et al., 2006; Striegl et al., 2007]. The other high-latitude rivers studied here also exhibit a large flux of DOC during the ice-out period. This is partially due to the large proportion of water that is discharged during the spring thaw (Figure 1). It is also due to the large increases in DOC concentration that occur during high flow (Figure 2). All watersheds showed a marked increase in DOC concentration with flow (Figure 2). Ice-out concentrations were ∼2–3× higher than winter concentration for the Lena, Yenisey, Mackenzie and Ob' rivers and ∼5× higher for the Yukon River (Figure 2).

[14] The additive effect of high discharge and elevated DOC concentrations during ice out result in a large proportion of DOC export in these rivers occurring in a small period of time. For the two month period following ice-out approximately 60, 62, 53, 39, and 63% of the annual DOC export occurs for the Yenisey, Lena, Ob', Mackenzie and Yukon rivers, respectively (Table 1). This is markedly different than temperate rivers. For the Mississippi River the two month period around the spring freshet contributes only ∼25% of the annual DOC export (Figure 3). In these high-latitude rivers the pulse also occurs later in the spring/early summer slightly before the Arctic summer. The high discharge and predominately cold river water temperatures likely cause this DOC to pass through coastal estuaries and the near coast and become broadly available in the offshore coastal zones during the Arctic summer. High DOC concentrations have been reported in the Arctic Ocean during May–June [Skoog et al., 2005] and these fluxes therefore have important implications to net metabolism and nutrient cycle of these coastal ocean systems [Skoog et al., 2005]. It should be noted that the more attenuated pulses in the Ob' and Mackenzie are likely due to a large dam and lake, respectively.

[15] Normalizing DOC export to watershed size indicates that a major control on DOC yield from these rivers is water throughput (Figure 4). According to this analysis, a 2× increase in water yield for these watersheds results in a ∼2.5× increase in DOC export. Thus as the hydrology of these watersheds continue to change in future years [Peterson et al., 2002; Walvoord and Striegl, 2007], large changes in DOC export may ensue, although changes in watershed temperature are also argued to be important to DOC export from these watersheds[Frey and Smith, 2005; Striegl et al., 2005] and may cause shifts in the relationship between water throughput and DOC yield.

[16] Annual river discharge for this two year study was close to the 20-year average and therefore the fluxes reported in Table 1 should approximate average fluxes. Compared to previous studies, these estimates are arguably the most direct because they represent a large number of samples taken directly at the mouth of the rivers and included spring-thaw samples. For the four rivers draining directly into the Arctic Ocean (Yenisey, Ob', Lena, and Mackenzie), the total annual DOC flux is ∼15 Tg. Except for the Lena, these estimates of DOC flux are on the high end of ranges published in recent pan-Arctic estimates [Dittmar and Kattner, 2003; Opsahl et al., 1999]. For the Lena our two years worth of data provide an estimate that is ∼20–40% larger than previous estimates. These yields are large compared to temperate rivers owing to the dynamics of the spring thaw. The Lena and Mississippi, for instance are similar in size and have very similar rates of annual discharge, yet the Lena exports ∼2.5× more DOC (Figure 3) owing to the sharp increase of concentration with flow. The relationship in Figure 4 allows for a first-order estimate of DOC export to the Arctic when coupled with estimates of total watershed area and river water export draining into the Arctic Ocean [McClelland et al., 2006]. For rivers draining directly into the Arctic Ocean the DOC flux is estimated to be ∼25 Tg yr⁻¹. If we include watersheds that drain into the Arctic Archipelago, the Hudson Bay region of the Arctic, and the Bering Strait [McClelland et al., 2006] the total DOC flux is 36 Tg yr⁻¹.

[17] Our estimates of 25–36 Tg of DOC, however, should be considered conservative. The annual average yield of DOC for the Lena River was ∼2.3 g C m² yr⁻¹. A recent study in high-latitude peat-rich subwatershed of the Lena, where approximately 60% of the Lena's water originates, reported yields of 4.75 g C m² yr⁻¹ [Suzuki et al., 2006]. Similar recent results were obtained for the Yukon River where yields increased dramatically after input from the high-peat Koyukuk River near the mouth of the Yukon [Striegl et al., 2007]. Therefore it appears that for large watersheds like the Lena, Yenisey, Yukon, and Mackenzie that have headwaters <60°N, a large proportion of the DOC exported is from the more northern peat rich regions. Thus the estimate of 25 Tg for watersheds draining directly into the Arctic Ocean is conservative because the undersampled high north watersheds flanking the Arctic Ocean that were extrapolated for likely have a higher DOC yield than those in Figure 4. That is, these unsampled watersheds have a greater percentage of their total area influenced by peat and probably have higher DOC yields.

3.2. Δ¹⁴C DOC

[18] Average Δ¹⁴C-DOC was 78, 14, 72, 7, and 8‰ for the Lena, Ob', Yenisey, Yukon, and Mackenzie, respectively. Flow-weighted averages were 94, 41, 106, 52, and 4‰ for the Lena, Ob', Yenisey, Yukon, and Mackenzie, respectively. The majority of samples were enriched in bomb carbon with respect to current day atmospheric values (Table 2), which is consistent with other recent measurements made on these rivers [Benner et al., 2005; Guo and Macdonald, 2006]. All rivers except the Mackenzie also showed a marked enrichment in Δ¹⁴C-DOC with increasing discharge and concentration (Figure 5). Ice-out DOC was Δ¹⁴C enriched for all watersheds excluding the Mackenzie (Table 2 and Figure 5). All rivers except the Lena have at least one depleted Δ¹⁴C-DOC value (Table 2), with the most Δ¹⁴C depleted value being winter during low flow.

[19] Smaller high-latitude watershed studies have demonstrated a Δ¹⁴C-depleted DOC pool with base flow, which is due to a higher percentage of DOC originating from water that has interacted with deeper soil profiles containing older OC [Neff et al., 2006; Schiff et al., 1998, 1997]. The relationship between flow and Δ¹⁴C (Figure 5) suggests that this is also occurring at the large watershed scale with a small, slightly Δ¹⁴C depleted component being exported with base flow. It is important to note, however, that owing to the very low flows and low DOC concentrations at base flow (Figure 2) this is a very small percentage of total export.

[20] The composite nature of bulk DOC Δ¹⁴C measurements make it difficult to ascertain ages and turnover times of DOC that is Δ¹⁴C-enriched [Raymond and Bauer, 2001c]. The dominance of young Δ¹⁴C-enriched DOC in the ice-out period which dominates the DOC flux, however, provides a unique opportunity to break DOC pools into distinctly aged fraction because it can be argued that there is no significant prebomb DOC contributing to these Δ¹⁴C values. A recent isotopic and compositional study of DOC in the Kolyma has argued that DOC exported from the ice out period originates from the leaching of surface soil layers [Neff et al., 2006] due to the frozen nature of deeper soil profiles. There is evidence for a Δ¹⁴C-depleted pool exported with base flow (Figure 5), although even these values are generally only slightly Δ¹⁴C-depleted. This slightly depleted base flow pool is presumably also present during the ice-out and the growing season periods. However, owing to the low concentrations of DOC in the low-flow base flow period (Figure 2), this pool will have little effect on the bulk Δ¹⁴C-age of the bulk DOC. Nevertheless the contribution of this pool to the bulk DOC can be removed assuming simple two end-member mixing. Removing the base flow component causes only small changes in Δ¹⁴C to the remaining DOC (Table 3), demonstrating the insignificance of the older base flow component during ice-out. Interestingly, removing the base flow component provides a fairly consistent Δ¹⁴C-DOC ratio for spring ice out DOC with an average of +100‰ (Table 3).

[21] Using the evidence that there is an insignificant amount of prebomb (<0‰) DOC in these base flow corrected Δ¹⁴C-DOC values (Table 3), provides an opportunity to constrain the relative contribution of different age classes to the bulk DOC pool. It is unlikely that riverine DOC in these rivers is represented by a single year when atmospheric Δ¹⁴C-CO₂ was last +100‰. The Lena River sample of +120‰ for instance is most likely not 100% DOC that originated from the early 1990's, when the atmosphere was last +120‰ (Figure 6). Previous studies in rivers and the coastal ocean have concluded that DOC comprises multiple fractions of varying lability [Hopkinson et al., 2002; Moran et al., 1999; Raymond and Bauer, 2000]. It has also been argued that there is a relationship between the age of material and its lability [Raymond and Bauer, 2001b].

[22] Soil studies have utilized simple diagenesis or decomposition kinetics to constrain the turnover time of soil DOC [Trumbore and Harden, 1997]. For riverine studies the approximate contribution of different age classes of DOC to the bulk pool is of interest. Information on the contribution of different age classes to DOC can also utilize decomposition kinetics if one is not constrained by the possibility of a significant old DOC component. The simplest model (Model 1, Table 4) is one that predicts the amount of DOC present from different years assuming that the DOC transported to rivers during ice-out is dominated by DOC from the previous year's primary production and that DOC contributions from subsequent years decreases on the basis of a single G decomposition constant. This type of model predicts that the mobile pool for lateral transport is constrained to the labile and semilabile pools (∼0–50 years in age) that fuel decomposition on land [Schuur and Trumbore, 2006] and does not include older pools that make up the stored peat, which we argue above is still frozen at ice-out and not evident in our samples. An assumption using this single G model is that the lateral transport of a given age class is proportional to the percentage of that age class remaining in the active layer of soils, so that

where %_t is the percentage of DOC from a proceeding year (e.g., 5 years ago), (%_t−1) is the percent remaining entering into that year, and k₁ is the decomposition constant in yr⁻¹. If k is 0.1, for instance, at year 4, 67% of the DOC is remaining (%_t−1; e^−0.1*4), leaving a 7% contribution from year 5 (0.67 − e^−0.1*5). Using this model in year 1, the percent remaining entering into that year (%_t−1) is equal to 1, or 100%.

[23] The radiocarbon data allows for the solving for the one unknown (k₁) since the inventory-weighted Δ¹⁴C content must provide the Δ¹⁴C-measured in the DOC so that

where Δ14C_atmt is the atmospheric Δ¹⁴C content at year t (Figure 6). Using this simple model and the average Δ¹⁴C of base flow corrected ice-out DOC (100‰; Table 3) predicts a k of ∼0.16 if the year 1 is 2004. Furthermore, this simple model predicts that 55% of the DOC was fixed within the last 5 years, 25% from 6–10 years ago, and ∼16% from 11–20 years ago (Table 4 and Figure 6).

[24] A more realistic model would consider multiple pools which follow equation 1, each with a different size and decay constant (i.e., a multiG approach). This however provides more than one unknown and therefore has an infinite number of answers. In order to utilize a multiG approach for three pools of DOC present in rivers (a labile, semilabile, and recalcitrant pool) two assumptions were made: (1) The proportion of the labile pool was the largest with the second and third pool making up 2/3 and 1/3 of the remaining DOC, respectively, and (2) the decay constants decrease so that k₂ is one half of k₁ and k₃ is one half of k₂.

[25] This provides a model with two unknowns, k₁ and the size of the most labile pool of DOC. It also sets a lower limit on the size of the most labile pool of 50%. Interestingly, using these simple assumptions provides a model that is fairly robust to the potential ranges of the two unknowns. At the lower limit of the size of pool 1 (50% of the total; Model 2, Table 4) in order to get a bulk Δ¹⁴C-DOC of +100‰, the multiG model predicts a k₁ of 0.29 and that 61% of the DOC was fixed within the last 5 years, 21% from 6−10 years ago, and ∼13% from 11−20 years ago (Table 4). These numbers are very similar to those obtained using the single G approach (Table 4). At an appreciably higher size of pool 1 (80% of the total; Model 3, Table 4) the multiG model predicts a k1 of 0.2 and that 57% of the DOC was fixed within the last 5 years, 23% from 6–10 years ago, and ∼14% from 11–20 years ago (Table 4).

[26] This constrained multiG model is also insensitive to how the relative size of the three pools and decomposition constants are decreased. For instance adapting model 2 so that k₂ is assumed to be 30% of k₁ and k₃ 30% of k₂, only changes the % contribution from the 11–20 age class by 3%. Altering model 2 so the second pool is assumed to make up 3/4 instead of 2/3 of the remaining DOC the % contribution from the 11–20 age class only changes by 0.1%. Finally, if there is a small pool of Δ¹⁴C depleted DOC that is not removed by removing base flow DOC, the % contributions from older age classes would be higher, making this a conservative model. It should be noted, however, that there are large differences to the first 3 years' age class using these different models.

[27] An alternative scenario would be that terrestrial ecosystems are highly efficient at recycling newly produced organic matter and that the majority of DOC that is exported is a fraction that slowly leaches from semilabile/aged pools of terrestrial organic carbon. This would provide a Gaussian-type distribution of ages that centers around an age class that is close to the historic atmospheric Δ¹⁴C-CO₂ ratio that is approximate to the Δ¹⁴C of DOC. The DOC samples, however, are associated with appreciable amounts of DON (PARTNERS unpublished data) and therefore it appears that the DOC exported during the ice-out includes labile/young fractions. Regardless, a Gaussian-type distribution would cause the largest DOC age class to be one that is approximate to Δ¹⁴C to the atmosphere. In the case of a +100‰ (Table 4) value this would mean the largest age class would be ∼10 years.

[28] It is important to note that these age class designations are approximate because they are based on composite Δ¹⁴C-DOC ratios. An argument has been made here for the application of simple diagenesis rules to DOC export during ice-out. Until a large enough data set exists and the decrease in riverine Δ¹⁴C-DOC over time can be determined and compared to the decrease in atmospheric CO₂ an exact age class distribution is not attainable. The fact that there is a significant decadal component to the DOC pool, however, is something that is necessary to obtain the enriched values measured by this study. What can be refined is the approximate amount from different age classes. The older (decadal) reworked nature of a significant fraction of the DOC (Table 4) is consistent with a study of large Russian rivers, including the Lena and Yenisey that determined that the DOM of these rivers is soil derived and has undergone a high degree of degradation [Lobbes et al., 2000]. Despite the reworked nature of a significant fraction of the DOC, we argue that the isotopic results also require a larger young age class that is labile on monthly to annual timescales and will therefore be utilized in nearshore coastal waters. This is consistent with recent studies from relatively pristine tropical rivers and salt marshes estuaries which found a young OC pool was fueling in situ decomposition [Mayorga et al., 2005; Raymond and Hopkinson, 2003], yet disparate from studies on the highly impacted Hudson which noted old DOC utilization [McCallister et al., 2004].