Sourcing and long-range transport of particulate organic matter in river bedload: Río Bermejo, Argentina

Fluvial transport of organic carbon from the terrestrial biosphere to the oceans is an important term in the global carbon cycle. Traditionally, the long-term burial flux of fluvial particulate organic carbon (POC) is estimated using river suspended sediment flux; however, organic carbon can also travel in river bedload as coarse particulate organic matter (POM_Bed). Estimates of fluvial POC export to the ocean are highly uncertain because few studies document POM_bed sources, flux, and evolution during long-range fluvial transport from uplands to ocean basins. This knowledge gap limits our ability to determine the global terrestrial organic carbon burial flux. In this study we investigate the flux, sources, and transformations of POM_Bed during fluvial transport over a ∼1300 km long reach of the Río Bermejo, Argentina, which has no tributary inputs. To constrain sourcing of POM_Bed, we analyzed the composition and stable hydrogen and carbon isotope ratios (δ²H, δ¹³C) of plant wax biomarkers from POM_Bed at six locations along the Río Bermejo and compared this to samples of suspended sediment, soil, leaf litter, and floating organic debris (POM_float) from both the lowland and headwater river system. Across all samples, we found no discernible differences in n-alkane average chain length or nC29 δ¹³C, indicting a common origin for all sampled POM_Bed. Leaf litter and POM_float nC29 δ²H values decrease with elevation, making it a useful proxy for POM_Bed source elevation. Biomarker δ²H values suggest that POM_Bed is a mix of distally derived headwater and locally recruited floodplain sources at all sampling locations. These results indicate that POM_Bed can be preserved during transport through lowland rivers for hundreds of kilometers. However, the POM_Bed flux decreases with increasing transport distance, suggesting mechanical comminution of these coarse organic particles and progressive transfer into the suspended load. Our provisional estimates suggest that the carbon flux from POM_Bed comprises less than 1 % of the suspended load POC flux in the Río Bermejo. While this represents a small portion of the river POC flux, this coarse, high-density material likely has a higher probability of deposition and burial in sedimentary basins, potentially allowing it to be more effective in long-term CO₂ drawdown relative to fine suspended particles. Because the rate and ratio of POM_Bed transport versus comminution likely vary across tectonic and climatic settings, additional research is needed to determine the importance of POM_Bed in the global carbon cycle.

1 Introduction

The burial of organic carbon (OC) in soils and sedimentary depocenters can remove carbon from the atmosphere over timescales of centuries to millennia (Hilton and West, 2020; Galy et al., 2015; Blair and Aller, 2012; Battin et al., 2009; Hayes et al., 1999; Stallard, 1998; France-Lanord and Derry, 1997; Berner, 1982). OC buried in sedimentary basins is mainly sourced from tectonically active environments, like mountainous areas where physical erosion mobilizes hillslope bedrock and soil, generating sediment (e.g., Galy et al., 2015; Blair and Aller, 2012; Stallard, 1998). Rivers play a key role as conduits in the carbon cycle, moving OC eroded from rock, soil, and vegetation from terrestrial to marine carbon reservoirs if the transported carbon is buried in a sedimentary basin (e.g., Schlünz and Schneider, 2000). The long spatial and temporal scales of this transport allow for carbon transformation during transfer and intermittent storage (Blattmann et al., 2019; Galy et al., 2008), for instance in floodplains, estuaries, and coastal mud belts (Repasch et al., 2022; Scheingross et al., 2021; Canuel and Hardison, 2016; Aller, 1998).

It is usually assumed that once particulate organic matter (POM) from terrestrial sources has been transferred into rivers, it is transported with the fine suspended sediment (Kao et al., 2014). Estimates of POM flux, as opposed to dissolved OM, suggest that rivers deliver 110–230 Mt C into the oceans annually (Galy et al., 2015). River bedload can comprise lithic fragments that contain fossil organic carbon (J. C. Smith et al., 2013; Hage et al., 2020; Kao et al., 2014), a phase of particulate OC that might be oxidized and emit CO₂ during floodplain transit (Dellinger et al., 2023). A second organic bedload component consists of plant debris.

Particulate OM includes coarse plant material and can be abundant, ranging from 10 % up to 80 % of the total fluvial OC flux (Seo et al., 2008; West et al., 2011; Kao et al., 2014; Turowski et al., 2016). Hillslope mass wasting, overland flow, and flooding of riparian zones can mobilize leaf litter and woody debris into the fluvial routing system (Turowski et al., 2016; West et al., 2011; Wohl et al., 2009). Coarse POM can float at the river surface, where it is visible and accountable (e.g., Ruiz-Villanueva et al., 2019; West et al., 2011; Wohl et al., 2009), or it can move along the riverbed where it is more difficult to observe and measure (e.g., Turowski et al., 2013). The physical, chemical, and biological breakdown of woody debris (Seo et al., 2008) resting in the landscape allows the material to become waterlogged within days or weeks (Hoover et al., 2010). This can increase the density of organic debris above the critical value of 1.0 g cm⁻³ so that it will sink to the riverbed (Turowski et al., 2016, 2013). There, it can be transported with the bedload (e.g., Schwab et al., 2022; Hage et al., 2020; Lee et al., 2019; Turowski et al., 2016; Liu et al., 2016) and has been observed to comprise up to 75 % of all coarse POM (Turowski et al., 2016).

Several studies describe fresh, coarse terrestrial organic debris transported to delta plains (Allen et al., 1979) and offshore (West et al., 2011) by turbidity currents (Hage et al., 2020; Liu et al., 2013; Tyson and Follows, 2000). When capped by muddy siliciclastic sediments, terrestrially sourced coarse POM is protected from fast degradation (Hage et al., 2020; Lee et al., 2019; McArthur et al., 2016; Sparkes et al., 2015), leading to preservation rates of up to 70 % (Hage et al., 2022; Kao et al., 2014). Coarse woody debris and litter fragments have been described in aged deep-marine fan deposits (Lee et al., 2019) and to represent up to 12 % of the mass in exhumed turbidite layers (Turowski et al., 2016; Tyson and Follows, 2000). Taken together, these observations suggest that transport of terrestrial POM in river bedload may be a relevant pathway in the global carbon cycle (Hage et al., 2020; Lee et al., 2019; Kao et al., 2014) both during source-to-sink transit and in depocenters.

However, organic debris at the riverbed is difficult to observe and quantify due to logistical difficulties in sampling and highly variable transport rates (Turowski et al., 2013). The occurrence, recruitment, sources, and fate of POM_Bed during transport have only been addressed for a few headwater streams (Turowski et al., 2016, 2013; Bunte et al., 2016; Iroumé et al., 2020; Fogel and Lininger, 2023) and even less for lowland systems (Hage et al., 2022; Schwab et al., 2022). The paucity of work on the origin of POM_Bed, and its endurance in long-range fluvial transport after erosion, makes it difficult to build a mechanistic model to predict bedload OC fluxes and quantify their role in the terrestrial OC cycle.

In this study, we evaluate the role of organic bedload in the OC cycle of a large lowland river, the Río Bermejo in northwest Argentina. This river is a major tributary to the Río Paraguay, draining a section of the eastern central Andes across an 800 km wide foreland without overwhelming human intervention (Repasch et al., 2023). The lowland portion of the Río Bermejo has no major tributaries or distributaries over a flow distance of almost 1300 km, ruling out sediment mixing complexities inherent to dendritic drainage networks. We address three questions designed to understand the role of POM_Bed in the terrestrial carbon cycle: (1) is POM transported with bedload in the lowland Río Bermejo? (2) If so, what are the source areas and mechanisms for POM_Bed recruitment? (3) Does POM_Bed survive long-range transport without transformation through the Río Bermejo? To answer these questions, we collected organic-rich material traveling at the riverbed, analyzed the geochemical composition of the collected organic-rich material, and tracked compositional changes with increasing distance downstream from the headwaters.

2 Study area, sampling methods, and analyses

2.1 Study area

The Río Bermejo, with a catchment area of 120 283 km², drains the northern–western Argentinian Andes before crossing the Gran Chaco alluvial plain and joining the Río Paraguay (Fig. 1a). The Río Bermejo headwaters extend to the eastern limit of the arid Puna Plateau at an elevation of ∼4 km. Headwater streams drain high-altitude dry grasslands, ultimately merging to form the upper Río Bermejo in the north and the Río San Francisco (a major headwater tributary of the Río Bermejo) in the south (Fig. 1). In the southern headwaters, mean annual rainfall is around 1000 mm yr⁻¹ and dense Yungas forests cover deeply weathered foothills. In the northern headwaters, the hinterland comprises the Eastern Cordillera and the sub-Andean zone with highly variable rainfall up to 1400 mm yr⁻¹ and Yungas montane forest (Fig. 1c). Steep relief and intense precipitation in these northern headwaters result in high sediment yields to the Río Bermejo.

Figure 1 Catchment overview of the Río Bermejo main stem and main headwater tributaries. (a) Topographic map (ASTGTM, NASA EOSDIS Land Processes DAAC, 2019) with sample types indicated at each sampling location. Top right: catchment location in Argentina. Bottom left: zoomed-in view of the tributary confluence at the mountain front. (b) C₄ vegetation cover as a percentage of total C₃ and C₄ vegetation cover (Powell et al., 2012). (c) Average annual precipitation rate (Hijmans et al., 2005).

At the Andean Mountain front, the upper Río Bermejo and the Río San Francisco merge to form the lowland Río Bermejo, a sand-bed river crossing the foreland basin. Just after the headwater confluence, the lowland Bermejo exhibits a braided morphology, with river width varying from 1–3 km. 175 km downstream of the confluence, the river transitions into a meandering channel, with migration rates from 40–80 m yr⁻¹. The channel width narrows to 170 m in the most downstream parts, towards the confluence with the Río Paraguay (Repasch et al., 2023, 2020; Sambrook Smith et al., 2016). River depth ranges between 4 and 10 m at high flow with channel depth increasing downstream (Sambrook Smith et al., 2016). The daily water discharge averages 432 m³ s⁻¹; however, the wet-season discharge of the austral summer accounts for ∼75 % of the annual flow, with daily discharges up to 2000 m³ s⁻¹ (Golombek et al., 2021; Sambrook Smith et al., 2016). Grain size analyses show downstream fining of the suspended river load (Repasch et al., 2020) from medium sand, on average 280 µm in the headwaters, to very fine sand to silty sand, on average 90 µm in the downstream part of the river (Sambrook Smith et al., 2016; McGlue et al., 2016). The Río Bermejo delivers ∼80 Mt yr⁻¹ of suspended sediment to the Río Paraguay. Suspended sediment input from the Andean headwaters at the Bermejo–San Francisco confluence is significantly higher at ∼103 Mt yr⁻¹, suggesting net deposition during foreland transit (Repasch et al., 2020). The Río San Francisco contributes ∼14 % to the total suspended sediment load of the Río Bermejo main stem, whereas the northern headwaters contribute ∼86 % (Repasch et al., 2020), likely reflecting the south-to-north gradient in precipitation, vegetation, and erosion rate.

Our work is focused on the Río Bermejo downstream of the eastern Andean Mountain front, from the confluence of the Río Bermejo and the Río San Francisco at 304 m a.s.l. to the Río Bermejo–Río Paraguay confluence at 50 m a.s.l. The linear distance between these points is ∼700 km, while the river has a channel length of ∼1300 km due to its high sinuosity. Over this distance, no significant tributaries join the Río Bermejo, making it an ideal setting for our source-to-sink study of organic bed material.

2.2 Sampling

During sampling campaigns in 2013, 2015, 2017, 2019, and 2020, we collected headwater and lowland floodplain leaf litter, coarse floating particulate organic matter (POM_float, >1 cm diameter), soils from dry paleo-channels, riverbank sediments, suspended sediment, and bedload material in the Río Bermejo catchment (Table 1, Fig. 1a). Sampling of bedload was curtailed by pandemic travel restrictions coming into force during fieldwork in March 2020.

Table 1 Overview of the bedload sampling sites, location, number of samples, and type of sampling.

The sampling sites at Río Colorado and Embarcación are located within the northern headwaters (HW_North), and Caimancito and Pichanal are sampling sites in the southern headwaters (HW_South) along the Río San Francisco (Fig. 1a). We sampled the Río Bermejo main stem at five separate locations between the most upstream location at Carboncito (14 km linear distance downstream of the mountain front confluence) and the most downstream site at General Mansilla (660 km linear distance downstream of the mountain front confluence) (Fig. 1, Table 1).

2.2.1 Bedload sampling: POM_Bed

In March 2020, we collected bedload material from cross-channel transects (Fig. 2a) at four locations upstream of the Bermejo–San Francisco confluence. At HW_South, we sampled the Río San Francisco at Pichanal (HW_South-1, n=13) and at Caimancito (HW_South-2, n=2). At HW_North, we sampled the upper Río Bermejo at Embarcación (HW_North-1, n=10) and the Río Colorado tributary (HW_North-2, n=4) (Fig. 1). Downstream of the confluence, we sampled the lowland Río Bermejo main stem at Puerto Lavalle (LL_-1, n=12), 481 linear kilometers (865 km streamwise) downstream of the confluence, and at El Colorado (LL_-2, n=8), 583 linear kilometers ( 1046 km streamwise) downstream of the confluence (Table 1).

Figure 2 (a) Schematic channel cross-section with transect sampling points (arrows), sample types used in this study, and indicated sample locations. (b) Cross-channel velocity profile using an acoustic Doppler current profiler at the LL_-2 station, with crosses indicating each transect sampling point at the station. Bridge pillars were at 30, 75, and 125 m, and measurements were taken ∼5 m downstream of the bridge.

We collected bedload samples with a Helley–Smith bedload sampler with a 100 µm mesh and an 8 cm×8 cm square opening (KC Denmark A/S, 2023). We ballasted the sampler with a 12 kg weight to enable sinking to the riverbed and attached a rope to both ends of the sampler. We deployed the device systematically from bridges across the Río Bermejo and Río San Francisco, where possible multiple times along a transect (Fig. 2a). We lowered the sampler to touch gently on the riverbed for 60 s and then manually pulled it up as quickly as possible, usually in less than 15 s, to minimize capture of material from shallower water column depths. We additionally collected bedload samples at one location per sampling site for a longer period, usually around 5 min, to attempt the collection of a sufficient amount of organic bedload material for compositional analysis. The wet samples were stored in airtight Whirl-Pak plastic bags and transferred to Potsdam, Germany, within 2 weeks. The samples were freeze-dried and dry-sieved over a stainless-steel sieve with 1 mm mesh size. Of the 49 bedload samples collected in 2020, 30 yielded no or less than 1 g of macroscopic organic material (POM_Bed), many from downstream sampling sites. For selected samples with abundant organic debris in both fractions, we processed and analyzed the organic bedload material following the protocol described in Sect. 3.2 for both the >1 mm and <1 mm size fractions separately. We found no significant differences between the two size fractions and therefore analyzed the remaining bedload samples as bulk.

In earlier sampling campaigns, we collected individual bedload samples throughout the catchment (Table 1): in March 2017, we collected single-point bedload samples with a self-built device at one headwater site, HW_South (Pichanal, Río San Francisco), and three main-stem Río Bermejo sites: Reserva Natural Formosa, El Colorado, and General Mansilla (Fig. 1a). In November 2019, six bedload samples were collected with a bedload grab that did not yield enough OM to permit compound-specific stable isotope measurements or significant amounts of POM_Bed. The sampling in these earlier campaigns was performed for the qualitative assessment of POM_Bed occurrence, and corresponding data were not used to quantitatively estimate POM_Bed.

For our purpose, we define POM_Bed as organic material that is entrained within the clastic bedload, transported as separate layer on top of the clastic bedload, or moves close to the riverbed. It is likely that the POM_Bed material is transported in a more extensive layer above the bed (Repasch et al., 2022; Schwab et al., 2022), also including saltating trajectories (Einstein et al., 1940; Turowski et al., 2010). The maximum particle size of the bedload samples was likely limited by the funnel opening width of 8 cm, as has been demonstrated for clastic bedload (Bunte et al., 2008), and our sample collection was restricted to the material transported within 8 cm above the bed. Our sample sizes may reflect these limitations of bedload sampling rather than a physical phenomenon. We assume, however, that these potential sampling biases do not affect the composition of the sample material.

2.2.2 Leaf litter and POM_float

During campaigns in 2015, 2019, and 2020, we collected 17 leaf litter samples within the floodplain ranging in elevation from 138 to 271 m a.s.l. and 11 leaf litter samples from HW_North and HW_South upstream of the Río Bermejo–Río San Francisco confluence between 317–854 m a.s.l. (Fig. 1). We typically collected ∼20 g of leaf litter; the samples were air-dried in the field and stored in paper bags upon arrival at the laboratory in Potsdam, Germany. There, the samples were oven-dried at 40 °C for up to 3 d and stored until further analysis.

We collected seven samples of floating, coarse particulate organic matter (POM_float, >1 mm) along the Río San Francisco (between 313 and 2458 m a.s.l.) during the high-flow season in March 2013 (Fig. 1a). These samples were collected with a net with mesh size 1 mm and held into the river water from the riverbanks until a sufficient amount of organic material was accumulated and contained organic and inorganic debris ranging from 1 mm to 10 cm. The sample material was stored in plastics bags and freeze-dried upon return to the laboratory. The dry sample masses ranged from 5 to 40 g.

2.2.3 Soil, bank, and suspended sediment

During campaigns in 2019 and 2020, we collected 15 sediment and 24 soil samples from the riverbanks and from a paleo-channel on the Río Bermejo megafan at elevations of ∼60 to 390 m a.s.l. We used a rinsed, stainless-steel shovel to collect material from 0 to 20 cm below the surface. The samples were stored in paper bags and air-dried in the field. Upon returning to the laboratory, we oven-dried all samples at 40 °C for up to 3 d and transferred them into brown glass bottles before analysis. We also include data on suspended sediment particulate organic carbon abundance and composition in the Río Bermejo from previous studies (Repasch et al., 2021, 2022). Briefly, 48 samples were collected from several locations spanning HW_South, HW_North, and the lowland Río Bermejo (Fig. 1a).

3 Methods: analysis and data treatment

3.1 Near-bed flow velocity and channel depths

During the 2020 campaign, we determined the flow velocity and channel depth at the sampling locations HW_North-1 and HW_South-1 near the confluence and downstream at LL_-1 and LL_-2. Surveys were conducted from bridges using an acoustic Doppler current profiler (ADCP; SonTek RiverSurveyor RS-M9). The ADCP was mounted on a floating board and towed by a rope from the bridges. We applied strong tension to the rope in an effort to prevent the ADCP from being submerged under the bridge. Where possible and necessary, the ADCP float was guided by a person in the water. The raw ADCP data were processed using the SonTek RiverSurveyor Live Software (Version 4.1). After quality assessment with the SonTek software, we further processed the data files using the velocity mapping toolbox (v.4.09) (Parsons et al., 2013) to calculate smoothed mean cross-sections with the river flow velocity determined at horizontal and vertical grid node spacings of 1 and 0.5 m, respectively. We extracted the river depth and flow velocity at each bedload sampling point by overlaying our GPS data points with these cross-sections. To approximate the bedload transport velocity, we multiplied the depth-averaged streamwise flow velocity from the ADCP velocity profiles by 0.7 (Chatanantavet et al., 2013).

The ADCP-estimated near-bed flow velocity ranged between −0.17–1.19 m s⁻¹ and was on average 0.5 ± 0.4 m s⁻¹ for the ensemble of sampling sites (n=4). We attribute occasional negative flow velocities to local flow patterns on large river bedforms (Allen, 1968) and assign no general significance to them (Fig. 2b). River depths measured during the wet season varied between 8.1–8.5 m at HW_South-1 and HW_North-1, 1.9–9.9 m at LL_-1, and 2.4–8.4 m at LL_-2. Due to the interrupted campaign in March 2020, these values do not cover the braided section of the river in the upper part of the eastern Andean foreland.

The high flow velocities and water depths of the Río Bermejo made clean sampling of river bedload difficult, and the composition of some samples may have been affected by admixture of suspended load during sampler recovery. This is most likely to occur in the sediment-laden upper reaches of the Río Bermejo and at downstream sites where bridge pillars cause additional turbulence.

3.2 Organic–geochemical analysis

To fingerprint source areas of the OM collected in this study, we use biomarker proxies, specifically long-chain n-alkanes, recalcitrant organic molecules that are often preserved in sediments (e.g., Thomas et al., 2021; Cranwell, 1972). Long-chain n-alkane stable carbon and hydrogen isotopes incorporate local environmental conditions at formation and can therefore help to track the source areas of POM in river sediment (e.g., Hemingway et al., 2016; Hoffmann et al., 2016; Bouchez et al., 2014; Ponton et al., 2014; Galy et al., 2011; Sachse et al., 2004).

We extracted n-alkanes from all samples and measured the compound-specific stable hydrogen and carbon isotope ratios. Soil, bank sediment, suspended sediment, and POM_Bed samples were ground with a mortar and pestle and lipid compounds were extracted with 9:1 methanol : dichloromethane using a Thermo Fisher Dionex Accelerated Solvent Extraction system (ASE 350). The non-polar n-alkane fraction was separated from the total lipid extract over silica gel columns with glass-fiber filters at the base and top (pore size 60 Å, 230–400 mesh particle size) by automated solid-phase extraction (SPE) with a Gilson ASPEC GX-271 and n-hexane as a solvent, following procedures described by Rach et al. (2020). n-Alkanes were extracted manually from leaf litter and POM_float by immersing between 0.2 and 10.0 g of the dried samples in a 9:1 methanol : dichloromethane mixture and placing them into an ultrasonic bath at 40 °C for 20 to 40 min. The n-alkanes were separated from the total liquid extract by solid-phase extraction (SPE) over a silica-gel-equipped 6 mL glass column (Macherey-Nagel, Düren, Germany) using n-hexane as a solvent (Rach et al., 2020). The solid parts stayed behind and samples were rinsed properly with 9:1 DCM$/$MeOH. The lipid-containing solvents were dried in the TurboVap and the lipids were transferred to a smaller vial in approximately 1.5 mL 9:1 DCM$/$MeOH. The solvents of the total liquid extract (TLE) were evaporated again with N₂ gas and the samples were dissolved in 1.5 mL n-hexane. To quantify n-alkane concentrations per sample, we added an internal standard, 5α-androstane (10 µg), and measured the samples in an Agilent gas chromatograph (GC 7890-A) with a flame ionization detector (FID) and a coupled single quadrupole mass spectrometer (MS 5975-C).

We quantified the abundances of n-alkane homologues relative to the internal standard using the FID chromatograms. We calculated the average chain length (ACL₂₅₋₃₃) of the most abundant n-alkanes with chain lengths between 25 and 33 as

$$\begin{array}{}\text{(1)}& {\text{ACL}}_{\mathrm{25}-\mathrm{33}}={\displaystyle \frac{\sum (C_n\times n)}{\mathrm{\Sigma }{C}_n}}.\end{array}$$

n is the number of carbon atoms of each n-alkane, and C_n is the concentration of each n-alkane with n carbon atoms. The subscripts in ACL₂₅₋₃₃ refer to the chain length range analyzed.

We determined the carbon preference index (CPI₂₅₋₃₃) after Bray and Evans (1961) for n-alkanes with 25 to 33 carbon atoms using

$$\begin{array}{}\text{(2)}& {\text{CPI}}_{\mathrm{25}-\mathrm{33}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\times \left({\displaystyle \frac{\sum C_{\mathrm{25}-\mathrm{33}\text{odd}}}{\sum C_{\mathrm{24}-\mathrm{32}\text{even}}}}+{\displaystyle \frac{\sum C_{\mathrm{25}-\mathrm{33}\text{odd}}}{\sum C_{\mathrm{26}-\mathrm{34}\text{even}}}}\right).\end{array}$$

$\sum C_{\mathrm{25}-\mathrm{33}\text{odd}}$ is the sum of the concentration of odd-chained n-alkanes with chain lengths between 25–33, $\sum C_{\mathrm{24}-\mathrm{32}\text{even}}$ is the sum of the concentration of even-chained n-alkanes with chain lengths between 24–34, and so forth.

To measure compound-specific hydrogen and carbon isotope ratios of the n-alkanes (expressed as δ²H, δ¹³C values), we used a Trace GC 1310 (Thermo Fisher Scientific) connected to a Delta V Plus isotope ratio mass spectrometer (IRMS) (Thermo Fisher Scientific), following the procedures described by Rach et al. (2020). n-Alkane δ²H and δ¹³C values were measured in duplicates. For each sample run, we measured the n-alkane standard mix A6 (with n-alkane chain lengths ranging from nC₁₆–nC₃₀) with known δ²H values obtained from Arndt Schimmelmann (Indiana University) for correction and transfer to the VSMOW scale. The H3+ factor from the ²H measurements was 5.9 ± 0.8 mV. The standard deviation for δ²H of the standard measurements was 2.3 ‰ ± 0.7 ‰, and the mean standard deviation for the samples was 1.1 ‰ ± 1.6 ‰. The mean standard deviation for ¹³C measurements of the standard measurements was 0.11 ‰ ± 0.06 ‰, and the mean standard deviation for the samples was 0.15 ‰ ± 0.23 ‰.

Figure 3 Left panels (a–f): average and standard deviation (bars and whiskers) of the long-chain n-alkane distribution (nC₂₅–nC₃₅) as the concentration per sediment dry weight (µg g⁻¹) of (a) POM_float, (b) leaf litter, (c) POM_Bed, (d) floodplain soil, (e) bank sediment, and (f) suspended sediment. Note the different y-scale ranges. Right panels (g, h): nC₂₉ versus nC₃₁ n-alkane (g) δ²H and (h) δ¹³C values of all OM sample types, as well as the linear regression equation and R² values.

Organic–geochemical raw results

We calculated the ACL₂₅₋₃₃ index for 201 samples collected in the Río Bermejo catchment, including POM_Bed, leaf litter, POM_float, soil, sediment deposits, and river suspended sediment. n-Alkanes with chain lengths ranging from nC₂₅–nC₃₅ had the highest abundance across all samples (Fig. 3a–f). The main components were nC₃₁, nC₂₉, and nC₂₇, followed by nC₃₃ and nC₂₅, while even-carbon-numbered n-alkanes were minor components in all samples. Leaf litter, POM_Bed, and soil samples had similar concentrations of nC₂₉ and nC₃₁. In riverbank and suspended sediment samples, nC₃₁ was the dominant n-alkane, whereas nC₂₇ was the main component of the POM_float samples. We measured the δ²H and δ¹³C values for the dominant chain lengths, nC₂₉ and nC₃₁. Because these n-alkanes showed a significant correlation for δ¹³C (R²=0.95, p<0.01) and δ²H values (R²=0.81, p<0.01, Fig. 3g and h), we will focus on the nC₂₉ values.

For a subset of 10 POM_Bed samples with abundant organic debris in both the >1 mm and <1 mm fractions (HW_South-1, HW_South-2, HW_North-1, HW_North-2, LL_-1, LL_-2), we analyzed the n-alkane composition of POM_Bed in both size fractions. In most of these samples, the two size fractions had distinct ACL₂₅₋₃₃, CPI₂₅₋₃₃, nC₂₉ δ²H, and nC₂₉ δ¹³C values (Fig. S1 in the Supplement). However, no discernible pattern emerged in these differences. For this reason, and because size-specific biomarker analyses were not possible for smaller sample amounts, we decided to work with n-alkane data from the bulk sample material. For samples analyzed as two separate fractions, we determined values for the bulk sample material as the weighted average of the two fractions. These values are dominated by the fraction <1 mm, which had the greater mass in all cases. In the following, we will refer to the results of analyses of the total particulate organic carbon of the sampled bedload material in the full size range as POM_Bed.

3.3 Data analysis

All data analyses were performed using R 4.1.2 GUI 1.77 High Sierra build (8007). We report the range of the values as average ± standard deviation and the number of samples. We used the Kolmogorov–Smirnov test to account for the non-normal distribution of our data and the Mann–Whitney U test to test for significant differences between independent sample groups. We used linear regression and associated R² values to test for significant trends in our data. Significance levels are reported at the 95 % confidence interval (p value<0.05).

4 POM_Bed presence, source, and ability to survive long-distance transport

We aim to answer questions regarding the source, transport, and survival of POM_Bed in the Río Bermejo. First, we want to observe when and where POM_Bed is in active fluvial transport and how POM_Bed composition varies throughout the catchment. We then investigate the sources of POM_Bed by determining the elevation and climate of the OM source using the n-alkane δ²H value, assess the biological variability of POM_Bed sources using δ¹³C values combined with ACL₂₅₋₃₃, and estimate POM_Bed maturity using CPI₂₅₋₃₃. Finally, we apply a mixing model using these geochemical data to define distinct POM_Bed sources. Collectively, this approach allows us to evaluate the survival and fate of POM_Bed during long-range fluvial transport in the Río Bermejo.

4.1 Is there transport of POM_Bed in the lowland Río Bermejo, and what is its geochemical composition?

4.1.1 Bedload mass and material composition

The mass of bedload samples comprising organic and clastic particles was up to an order of magnitude greater at the headwater sites than at the downstream sites for the same sampling procedure. Sampling at HW_South-1 and HW_North-1 yielded totals of 560 g min⁻¹ (n=13) and 259 g min⁻¹ (n=10), respectively, while at LL_-1 we collected a total of 51 g min⁻¹ (n=12) and at LL_-2 77.1 g min⁻¹ (n=8), respectively (Fig. 4d). The mass of dry POM_Bed in individual bedload samples varied from 0 to ∼20 g, without clear spatial patterns. The mass of organic bedload scaled loosely with the amount of clastic sediment collected (Fig. 4d), but there was no correlation of the collected material and near-bed velocities. The material collected in POM_Bed samples ranged from fragile, intact leaves and twigs to robust wood fragments with frayed edges and organo-clastic aggregates in the size range from <1 mm to ∼10 cm, often mixed with fine organic debris <1 mm (Fig. 4). The organo-clastic aggregates were easily dissociated during sieving, separating into fine sand and silt, and apparently degraded OM particles ranging from <1 mm to >1 cm (Fig. 4).

Figure 4 Examples of captured POM_Bed as a (a) bulk fraction from the Río San Francisco at HW_South-1 and (b) Río Bermejo at LL_-1 in particle size separates: >1 mm aggregated (left) and dissociated (middle) as well as <1 mm mixed with clastic material. (c) Bulk at the Río Bermejo at HW_North-1 and (d) sampled bulk bed material (empty symbols) and POM_Bed>1 mm (filled symbols) per sampling point at all sampling locations. Note that POM_Bed>1 mm is shown in g min⁻¹ and bulk bed material in $\mathrm{g}{\min }^{-\mathrm{1}}\times {\mathrm{10}}^{-\mathrm{1}}$ to account for the mass difference of the samples.

Chemical maturity of an OM sample can be expressed with the CPI₂₅₋₃₃, where lower values indicate a higher sample maturity due to chemical alterations of the OM. The strong predominance of odd-over-even n-alkane chain lengths in POM_Bed (CPI₂₅₋₃₃ average: 7.4 ± 3.0, n=39, Fig. 5) and the ACL₂₅₋₃₃ averaging 29.6 ± 0.9 (range: 27.4–31.6, n=39, Fig. 5) throughout the catchment indicate relatively fresh or little-degraded vascular plant material in the Río Bermejo bedload (Eglinton and Hamilton, 1967; Bray and Evans, 1961). The occurrence of low CPI${}_{\mathrm{25}-\mathrm{33}}<\mathrm{1}$ values in our dataset may be due to the occasional presence of burned OM, such as charcoal particles, in the river bedload. While it is difficult to evaluate the volumetric contribution of this low-CPI₂₅₋₃₃ material, it is negligible relative to other POM_Bed sources.

Figure 5 Summary of (a) ACL₂₅₋₃₃, (b) CPI₂₅₋₃₃, (c) nC₂₉ δ¹³C, and (d) nC₂₉ δ²H values of POM_float from headwater and floodplain leaf litter. POM_Bed from the northern headwater, southern headwater, downstream floodplain, floodplain soil, bank sediment, and suspended sediment. The box plot width shows the interquartile range, with the black line showing the median and the whiskers the minimum and maximum range of the data without outliers. Black dots indicate outliers with a 0.75 quantile + 1.5 × interquartile range and a 0.25 quantile − 1.5 × interquartile range, respectively.

The POM_Bed fraction <1 mm consisted of mixed organic debris and clastic sediment, and the coarser fraction >1 mm was composed predominantly of organic debris. Only at location HW_South-1 did the Río San Francisco carry substantial quantities of pebbles in the fraction >1 mm, while organo-clastic aggregates were more abundant at downstream sites. The geochemical analysis also revealed substantial differences in the size fractions >1 mm and <1 mm POM_Bed in nearly all samples with sufficient mass for an articulated analysis (Fig. S1). The compositional differences between size fractions, as measured in our limited sample set, lack any discernible systematics, suggesting that variable sourcing affects both size fractions. Variable proportions of woody debris versus leaf fragments, biogeochemical quality, and size distribution indicate significant variability in the sourcing of materials of different sizes at each sampling location of POM_Bed even within a single site (Fig. S1b and f). This implies that OM fragments in the Río Bermejo bedload have diverse histories and indicates incomplete mixing both across the channel and along the length of the river, possibly due to highly active channel migration in some parts of the foreland. We did not measure the particle size distributions within individual samples to have sufficient sampling material for geochemical analysis, but we detected a clear reduction in grain size between headwater and downstream locations (Fig. 4).

Geochemically, the similarity of POM_float, leaf litter, soils, sediments, and POM_Bed suggests that all of the former can be sources of the sampled POM_Bed and no major transformations have occurred during fluvial transit (Fig. 5). Leaf litter and soil inputs from hillslope and riparian areas have potential to be important sources adding to the POM_Bed from the headwaters, yet the sampled soils contained little to no coarse POM, suggesting that eroded soil is more likely to become river suspended sediment than POM_Bed. Therefore, POM_float and leaf litter are likely more important contributors to the POM_Bed load. Another presumable POM_Bed source is abundant woody debris in the catchment; however, this material contains insufficient n-alkanes to perform the same analysis. POM_Bed CPI₂₅₋₃₃ values (average: 7.4 ± 3.0, n=39) were not significantly different from leaf litter (7.4 ± 4.0, range:1.0–19.8, n=28) and riverbank sediments (6.5 ± 3.7, range: 0.3–13.7, n=18). However, on average POM_Bed CPI₂₅₋₃₃ values were significantly higher than in soils (5.9 ± 3.6, range: 0.2–16.3, n=29) and suspended sediments (5.5 ± 1.0, range: 1.1–7.8; n=41), indicating a lower maturity of POM_Bed. CPI₂₅₋₃₃ values of POM_float were substantially higher than values in all sediment samples (15.0 ± 8, range: 2.4–27.9, n=6), indicating it was, on average, less degraded than other sampled OM.

Before we analyze the details of POM_Bed sourcing, we discuss the mechanisms of POM_Bed recruitment.

4.1.2 Mechanisms of recruitment and transport of POM_Bed

POM_Bed is ubiquitous in the bedload of the Río Bermejo throughout the catchment during the high-flow season (December–April), when the South American monsoon drives intense precipitation events throughout the study area. In contrast, we did not recover significant amounts of POM_Bed during the low-flow season. The high variability in δ²H values from POM_Bed sampled in 2017 and 2020 and the lack of significant amounts of POM_Bed during the dry season of 2019 suggest that the processes of recruitment and transport occur quickly, possibly on (sub-)seasonal timescales. The significantly higher CPI₂₅₋₃₃ values of POM_float (Fig. 5b) suggest that the degradation of organic debris occurs after it becomes waterlogged. However, aerobic decomposition seems unlikely during active bedload transport due to the high turbidity and depth of the river water. Instead, fresh plant debris is likely stored intermittently on hillslopes, on floodplains, or in-stream during the low-flow season (May–November), where it can be degraded to various extents and waterlogged. With the onset of strong precipitation and high water levels, it is mobilized by overland flow, hillslope mass wasting, or lateral channel erosion (e.g., Wohl et al., 2019; Turowski et al., 2016; J. C. Smith et al., 2013; Hilton et al., 2012, 2008; Selva et al., 2007) and subsequently transported as POM_Bed (Turowski et al., 2016) during the high-flow season.

In the headwaters, the erosive potential during the high-flow season may exceed the production of organic debris in the riparian corridor on a seasonal timescale, limiting the supply of organic debris to the channel (Hilton et al., 2012; Yager et al., 2012; Garcia et al., 1999) and, consequently, the recruitment of POM_Bed. In this case, POM_Bed export may peak early in the high-flow season and diminish over the course of the season as the supply of organic debris is progressively reduced with each subsequent rain storm and POM_Bed travels steadily downstream, similar to findings of seasonal sourcing of suspended OC at the Río Bermejo (Golombek et al., 2021).

In the lowland floodplain, POM_Bed recruitment is also highly seasonal, but the mechanism of recruitment differs from the steep headwater portion of the catchment (Wohl et al., 2019). Lateral channel migration eats into the floodplain forests, mobilizing large volumes of sediment, soil, leaf litter, and standing biomass with each bank failure. Bank erosion is most active at the peak of high-flow season, providing a source of fresh OM to POM_Bed and contributing to the poorly mixed nature of POM_Bed. However, poorly mixed samples could be expected at high rates of lateral migration from the upstream confluence to ca. 400 km downstream, but enhanced mixing seems more likely in segments with less active migration in the furthest downstream reach of the Río Bermejo. Our results do not indicate enhanced mixing downstream, possibly because of the quick transport timescales compared to the production and recruitment of POM_Bed.

Both upland and lowland lateral erosion processes are markedly reduced during the low-flow season, aligning with our observations of ubiquitous POM_Bed transport in December–April and negligible POM_Bed transport in May–November.

4.2 What are the source areas for recruitment of POM_Bed, and does it survive long-range transport?

In the previous section, we concluded that POM_Bed is a heterogenous mixture of OM from various sources in the catchment. In this section, we first aim to determine the stable carbon and hydrogen isotope composition of potential source materials. Subsequently, we apply a model to define the mixing space using source and POM_Bed geochemical compositions. In this, our aim is to determine the sources of Río Bermejo POM_Bed in order to determine its transformation and fate during long-range transit.

4.2.1 Biomarker stable isotope insights into POM_Bed source areas

The stable carbon and hydrogen isotope composition of long-chain n-alkanes is used in paleoclimate research to reconstruct continental paleo-vegetation (e.g., Sachse et al., 2012; Huang et al., 2007; Schefuss et al., 2005; Freeman and Colarusso, 2001) or to deduce environmental conditions and potential sources of OM (e.g., Hemingway et al., 2016; Bouchez et al., 2014; Galy et al., 2011; Sachse et al., 2004). Building on these approaches, we use ACL₂₅₋₃₃, nC₂₉ δ¹³C, and δ²H isotope compositions (Fig. 5) to determine the source areas of POM_Bed in the Río Bermejo catchment.

n-Alkane δ¹³C values in plants are controlled by the plant metabolism type, with C₃ plants around −37 ‰– −27 ‰ and C₄ plants around −22 ‰– −10 ‰ (e.g., Garcin et al., 2014; Huang et al., 2007; Freeman and Colarusso, 2001; Collister et al., 1994). Río Bermejo POM_Bed nC₂₉ δ¹³C values averaged −32.4 ± 2.9 (range: −36.9 ‰– −24.4 ‰, n=35, Fig. 5d), reflecting an overall C₃-dominated plant input (Garcin et al., 2014). Silva et al. (2011) identified a shift towards C₃ vegetation in our study area in the last ∼3000 years. This imposes an age cap for the organic bedload material, below which our observation of the predominance of fresh to degraded organic material can be accommodated.

There were no distinct differences between the POM_Bed δ¹³C values and those for catchment soil (−34.5 ‰ ± 2.2 ‰, range: −40.2 ‰– −29.1 ‰, n=22), deposited sediment (−33.7 ‰ ± 4.5 ‰, range: −38.0 ‰– −18.9 ‰, n=16), and floodplain leaf litter (−33.9 ‰ ± 2.2 ‰, range: −39.1 ‰– −30.6 ‰, n=15), signaling that all these reservoirs are dominated by C₃ inputs (Fig. 5a and c). However, local findings of more positive δ¹³C values and higher ACL₂₅₋₃₃ (Fig. 7) suggest C₄ plant contributions to POM_Bed, possibly from maize cultivation in the HW_South catchment area (Powell et al., 2012) and to a lesser extent in the lowland reach near LL_-2 (Fig. 1b). Because lateral channel migration rates are low in the downstream reaches of the Bermejo (Repasch et al., 2023, 2020) it is unlikely that local C₄ plants will make it into the bedload in any significant quantities. Finding signatures of C₄ plants in downstream lowland bedload samples provide evidence for long-range transport of POM_Bed from the Bermejo headwaters to nearly 1000 km downstream.

nC₂₉ δ²H values can serve as a proxy for the local meteoric water composition taken up by plants, as it is influenced by the water incorporated into plants during photosynthesis (Sachse et al., 2012; Hou et al., 2008; Chikaraishi et al., 2004) and is strongly related to temperature, humidity, rainfall amount, moisture source, and elevation (e.g., Walker and Richardson, 1991; Allison et al., 1984; Stewart and Taylor, 1981). In our study area, increasing altitude causes decreasing δ²H values in meteoric water captured in plant tissue: more negative δ²H values corresponded to higher elevations. Nieto-Moreno et al. (2016) measured soil nC₂₉ δ²H values ranging from −150 ‰ to −110 ‰ in samples collected along a valley transect ranging from ∼300 to ∼4000 m in elevation at a headwater tributary of the Río Bermejo between 22 and 24° S, and the same pattern was described for stream water δ²H values in this region (Rohrmann et al., 2014). Our samples follow these systematics, with δ²H values averaging −136 ‰ ± 15 ‰ (range: −155 ‰– −112 ‰, n=13) for leaf litter sampled at ∼270–320 m elevation in the Chaco lowland and more negative δ²H values in both leaf litter (−168 ‰ ± 21 ‰, range: −187 ‰– −130 ‰, n=9) and POM_float (−160 ‰ ± 7 ‰, range: −171 ‰– −153 ‰, n=5) collected upstream of the confluence.

²H depletion correlates significantly with increasing elevation in floodplain and headwater leaf litter samples (Fig. 6, $y=-\mathrm{640}-\mathrm{6.4}x$, R²=0.557, p<0.01). This is true also for the POM_float samples but not significantly ($y=-\mathrm{3606}-\mathrm{30}x$, R²=0.14, p>0.05). The sampling transect for POM_float covers a rapid westward precipitation decline, where samples from >1000 m a.s.l. are likely sourced from arid areas, whereas the samples from lower elevations receive orographic precipitation, suggesting that the gradient in precipitation amount may cause the observed ²H depletion and the worse fit in POM_float samples.

Figure 6 Sampling elevation and nC₂₉ δ²H values of floodplain and headwater leaf litter, as well as POM_float. Colors represent the organic matter type; symbols represent the catchment source area. The linear regression was done using floodplain and headwater leaf litter data. Linear regression of POM_float samples is not shown and was not significant. Measurement uncertainty is plotted inside the symbol.

Nevertheless, the elevation-dependent δ²H composition is caused by a precipitation-dependent change in δ²H values of the areas covered by our study, as shown by Nieto-Moreno et al. (2016) and Rohrmann et al. (2014): upland areas are depleted in ²H compared to the Río Bermejo lowland, where convective precipitation is the main source of moisture (Rohrmann et al., 2014). We take advantage of this elevation-dependent trend in source POM δ²H values to identify source elevations of the POM_Bed samples, where more negative δ²H values indicate a higher-elevation origin, and less negative δ²H values signal a low-elevation floodplain origin. POM_Bed samples do not show such a correlation of sampling elevation with δ²H values, and we suggest this is due to transport between production and sampling of this material.

4.2.2 Mixing model analysis

We defined three potential POM_Bed sources from coarse organic debris we sampled at distinct elevations in the catchment: floodplain leaf litter (<320 m), headwater leaf litter (320–1000 m), and headwater POM_float (>320 m). Significant end-member unmixing of the respective sources to the mixed POM_Bed signal is not expedient using the geochemical proxies applied in this study. Instead, we aim to understand the mixing range of the widely spread POM_Bed. We determine the range of a possible POM_Bed mixing signal of the sources within the geochemical parameters and, in addition, determine potential missing POM_Bed sources. We use a mixing space model developed by J. A. Smith et al. (2013). In short, the model uses Monte Carlo simulations to iterate mixing spaces (“convex hulls”), outlining the probability that observed POM_Bed sample compositions can be explained by a mixing model of the proposed sources. The model utilizes resampled source averages and their standard deviations and considers the distribution of the mix data (POM_Bed data) within a pre-defined mixing space.

For our purpose, we assume uniform source mixing of the POM_Bed samples without fractionation from the source composition of POM_float, floodplain, and headwater leaf litter $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/{\text{ACL}}_{\mathrm{25}-\mathrm{33}}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ to POM_Bed. We use minimum and maximum boundary conditions based on our source data to define the extent of the initial mixing space: 25 to 35 for the ACL₂₅₋₃₃, −190 ‰ and −110 ‰ for δ²H values, and −40 ‰ and −20 ‰ for δ¹³C values. The initial mixing space is iteratively adapted over 2000 iterations using source data average and standard deviation, resampled from a normal distribution. Through each iteration, a point-in-polygon algorithm tests if the mixed POM_Bed data remain within the iteratively adapted mixing space. This ensures an “optimization” of the mixing space according to our POM_Bed data, permitting us to evaluate if the actual POM_Bed sources are not fully represented by our included sources. The point-in-polygon algorithm is applied to a testing grid within the mixing space. With a grid resolution of 500 ‰², the point-in-polygon algorithm is tested on 500×500 values per iteration within the mixing region. Simultaneously, the area of the mixing space is assessed, and the variance between all previous iterations is calculated. The stabilized variance value, then, represents the optimized mixing space area for the available source data and the mixed POM_Bed data.

The variance of the convex hull area stabilized after ∼1000 iterations for the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/{\text{ACL}}_{\mathrm{25}-\mathrm{33}}$ model at a variance of 40 ‰² and after ∼1000 iterations for $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ at a variance of 60 ‰². The resulting mixing regions were not sensitive to variations of the boundary conditions. The results are plotted as derived mixing regions, with different levels of confidence representing the likelihood with which the observed data can result from mixing of the source data (Fig. 7).

Figure 7 Upper panels (a, b): headwater and floodplain leaf litter and POM_float samples with the average ± standard deviation depicted as potential POM_Bed sources using (a) nC₂₉ δ²H versus nC₂₉ δ¹³C and (b) nC₂₉ δ²H versus ACL₂₅₋₃₃. Lower panels (c, d): grayed-out area symbols and colored averages correspond to the OM sources from panels (a) and (b). Colored symbols are POM_Bed samples from headwaters and floodplains for (a) nC₂₉ δ²H versus nC₂₉ δ¹³C and (b) nC₂₉ δ²H versus ACL₂₅₋₃₃. Colored lines are probability contours of the simulated mixing area (J. A. Smith et al., 2013) using the organic matter source. The outermost contour represents the 5 % confidence level and the innermost contour the 95 % confidence level. Labeled POM_Bed samples are those that fall outside the mixing area and the source area of both plots. Uncertainty is plotted inside the symbol. Linear regression in panels (a) and (c) was conducted using lowland leaf litter data.

4.2.3 Mixing model insights into POM_Bed source areas

We proceed to use floodplain leaf litter, headwater leaf litter, and headwater POM_float ACL${}_{\mathrm{25}-\mathrm{33}}/n{\mathrm{C}}_{\mathrm{29}}$ δ²H and nC₂₉ $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}/n{\mathrm{C}}_{\mathrm{29}}$ δ²H values as potential sources to create mixing regions for our POM_Bed samples (Fig. 7a and b). Half of the samples follow the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ compositional trend of leaf litter ($y=-\mathrm{234.8}x+\mathrm{10.8}$, R²=0.94), with nC₂₉ δ²H values ranging from −174 ‰– −112 ‰ and δ¹³C values ranging from −39.0 ‰– −30.0 ‰ (Fig. 7c). This indicates that locally supplied floodplain leaf litter is an important source to the sampled POM_Bed, as already suggested earlier. In these samples, this lowland-derived plant material dominates over any existing bedload sourced from upstream areas.

A total of 20 % of the POM_Bed samples fall within the region of the nC₂₉ $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}/n{\mathrm{C}}_{\mathrm{29}}$ δ²H mixing space defined by headwater sources. This suggests input from at least one high-elevation upland source. Evidence for this high-elevation OM source in lowland POM_Bed samples confirms that mountain-derived OM can survive long-range bedload transport in the Río Bermejo.

One-third of the POM_Bed samples fall outside the 95 % confidence interval of the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ mixing space, while some POM_Bed samples collected from both the headwaters and the lowland are not within the wider mixing space constrained by $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}/{\text{ACL}}_{\mathrm{25}-\mathrm{33}}$ (Fig. 7c and d). These samples have nC₂₉ δ²H values ranging from −125 ‰– −150 ‰ and nC₂₉ δ¹³C values ranging from −27 ‰– −33 ‰ (see labeled POM_Bed samples in Fig. 7). This suggests input from at least one additional, herein unconstrained, source, and the comparably heavy δ¹³C values suggest this could contain C₄ plant contributions. C₄ plants grow only sparsely in the catchment (Fig. 1b) and mostly on agricultural land in HW_South. We suggest that the missing source in our samples is farmland OM, which would indicate that the POM_Bed carbon flux can be directly influenced by anthropogenic land use.

Sediment load data from gauging stations on the upper Río Bermejo and Río San Francisco indicate that HW_North contributes 6 times more suspended sediment load to the lowland Río Bermejo than HW_South (Repasch et al., 2020). There is a north-to-south rainfall gradient with almost 3 times more precipitation in HW_North (Fig. 1c, Hijmans et al., 2005), indicating higher erosion potential and possibly causing more recruitment of organic bedload from this area (Galy et al., 2015). This suggests that the HW_North source area should dominate organic input from the headwaters. However, this cannot be constrained with the geochemical proxies adopted in this study.

4.3 Long-range transport of POM_Bed

Recent studies show that fresh, coarse organic debris can generally be found near the riverbed (Repasch et al., 2022; Schwab et al., 2022; Feng et al., 2016), but the implications of this transport for the OC cycle are underconstrained. Using our new understanding of the POM_Bed source areas and the evolution downstream, we now consider the fate of POM_Bed during long-range fluvial transport through the Río Bermejo lowland river system.

The proportion of POM transported at the bed versus in suspension is a function of particle size, density, and shape, the recalcitrance of the POM, the river's turbulence, sediment load, the flow close to the bed (Turowski et al., 2016; Nichols et al., 2000), and secondary flow motions (Schwab et al., 2022). We assume that POM_Bed moves as clastic bedload, with a pace of about 0.7 times the depth-averaged flow velocity (Chatanantavet et al., 2013), which is around 0.46 m s⁻¹ (average of positive flow velocities, Dosch et al., 2023) at both the Bermejo–San Francisco confluence and downstream locations. At this velocity, a POM_Bed parcel could be transported 1300 km through the Río Bermejo floodplain in ∼45 d. Moreover, slightly higher CPI₂₅₋₃₃ values at downstream sites (average: 8.2 ± 3.0, range 4.9–15.2, n=14, Fig. 5) versus upstream (average: 6.9 ± 3.0, range 0.1–11.7, n=25, Fig. 5) suggest there is no systematic, progressive degradation of POM_Bed during long-range fluvial transport, likely because transport timescales are too short to produce significant chemical degradation. This demonstrates that POM_Bed can travel from the Andean headwaters into the Río Paraguay within one single high-flow season in the Río Bermejo. It explains the absence of POM_Bed during the low-flow season when little OM is likely to come in from outside the river channel. Moreover, it elucidates the occurrence of coarse organic debris in turbidity currents and turbidite deposits, where it is buried under finer clastic sediments (Hage et al., 2020; Lee et al., 2019; McArthur et al., 2016; Sparkes et al., 2015).

Once entrained, POM_Bed likely moves as a semi-separate layer in deep and fast-flowing parts of the channel cross-section where bed shear stress is greatest. Without substantial mixing, parcels of POM_Bed may shuffle downstream due to discrete erosion events at high flow and deposition at low flow. Episodic flushing events of the channel can transport parcels efficiently (Heijnen et al., 2022) and facilitate waterlogging of the organic debris (West et al., 2011). Waterlogged coarse organic debris is not prone to spilling onto the riverbank or across levees into adjacent flood basins and can progress downstream. Thus, the downstream advection of POM_Bed in the Río Bermejo may be sustained by seasonal high flow driven by monsoonal rainfall and quasi-uninterrupted by intermittent deposition. The helical flow induced by the constant meandering of the lower Río Bermejo could augment this process (Schwab et al., 2022). As such, significant amounts of OM could be efficiently transported with bedload through large drainage basins from mountainous uplands to marine basins on sub-seasonal timescales.

Mechanical comminution of POM_Bed during long-range transport

Despite evidence for survival during long-range fluvial transport, we might expect that the physical interaction between bedload OM and clastic particles can cause comminution of organic particles (Dosch et al., 2021; Scheingross et al., 2019; Turowski et al., 2016; Hilton et al., 2012; Attal and Lavé, 2009; Nichols et al., 2000) and that comminuted POM_Bed transfers into the river suspended load. The observed decrease in POM_Bed sample size with increasing distance downstream from the confluence suggests a progressive loss of POM_Bed (Fig. 4d), despite the input of OM from lateral erosion of the lowland floodplain. The similarity between n-alkane δ¹³C, δ²H, and ACL₂₅₋₃₃ values of the Río Bermejo suspended sediment and POM_Bed (Fig. 5) implies that the two fractions share an origin. River suspended sediment samples yielded similar CPI₂₅₋₃₃ values, on average 5.5 ± 1.0 (range: 1.1–7.8; n=41), but with less variability and values only as high as 7.8, suggesting advanced mixing and maturity compared to POM_Bed (average: 29.6 ± 0.9, range 27.4–31.6, n=39).

The particle size distributions of POM_Bed samples from downstream locations were considerably finer than samples from the confluence, even if the recovered sample sizes did not permit separation of aliquots for grain size measurements. As flow velocities increase, POM_Bed particles may increasingly transfer into suspension rather than traveling at the riverbed (Turowski et al., 2016). This transition is smooth, and POM_Bed may move in saltation at the transition between the two transport modes (e.g., Turowski et al., 2016; Nichols et al., 2000). However, flow velocities along the Río Bermejo main stem show no significant variability; therefore, transfer from POM_Bed to the suspended load is more likely to result from the comminution of coarse particles rather than changes in river hydrodynamics. To our knowledge, there are no experimental or field data showing comminution of particulate organic matter. However, Merten et al. (2013) suggested that physical breakage of large woody debris in streams is likely dominantly controlled by the structural properties and the in-stream position of the organic matter as opposed to hydraulic and geomorphic variables, concurring with our interpretation. Further research is needed to determine the scope and controlling factors of physical decay of coarse organic matter in the water column.

5 Synthesis: bedload carbon fluxes at the Río Bermejo

We close with a provisional examination of the organic carbon flux associated with bedload transport, its relation to suspended transport, and its role in the terrestrial carbon cycle. We extrapolated our local, high-flow season point measurements of POM_Bed across the respective river transects, ignoring any OM particles <1 mm in near-bed transport, and used a simple upscaling approach to estimate the flux of organic carbon with bedload POC_Bed (in t C yr⁻¹).

$$\begin{array}{}\text{(3)}& {\mathrm{\text{POC}}}_{\text{Bed}}={\displaystyle \frac{\sum {\text{POM}}_{\text{Bed}}\times \frac{\text{Transect width}\times \mathrm{0.5}}{\text{Funnel width}}\mathrm{0.58}}{t_{\text{sampling}}}}\times t_{\text{transport}}\end{array}$$

We sum the sampled POM_Bed (POM_Bed>1 mm, in mass × min⁻¹), apply the van Bemmelen factor of 0.58 to represent the carbon content in (soil) organic matter (Allison, 1965), and divide by the sampling time t_sampling (in min) along each transect using the respective standard deviation of POM_Bed to define an upper and lower boundary of the estimated bedload carbon flux. To do this, we extrapolated the samples obtained with a funnel opening of 0.08 m to the central 50 % of the river channel (transect width × 0.5, in m²), where it is assumed that POM_Bed transport is concentrated. Since we did not capture significant amounts of POM_Bed during the dry season, we assumed that POM_Bed transport only occurs during the 6 months of the high-flow season (t_transport=182.5 d). This approach allowed us to estimate POC_Bed without using near-bed velocities that were not available for all locations.

HW_North and HW_South both show an increase in the POM_Bed flux from the upper headwater locations (HW_North-2 and HW_South-2, respectively) to the lower headwater locations (HW_North-1 and HW_South-1, respectively), demonstrating the possibility of fast recruitment of POM_Bed in short distances. To determine changes in mass flux from the mountain front to the downstream reaches, we compare the combined headwater fluxes to those estimated at downstream sampling sites. The fluxes of POM_Bed at the upper Río Bermejo at HW_North-1 and Río San Francisco at HW_South-1 define the flux into the low-gradient portion of the river. We estimate the POM_Bed fluxes from HW_South and HW_North at 926–1138 t C yr⁻¹ and 112–188 t C yr⁻¹, respectively, for a total of 1038–1326 t C yr⁻¹ exported from the headwaters to the lowland reach (Fig. 8, Table 2). Downstream at LL_-1, we calculate a POM_Bed flux of 155–351 t C yr⁻¹, while at LL_-2 we estimate 19–27 t C yr⁻¹.

Figure 8 (a) Bar plot showing the order-of-magnitude estimates (± standard deviation) of the POM_Bed flux (t C yr⁻¹) at HW_South and HW_North. (b) Illustration of the source areas, transport, and fate of POM_Bed at the Río Bermejo from upstream recruitment downstream: sourcing from the headwater and the floodplain during wet-season erosion, partial attrition and transfer to suspended load causing a net decrease in POM_Bed, and onward transport of the remaining POM_Bed leading to (c) order-of-magnitude estimates of the lowland POM_Bed flux (t C yr⁻¹).

Table 2 Bedload sampling locations and yields from the field campaign in 2020, as well as the estimated flux of particulate organic carbon on the riverbed.

NA: not available. ^a Measured using ADCP. ^b Calculated using Eq. (3).

It is difficult to reconcile the minimum 66 % loss of C in bedload between the Bermejo–San Francisco confluence and LL_-1 (865 km transport distance) with the subsequent large apparent 10-fold reduction of the bedload C flux from LL_-1 to LL_-2 over a transport distance of only ∼220 km. The downstream increase in POM_Bed aggregates that dissociated into POM_Bed<1 mm could cause an underestimation of the POM_Bed>1 mm flux at LL_-2. The discrepancy may further be due to sampling bias (Turowski et al., 2013), bedload flux variability due to bedforms, the formation of OM waves at the channel bed, interspersed with relatively barren intervals (Heijnen et al., 2022), strongly sustained discharge (e.g., Rickenmann, 2018; Turowski et al., 2016; Reid et al., 1998), and channel geometry (Fogel and Lininger, 2023). Bedload sampling, particularly using Helley–Smith samplers, can be prone to high variability, and we suggest that our estimates are an order-of-magnitude approximation (Bunte et al., 2008) that likely underestimates the maximum POM_Bed transport rate. This could also cause the flux rates of HW_South that exceed the fluxes at HW_North, despite the higher erosive potential at HW_North. More agricultural activities in HW_South could also enhance surface erosion and with that OM input locally. Our approach assumes that the dimension of the POM_Bed layer and its individual particles are within the constraints of the funnel height of the sampler and that the samples and sampling points across each transect represent the cross-section of the channel. A larger sample size and sampled surface area, longer sampling times, and better understanding of distribution and dynamics of the POM_Bed layer could greatly enhance the accuracy of the sample set and flux estimates. Despite these uncertainties, this is the first estimate of its kind and shows that river POM fluxes may be underestimated without considering the bedload OM flux.

The inferred POM_Bed grain size reduction during transport, similar to clastic sediment (Attal and Lavé, 2009), contributes to the overall suspended sediment yield of the river. The total suspended organic carbon flux is $\sim \mathrm{1.85}\times {\mathrm{10}}^{\mathrm{5}}$ t C yr⁻¹ at the Bermejo–San Francisco confluence (Repasch et al., 2021), suggesting that the estimated POM_Bed carbon flux near the confluence is less than 1 % of the total carbon load. The Río Bermejo exports $\sim \mathrm{2.24}\times {\mathrm{10}}^{\mathrm{5}}$ t C yr⁻¹ in suspension downstream to the Río Paraguay, implying that about 0.39×10⁵ t C yr⁻¹ of suspended organic carbon is delivered to the lowland channel by lateral erosion (Repasch et al., 2021). If we take the downstream estimates of bedload C flux at face value and assume this loss transfers to the suspended load, a mass balance suggests that less than 1 % of the suspended load gain between the confluence and LL_-1 and LL_-2 could be due to grain size reduction of the coarse organic load. While our bedload carbon flux estimates are tentative, it is clear that this eye-catching mode of organic carbon transfer is small in comparison with the fluvial export of organic carbon in the suspended load of the Río Bermejo.

However, the Río Bermejo's suspended sediment yield is exceptionally high (Sambrook Smith et al., 2016). Assuming the loss between the HW_South-1+HW_North-1 and LL_-1 as well as LL_-2 transfers completely to the suspended load, 79 % and 98 %, respectively, of POM_Bed would transfer into suspension while simultaneously recruiting additional bedload from the floodplain. POM_Bed in rivers and sedimentary deposits could contribute substantially to the overall flux in river systems with lower suspended sediment yield and where bedload dominates the fluvial sediment flux (Turowski et al., 2016) or in highly erosive headwater streams with short transport distances from recruitment to subsequent deposition and burial (Blair and Aller, 2012; Hilton et al., 2011).

This coarse particulate OM may also have a higher probability of preservation and rapid burial in depositional basins, as its particle settling velocity is higher than smaller particles and may be less prone to resuspension and oxidation. Future work should aim to enhance our understanding of the significance of POM_Bed export and burial over varying transport length scales.

6 Conclusion

In this study, we investigated the occurrence, recruitment mechanisms, source areas, and survival of POM_Bed during long-range transport to understand the implications for the terrestrial organic carbon cycle. We found the persistent occurrence of POM_Bed along a 1300 km section of the Río Bermejo in northern Argentina from the headwaters to its confluence with the Río Paraguay.

Our results provide evidence that the POM_Bed originates from erosion of fresh, terrestrial organic debris from the local floodplain, as well as distal headwater sources, and is composed of a heterogenous mix from catchment-wide sources, dominated by C₃ plant input, and local C₄ point sources. POM_Bed can remain geochemically unaltered over long fluvial transport distances due to fast transport. The high geochemical variability within and between the bedload sampling transects (Fig. S2) does not allow a quantitative unmixing of POM_Bed sources; however, this variability indicates that POM_Bed may travel in parcels, each representing an erosive event that delivered fresh OM to the river, where it was subsequently waterlogged and translocated to the channel bed. Our data suggest coherent and continuous transport of individual OM parcels at the riverbed such that POM_Bed is never fully mixed longitudinally in a river system. Although POM_Bed appears to survive long-range fluvial transport, it can mechanically break down during transport, contributing to the overall suspended sediment load of the river. After physical breakdown, the fate of POM_Bed does not lay exclusively in the suspended fraction. The geochemical proxies for the Río Bermejo soil and bank sediment suggest it can also be stored on floodplains and undergo partial mineralization. Changes in hydraulic settings can lead to particle settling, and high suspended sediment yields enhance the burial efficiency and promote drawdown of atmospheric CO₂ over longer timescales.

The decreasing POM_Bed flux downstream also indicates that the loss of POM_Bed exceeds the recruitment from the floodplain and therefore that local ecosystem productivity is likely not the main factor controlling POM_Bed genesis at the Río Bermejo. Rather, the main control is likely to be the preconditioning of OM such that it can rapidly absorb water, increasing its density and its settling velocity, allowing it to sink (Hoover et al., 2010). Overland flow is also necessary to facilitate OM erosion and transport to the active channel, and sufficient river flow velocity is required to maintain transport at the bed (Galy et al., 2015). Likewise, the ratio of POM_Bed export and potential burial versus abrasion and transport as suspended sediment likely depends on plant type, climate, river morphology, hydrodynamics, and transit time (Hoover et al., 2010).

We found that transport of POM_Bed is not a major contributor to carbon burial on short timescales at the Río Bermejo, but ongoing floodplain recruitment contributes to the genesis of POM_Bed, and fluvial transport could export POM_Bed from the headwaters to the ocean on short timescales, representing a hitherto under-investigated mechanism for fluvial organic carbon export and burial. Bedload transport can convey OM to downstream basins with subsequent burial over millennial timescales; however, this is possibly more significant for shorter transport distances. The abundance and magnitude of bedload is in general highly variable and notoriously difficult to measure, and it remains challenging to quantify the total amount of POM_Bed due to stochasticity of bedload transport and the heterogeneity of POM_Bed abundance. Our approach is a first step to evaluate the origin and fate of POM_Bed during long-range fluvial transport in a natural setting. Further experimental and field studies are necessary to improve our understanding of the fraction of POM transported as bed material versus in the water column.