Integrating Suspended Sediment Flux in Large Alluvial River Channels: Application of a Synoptic Rouse‐Based Model to the Irrawaddy and Salween Rivers

A large portion of freshwater and sediment is exported to the ocean by a small number of major rivers. Many of these megarivers are subject to substantial anthropogenic pressures, which are having a major impact on water and sediment delivery to deltaic ecosystems. Due to hydrodynamic sorting, sediment grain size and composition vary strongly with depth and across the channel in large rivers, complicating flux quantification. To account for this, we modified a semi‐empirical Rouse model, synoptically predicting sediment concentration, grain‐size distribution, and organic carbon (%OC) concentration with depth and across the river channel. Using suspended sediment depth samples and flow velocity data, we applied this model to calculate sediment fluxes of the Irrawaddy (Ayeyarwady) and the Salween (Thanlwin), the last two free‐flowing megarivers in Southeast Asia. Deriving sediment‐discharge rating curves, we calculated an annual sediment flux of 326−70+91 Mt/year for the Irrawaddy and 159−51+78 Mt/year for the Salween, together exporting 46% as much sediment as the Ganges‐Brahmaputra system. The mean flux‐weighted sediment exported by the Irrawaddy is significantly coarser (D84 = 193 ± 13 μm) and OC‐poorer (0.29 ± 0.08 wt%) compared to the Salween (112 ± 27 μm and 0.59 ± 0.16 wt%, respectively). Both rivers export similar amounts of particulate organic carbon, with a total of 1.9−0.9+1.4 Mt C/year, 53% as much as the Ganges‐Brahmaputra. These results underline the global significance of the Irrawaddy and Salween rivers and warrant continued monitoring of their sediment flux, given the increasing anthropogenic pressures on these river basins.


Introduction
Rivers are the main conduits of dissolved and particulate matter from the continents to the oceans. Accurate quantification of material exported by rivers is thus often the most reliable and efficient way to constrain such key processes as continental erosion, chemical weathering, and organic carbon cycling (e.g., Gaillardet et al., 1999;Galy et al., 2015;Horan et al., 2019;Meybeck, 1987;Viers et al., 2013;West et al., 2005), leading to an improved understanding of the long-term controls on Earth surface conditions (e.g., Berner & Kothavala, 2001;France-Lanord & Derry, 1997;Godderis et al., 2009;Hilton et al., 2015;Mackenzie & Garrels, 1966;Maher & Chamberlain, 2014), as well as the anthropogenic perturbation of these processes (e.g., Allison et al., 2007;Syvitski & Kettner, 2011;Wilkinson & McElroy, 2007). On a global scale, the world's 30 largest rivers by discharge are estimated to account for ∼50% of all freshwater and ∼25% of all particulate matter export to the ocean (Milliman & Farnsworth, 2011). Southeast Asian rivers in particular dominate the global sediment flux, delivering about 2/3 of the supply to the ocean, due to a combination of active tectonics and monsoonal climate (Milliman & Farnsworth, 2011). The sediment fluxes of the Ganges-Brahmaputra, Mekong, Irrawaddy, and other major Southeast Asian rivers maintain extensive and fertile deltas, supporting large natural and agricultural ecosystems-the primary food source for several hundred million people. In addition, the tropical monsoonal climate enables high net primary productivity and efficient export and oceanic burial of biospheric carbon-an important sink for atmospheric CO 2 (e.g., Galy et al., 2007Galy et al., , 2015Hilton et al., 2008). Constraining the sediment and particulate organic carbon flux of large Southeast Asian rivers can significantly reduce uncertainties in the global exogenic carbon cycle, helping determine both the importance of natural feedback processes, as well as the scale of human perturbation in these river basins.
Accurately measuring the total sediment flux and its mean physicochemical composition is difficult in large alluvial river channels due to hydrodynamic sorting of sediments, which results in strong gradients in sediment grain size, concentration, and mineral type with depth (Jordan, 1965;Meade, 1994;Rouse, 1950;Vanoni, 1980). Although turbulent shear forces affect all particles equally, heavier (larger and denser) particles have higher settling velocities (Dietrich, 1982;Rouse, 1950). Suspended sediment concentration (SSC) at the surface is therefore not representative of the total sediment flux, which may be assessed by collecting discrete instantaneous samples at different depths, or by collecting a single depth-integrated sample, where the sampler is filled at a constant rate while being vertically lifted through the water column. However, it is often unclear how representative such single depth-integrated samples are, as the quality of integration strongly depends on sampler geometry, the speed at which the sampler is lifted through the water column, and the ability to maintain isokinetic sampling conditions (e.g., Murray Hicks & Gomez, 2016). The point-sampling approach has a major advantage, in that it allows an empirical calibration of sediment concentration as a function of flow conditions specific to each sample in the river reach of interest, potentially enabling the mapping of sediment load synoptically (with depth and across the river channel).
To date, most sediment flux and composition estimates of large rivers still rely on surface samples, with the notable exceptions being the Amazon and its major tributaries , Ganges , Changjiang (Guo & He, 2011), Mekong (Darby et al., 2016), Huanghe (Wang et al., 2007), Orinoco (Meade, 1994), and Mississippi (Meade & Stevens, 1990) rivers, which all have estimates derived via depth-and cross-channel sampling. A previously reported Irrawaddy River flux is also based on depth sampling, however, primarily from data collected in the 19th century using techniques that have since been significantly refined (Gordon, 1880;Robinson et al., 2007); see discussion below. All of the above-mentioned point-sampling studies of large rivers have revealed large variations in sediment concentration and composition in the river channel, indicating the need for depth (and lateral) sampling to obtain accurate estimates of sediment concentration and flux.
With the advent of Acoustic Doppler Current Profiler (ADCP) technology, it is now relatively simple and routine to measure flow velocity distribution in two dimensions (laterally and with depth) with submeter resolution in large river channels (e.g., Parsons et al., 2013;Thorne & Hanes, 2002;Yorke & Oberg, 2002). As a result, a number of attempts have been made to obtain a fully parametrized law for hydrodynamic sorting, which would allow the use of flow velocity data to predict sediment distribution across a river channel, with the need of just a few reference point samples. These attempts have revealed that the original Rouse model (Rouse, 1950) is unable to properly parametrize sediment distributions as function of velocity and depth, whether in large rivers Lupker et al., 2011) or in flume experiments (Muste et al., 2005, and references therein). The possible reasons are the complex distribution of particle sizes and shapes , particle aggregation due to organic matter , and the complex variation of the water and sediment diffusivity coefficients with sediment concentration (Muste et al., 2005;Pal & Ghoshal, 2016).
As an alternative, a number of indirect (surrogate) methods to determine riverine suspended loads, relying on optical and acoustic detection of sediments, have been tested (e.g., Armijos et al., 2017;Gray & Gartner, 2009;Vergne et al., 2020). In particular, ADCP instruments determine water flow velocity by using the acoustic echo from suspended particles, potentially allowing the simultaneous quantification of SSC with depth and across the river channel with high resolution (e.g., Thorne & Hanes, 2002). ADCP backscatter signal was successfully calibrated to calculate sediment flux of the Mekong River (Darby et al., 2016) and, more recently, the Paraña River (Szupiany et al., 2019). A number of complications have so far limited the applicability of this approach, however. First, acoustic instruments have variable sensitivity to different particles, most strongly impacted by grain size. Therefore, a single-frequency instrument is often unable to capture SSC variations in large rivers with complex, often multimodal particle size distributions and/or variable hydrodynamic conditions (e.g., Latosinski et al., 2014). Second, the calibration is typically instrument-specific such that raw data between two instruments (even of the same model) may not be comparable, requiring individual calibration for each acoustic instrument. Emerging multifrequency acoustic backscatter methods that are sensitive to a range of particle sizes show great promise but are still in their infancy (Vergne et al., 2020).
As a result, a hybrid empirical-theoretical approach based on the Rouse equation (Rouse, 1950; see section 4.1) has emerged as a robust way to quantify suspended sediment flux and chemical composition in large alluvial rivers with complex particle size distributions and/or highly variable hydrodynamic conditions Lupker et al., 2011). Instead of attempting to calibrate acoustic or optical sensing instruments, or to determine particle settling velocities for a fully theoretical prediction of SSC, point depth samples are collected to empirically calibrate the SSC-depth relationship under known hydrodynamic conditions (determined using ADCP). This approach assumes that instantaneous point samples are representative of equilibrium conditions (i.e., there is no net sediment suspension/deposition within the immediate channel reach). Any resulting error due to short-term turbulent fluctuations (e.g., Diplas et al., 2008) can be mitigated by collecting and averaging a larger number of samples (keeping in mind logistical constraints). It is also important to avoid reaches with significant changes to local hydrodynamic conditions, such as confluences, bifurcations, and meanders. This empirical calibration is repeated under different hydrodynamic conditions, which enables the construction of a SSC-discharge rating curve. Lupker et al. (2011) have demonstrated how point depth-sampling coupled with ADCP velocity measurements can enable more robust estimates of sediment flux, especially in kilometer-scale wide river channels with complex hydrodynamics and large lateral variations in flow velocity and sediment flux.
Here, we present an alternative approach to empirically calibrating the Rouse equation describing the SSC versus depth versus flow velocity relationship and apply this framework to the Irrawaddy and the Salween rivers in Myanmar. In contrast to previous efforts, this method makes fewer averaging assumptions and allows us to synoptically map high-resolution spatial variations in sediment concentration and composition both across the river channel and with depth. We use this approach to provide new estimates of the sediment and particulate organic carbon export flux by the Irrawaddy-Salween river system and compare them to values obtained using simple averaging approach, as well as previously published estimates.

Study Site
The Irrawaddy, also known as Ayeyarwady, and the Salween, also known as Thanlwin (Figure 1a), are believed to be among the largest rivers in terms of water and sediment flux globally, although previous data are scarce (Chapman et al., 2015;Furuichi et al., 2009;Robinson et al., 2007). The headwaters of the Irrawaddy originate in the southern margin of the eastern Himalayan Syntaxis. It runs for about 2,000 km, spanning the whole length of Myanmar and forming a large delta distributary network in the south prior to discharging into the Andaman Sea, with a basin surface area (taking topographic roughness into account using a 90 m resolution DEM) of 437,000 km 2 . The Salween originates in the Tibetan Plateau, traverses the eastern Himalayan Syntaxis, and flows south across the Shan Plateau in southeastern Myanmar. It has a length of ∼2,800 km and a basin surface area of 283,000 km 2 (Figure 1a). The Irrawaddy basin has a large central (relatively dry) valley, with a mean and maximum elevation of 862 and 5,798 m, respectively, and a median slope of 7.1°. In contrast, the Salween catchment is steep and narrow for such a large basin, with a mean and maximum elevation of 3,515 and 6,857 m, respectively, and a median slope of 16.4°.
Both river basins are composed of a wide variety of sedimentary, igneous, and metamorphic rocks, ranging from Pre-Cambrian to Cenozoic in age and transposed by a complex network of sutures and faults (e.g., Licht et al., 2013;Mitchell et al., 2012;Najman et al., 2020;Searle et al., 2007;Westerweel et al., 2019;Zaw et al., 2017;Zhang et al., 2018). The climate of both basins is dominated by the southwest Asian monsoon (and to a lesser degree the northeast monsoon), with most precipitation and discharge taking place in June through September (Zaw et al., 2017). Mean annual precipitation rates vary from <800 up to >4,000 mm/year within the Irrawaddy basin, depending on the location (e.g., Chen et al., 2017;Sein et al., 2018). Most water to both rivers is supplied by the monsoon precipitation, with additional (unquantified, but likely minor and further diminishing) inputs from mountain glacier melt and snowmelt in the north.
In terms of water and sediment flux and their chemical composition, the Irrawaddy and the Salween have very little data available compared to other Asian megarivers, largely due to historically difficult access to the country of Myanmar, which contains the major portions of both catchments (Figure 1a). The little data that are available point to the Irrawaddy-Salween system being a globally significant source of sediment and particulate organic carbon (POC) to the ocean, but these estimates have a large uncertainty (Bird et al., 2008;Furuichi et al., 2009;Robinson et al., 2007). Recently, Garzanti et al. (2016) carried out a study of detrital composition of bank sands across the length of the Irrawaddy River, providing a first estimate of relative sediment source-partitioning within the basin. Both Irrawaddy and Salween have remained largely undammed, with free-flowing mainstems (Grill et al., 2019) and only several small dams on minor tributaries within Myanmar (Hennig, 2016). However, the relatively small Chinese portion of the Irrawaddy basin has been extensively dammed, with a total generation capacity of around 6 GW. Additionally, over 40 dams, ranging from small to very large (>5 GW), have been announced and are either in planning or construction stage on the two rivers, with a total capacity of more than 45 GW (Hennig, 2016;Kirchherr et al., 2017;Lazarus et al., 2019), which will significantly alter their water and sediment discharge dynamics. In addition, Southeast Asian river sand is a major construction resource that is often unsustainably dredged and becoming increasingly scarce, resulting in bank erosion and collapse downstream and condemning low-lying river deltas to seawater intrusion and inundation Bendixen et al., 2019;Hackney et al., 2020;Kondolf et al., 2018;Xiqing et al., 2006). Collectively, damming, sand mining, and climate change have already impacted sediment delivery to the Irrawaddy delta, altering its extent and morphology (Chen et al., 2020). These pressures are set to only increase in the future, with potentially severe negative consequences for downstream ecosystems and communities. It is therefore crucial to establish a baseline of the current sediment flux and composition, so that any impact from potential future environmental change can be accurately assessed.

Discharge Measurements Using ADCP
Flow velocity measurements and sediment samples of the Irrawaddy-Salween rivers were collected during two monsoon seasons, in August 2017 and 2018, and two dry seasons, in February 2018 and May 2019. Both rivers were sampled just upstream of their delta distributary networks ( Figure 1). Flow velocity was measured using an Acoustic Doppler Current Profiler (ADCP) Rio Grande II (1,200 kHz) made by Teledyne Instruments, deployed on a moving boat. The ADCP was attached on a rigid frame close to the bow, in a down-facing orientation, and the transducer submerged at 40-60 cm depth. Data were collected while the boat crossed the river perpendicular to the flow. Boat position during the transect was recorded using an external GPS unit with horizontal accuracy of <5 m. Between 1 and 5 such transects were collected, depending on the site, with discharge reproducibility typically better than 6%, in agreement with previous applications of moving-vessel ADCP (e.g., Szupiany et al., 2007). Additional flow velocity data were collected during suspended sediment sampling in August 2018 and May 2019. Flow velocity profiles were recorded while the boat was drifting with the current, and 10 s of velocity data immediately before sample collection was averaged to represent the local, instantaneous hydrodynamic conditions affecting each suspended sediment sample. The averaged flow velocity profiles were used to calculate the shear velocity associated with each suspended sediment sample (Text S1 in the supporting information).
ADCP data were collected and initially processed using WinRiver II software. The data were then exported and further processed using Velocity Mapping Toolbox (Parsons et al., 2013). Using multiple river cross-sectional transects, a mean cross-section (MCS) was created for each sampling date (Figure 2), ensuring it was perpendicular to river channel, and calculating the average stream-wise flow velocity field across the river channel (Figures 1b and 1c). The data were then additionally processed in MATLAB 2019b, interpolating data gaps and removing erroneous outlier data (e.g., due to excessive pitch and roll) and extrapolating to the river surface (above ADCP transducer) and bottom (below sidelobe interference) using inpaint_nans function (D'Errico, 2018). This method uses a spring model to solve a sparse matrix of the measured velocity values and aims to extrapolate using a constant function. It was found to yield a more realistic velocity structure in the extrapolated regions, with high similarity to nearby bottom regions where ADCP data were available (Figure 2), in contrast to Laplacian methods that provide a smooth linear extrapolation (not shown). In regions with available data, it is apparent that flow velocity decreases significantly only in the near-bottom boundary region (i.e., following a power law with depth). In most cases, this low-velocity region is significantly thinner than the extrapolated region, and therefore, the model results presented are relatively insensitive to the extrapolation method. Both the raw and the extrapolated flow velocity data are freely available in an online data repository (see the Data Availability Statement).

Sediment Sample Collection and Processing
Depth profiles of suspended sediments were collected in August 2017, August 2018, and May 2019 at both sites. Only surface sediment samples were collected in February 2018. Sediment samples were collected at various depths using a modified 8.5 L capacity Van Dorn-style depth sampler (a horizontally oriented Perspex acrylic tube open at both ends, with pneumatically triggered doors, modified from Wildco, USA). Depth was determined either from measured rope length (August 2017) or a pressure transducer (August 2018 and May 2019). Approximately 30 kg of metal weights (hammer heads) were attached below the sampler to ensure vertical position of the sampler relative to the boat. The samples were collected while the boat was drifting with the flow, which, together with the open, wide (15 cm diameter) tube sampler design, ensures approximately isokinetic conditions. However, the sampler is not hydrodynamically shaped, and possible horizontal rotation could result in sampled sediment not always being representative of isokinetic conditions. Once at the required depth, the sampler doors were pneumatically shut using a bicycle pump. Additional bedload samples were collected by dredging river bottom sediments using a weighted metal bucket.
Samples were collected into 10 L sterile polyethylene bags, ensuring complete transfer of all sediment particles. The bags were weighed, and the samples filtered within 24 hr using 0.2 μm polyethersulfone (PES) membrane. The sediments were immediately washed off the filter and into an opaque glass jar, using filtered river water collected at the same site. The samples were kept sealed in the dark during transport back to the lab (between 1 and 2 weeks). They were then allowed to settle and were decanted (except very clay-rich samples), followed by freeze-drying using a Thermo Scientific ModulyoD freeze dryer. Suspended sediment concentration was calculated by dividing the dried sample weight by the weight of the total water sample prior to filtration, ignoring the <1% error due to sediment mass (<10 g/kg) in the original sample.
Particle size distributions of dried samples were measured using a Malvern Mastersizer 2000 laser diffractometer, at a 20 bin resolution ranging between 0.35 and 2,000 μm. Each sample amount was adjusted to achieve 10-20% obscuration and ranged from 50 to 5,000 mg, depending on the coarseness. Each sample was dispersed in tap water and sonicated for 2-5 min until grain-size distribution appeared stable. Each measurement was repeated three to five times. Typical uncertainty was less than 10% for each grain size bin, with most of the uncertainty due to subsampling errors of the coarse particles.
To measure the organic carbon concentration (weight %OC), carbonate was removed from the samples by a liquid HCl phase, within capsules with no rinse step (Komada et al., 2008). In detail, crushed sediment powders were weighed (∼5-10 mg sample for suspended sediments and 20 mg for bedload, attempting a target mass of organic carbon of ∼100 μg C) into 8 × 5 mm silver capsules that had previously been combusted (450°C for 4 hr, within 3 days of processing) and loaded open into a PTFE sample tray. Around 50 μl of 1 M HCl was added to each capsule, with the liquid reactant evaporated at 65°C to dryness in an oven. Acid addition and drying were repeated three times in total. Capsules were folded closed and analyzed by EA-IRMS at Elemtex with a range of international calibration standards and external standards (IAEA 600, IAEA CH3) and to check for full carbonate removal (NCS-DC73319). Measured %OC values were corrected for a full procedural blank (<5% of the sample carbon mass), and repeat measurements of samples and external standards had a precision of 0.05%.

Results
Water discharge measurements were performed for both rivers at the peak of the monsoon season and in middle and late dry season and therefore span about an order of magnitude range in discharge (Table 1; 3,000-42,100 m 3 /s for the Irrawaddy and 1,800-14,300 m 3 /s for the Salween). Importantly, these values bracket almost the full range of monthly mean discharge for both rivers (Table S4), allowing us to interpolate the results of this study for each month, yielding long-term average sediment composition and annual flux (see discussion below and Text S3).
The measured suspended sediment concentrations (SSC) ranged from 55 to 5,500 mg/L in the Irrawaddy and 47 to 10,200 mg/L in the Salween (all individual sample details and measured values are given in Table S1). The median grain size (D 50 ) ranged from 5 to 150 μm in the Irrawaddy and 8 to 130 μm in the Salween. The most concentrated (and coarsest) samples were collected during the monsoon and typically closer to the channel bottom, indicating the influence of hydrodynamic sorting. However, a significant number of coarse, high-concentration samples in both rivers were collected at mid-depth ( Figure 3). Because our depth sampler collects instantaneous samples without time averaging, the variable vertical dispersion of sand in our samples reflects the complexity of hydrodynamics in these rivers (e.g., non-steady-state turbulent sediment suspension events, secondary flow, and bedform effects). It is also possible that some of the scatter is due to a departure from truly isokinetic conditions during sampling, for example, due to sampler rotation. As discussed above, this complexity prevents simple spatial averaging with depth or across the river channel to calculate the total sediment flux and requires a fully spatially resolved sediment transport model (section 4).
The suspended particulate organic carbon concentration (wt%OC) ranged from 0.10% to 0.97% in the Irrawaddy and from 0.23% to 1.08% in the Salween (Table S1). As in many other rivers, most organic carbon is associated with finer particles, and sediment %OC is closely correlated with median sediment grain size ( Figure 4). This relationship can be used to convert the spatial D 50 distribution into %OC, and subsequently, the POC flux variations can also be calculated synoptically (with depth and laterally) across the river channel using the sediment transport model (section 4).

Model Description
River sediment is transported in suspension when turbulent shear stress (which can be expressed as shear velocity, denoted as u * ) is sufficient to overcome the particle settling velocity (e.g., Miller et al., 1977, denoted as w). Because turbulent shear stress affects all particles equally, whereas settling velocity depends on particle size and shape, the ratio of these two parameters can theoretically predict how the concentration of particles of a given size, shape, and density would vary with depth (Rouse, 1950): where C i is the sediment concentration in grain size class i, and z 0 is a reference height, defined here as fixed fraction of total water depth 0.001 · H, following Lupker et al. (2011), although we have also tested higher values (see Text S1). The sediment concentration at this reference height is C i 0 . The "Rouse depth," z r , is the sample depth z, nondimensionalized relative to the reference height z 0 , and total water column height H.
The power exponent in Equation 1 is commonly referred to as the Rouse number: The value of R i is dependent on particle settling velocity w i of sediment grain size i, the ratio of sediment and water momentum diffusion coefficients β, and shear velocity u * , which is representative of the boundary shear stress and can be calculated from depth-averaged flow velocity (see Text S1 and Equation S2); κ ¼ 0.41 is the von Kármán constant. The higher R i , the stronger the increase in sediment concentration with depth.
Attempts to obtain R i from fully theoretical considerations have so far been unsuccessful, due to a number of reasons. First, it is difficult to accurately determine particle settling velocity, especially for natural sediments composed of mixtures of mineral and organic matter of variable density and shapes (Dietrich, 1982), with potential particle aggregation adding further complication . Second, while many simpler treatments take β to be equal to 1, experimental data have shown it to vary considerably  (Table S1).

10.1029/2020JF005554
Journal of Geophysical Research: Earth Surface with sediment concentration (Muste et al., 2005), likely the reason for the complex variations in β observed in real rivers . For these reasons, previous workers were unable to apply Equation 3 to large rivers, instead turning to empirical calibration of R i using measured variations in sediment concentration with depth (Equation 1) Lupker et al., 2011).
In these previous applications of the Rouse model to large rivers, Equation 1 was used to either obtain one average R i across a river channel, effectively averaging laterally Bouchez, Métivier, et al., 2011), or applied to depth profiles collected under varying hydrodynamic conditions and establishing an empirical fit between depth-averaged sediment flux and u * . In other words, Bouchez, Métivier, et al. (2011) and  applied a single shear velocity value per cross section, therefore only integrating the geometry of the channel to calculate the flux, without modeling the lateral variation in hydrodynamic conditions. This approach worked well for Bouchez et al. because they were modeling very deep (up to 60 m) river channels in relatively straight sections of the Amazon River and its major tributaries, where the lateral variation in shear velocity was minimal. This, however, is not the case for many rivers with more complex channel cross-section morphologies, such as the lower Irrawaddy and Salween rivers studied here (Figure 2).
In contrast, Lupker et al. (2011) collected eight sediment sample depth profiles (n ¼ 3-9 per profile) at the same site on the Ganges River, but under strongly varying hydrodynamic conditions over the course of several years. They then applied Equation 1 individually to each depth profile, obtaining a vertically integrated sediment flux, relating it to local u * , and then using this relationship to laterally and temporally extrapolate the vertically integrated sediment flux. While robust, this approach requires a large number of suspended sediment samples and was enabled by a continuous field effort over the period of several years and is therefore not ideal for the smaller sample set of our study.
Here, we employ a different approach from these previous studies to address the highly dynamic flow conditions of the rivers studied here, while using a smaller number of sediment depth samples. We do this by explicitly factoring u * out of the fitted exponent in the Rouse equation: 10.1029/2020JF005554

Journal of Geophysical Research: Earth Surface
where z r is calculated from sample depth recorded during collection (Equation 2), u * is calculated from the depth-integrated flow velocity during sample collection (see Text S1 for details), and C i 0 and b i are fitted parameters (obtained separately for each grain size bin i).
Because b i is strongly dependent on sediment grain size, and grain-size distribution is known to vary widely with depth and hydrodynamic conditions in large rivers, measured sediment concentrations are divided into five grain size bins (i ¼ 0.2-4, 4-16, 16-63, 63-250, and 250-2,000 μm, chosen to minimize the number of C i values below detection, while retaining separation between particles with significantly different settling velocities). Equation 4 is then fitted individually to each set of C i values ( Figure 5; see Text S1). The near-bed equilibrium concentration C i 0 also varies as a function of u * (García, 2008). Even though this relationship is not explicitly represented in Equation 4, this equation constrains the C i =C i 0 ratio as a function of both depth and u * , and therefore, the relationship between C i 0 and u * is folded into the empirically determined parameter b i .
The empirically calibrated C i 0 and b i values can then be applied to ADCP-measured primary streamwise velocity data to calculate and map high-resolution variations in sediment concentration C i with depth and across the river channel ( Figure 6). Combining the five C i values also yields the variation in sediment grain size across the channel (Figure 6). The suspended sediment flux q s (kg m −2 s −1 ) distribution across the channel is then calculated for each ADCP data bin as q s ðz; xÞ ¼ ∑ i C i ðz; xÞ · u p ðz; xÞ: Integrating in both the lateral and vertical dimensions (i.e., summing up all ADCP bins) yields the total instantaneous suspended sediment flux Q S (kg s −1 ): q s ðz; xÞ · Aðz; xÞ: Here z and x are the bin coordinates in vertical (depth) and horizontal (lateral distance across the channel) direction, respectively; u p is primary streamwise flow velocity component; and A is the cross-sectional area of a given ADCP bin (e.g., 0.25 m × 0.5 m; variable depending on ADCP data resolution).
In summary, the method described here has certain advantages over previous applications of the point sampling approach to integrate sediment variation with depth in large rivers:   Table S1. 3. It enables a two-dimensional synoptic map of sediment concentration, flux, and grain-size distribution across morphologically complex river channels, where depth and flow velocity often show significant lateral variations ( Figure 6) and where averaging across the channel Morin et al., 2018;Santini et al., 2019) would likely result in significant errors of the calculated sediment flux and mean composition.
The above model applies only to sediment transported in suspension and does not include sediment carried as bedload below the reference height z 0 . To calculate the bedload flux, we adopted the semiempirical bedload transport equation The sediment modeling procedure described above was applied to the Irrawaddy and the Salween rivers separately, calculating the mean sediment concentration, grain size, and %OC distribution, as well as the total instantaneous sediment and POC flux for each of the four ADCP cross sections measured at each site. The results are summarized in Table 1, and the figures equivalent to Figure 6 for the other seven cross-sections are given in the supporting information.

Model Results: Instantaneous Sediment Flux and Composition
The calculated total instantaneous sediment flux ranged from 700 to 56,000 kg/s and from 400 to 25,000 kg/s for the Irrawaddy and the Salween, respectively (Table 1). The grain-size distribution was generally coarser and more variable in the Irrawaddy (D 50 range 10-43 μm) relative to the Salween (D 50 range 11-32 μm). Although the Irrawaddy discharge and sediment flux are about twice that of the Salween, due to the higher %OC of Salween sediments, the POC fluxes were similar in both rivers, ranging from 0.3 to 12 · 10 9 g C/day. The calculated bedload sediment flux ranged from 11 to 1,500 kg/s in the Irrawaddy and 6 to 740 kg/s in the Salween, representing only 1-3% of the total sediment flux in each case, regardless of the hydrodynamic conditions. These results agree well with the similarly small proportion (∼1.5%) of total sediment flux carried in the bedload in the Ganges River , as well as the Mekong River (Hackney et al., 2020), both similar in size to the Irrawaddy in their lower reaches. The total instantaneous (Table 1), monthly, and annual (see section 4.5) sediment flux values are all given as the sum of the suspended and the bedload sediment fluxes. The bedload POC flux is ignored, given that coarse sand contains low %OC ( Figure 4) and that the majority of sediment is carried as suspended load; this approximation should result in a negligible underestimation of the total POC flux.

Model Performance
To assess the performance of the model, the measured sample compositions can be compared to values calculated using the model at the equivalent locations (depth and lateral) in each channel cross section (see an example in Figure 6). The degree of misfit between measured and modeled values (represented as a mean relative standard error) was less than 5% for SSC and D 50 in both rivers, while the %OC relative standard  (Table 2), using the SSC rating curve determined by Robinson et al. (2007) and the season-average wt%OC determined by Bird et al. (2008). b Calculated as product of discharge and a simple mean of SSC and POC for all samples collected and analyzed on a given date, where n > 1 (Table S1). c Calculated using the Rouse modeling approach described in section 3. d Calculated as the relative difference between the simple mean-calculated value and the Rouse model-calculated value.

10.1029/2020JF005554
Journal of Geophysical Research: Earth Surface error was −35% for the Irrawaddy and +30% for the Salween (Figure 7). The higher and more systematic misfit of %OC is likely due to the considerably smaller number of data available to calibrate the model (Figures 4 and 7c) compared to SSC and D 50 and should be improved with additional analyses. We also note that this is not a strict test of the model, as it uses the training data set to assess the performance. A more rigorous assessment can be performed in the future against similar additional data sets (i.e., sediment samples coupled to ADCP flow velocity measurements) at these sites.
We propose that there are three main reasons for the misfit between the modeled and the measured values: 1. In some cases, large deviations from expected sediment sorting were observed, with several coarse, high-SSC samples collected at mid-depth (Figure 3), likely due to non-steady-state suspension events during sampling as discussed above. 2. There is some degree of mismatch between the ADCP velocity measurements (which integrate over an increasingly larger horizontal area with increasing depth) and the exact location of the collected sediment samples. 3. The location and the shape of the channel cross section varied slightly from year to year at both sites (Figures 1b and 1c; figures in the supporting information).
These factors inject substantial noise into our sample set, resulting in an offset between the sampled sediments and the local hydrodynamic conditions (represented by shear velocity) assigned to each sample (see Text S1). Finally, an additional source of uncertainty is the possible change in sediment supply to each river (e.g., seasonal hysteresis or interannual variations caused by landsliding or land use changes upstream) during the time span over which samples were collected for this study. However, such effects are typically local, and we expect them to be minor compared to the immediate turbulence-induced noise (Point No. 1 above), and to be mostly averaged out on the large basin-scale considered here. Ultimately, the spatial distribution of sampled sediment composition cannot be fully reconciled using a model that implicitly assumes constant sediment supply, constant channel structure, and equilibrium hydrodynamic conditions. Despite these complications, the sediment transport model presented here recovers the initial sample sediment composition for both the Irrawaddy and the Salween, without any large systematic errors (Figure 7).
Assessing the overall degree of agreement between measured and modeled values, the root-mean-square error (RMSE) was 23 μm for D 50 and 0.19 for %OC. For SSC, RMSE (calculated on a log scale) was a factor of 2.6, just a fraction of the total range, which spans a factor of >200 (Figure 7a). Additionally, the Nash-Sutcliffe efficiency (R 2 ) can be used to assess the performance of the model. Similarly to the coefficient of determination, R 2 ¼ 0 indicates that the model does no better than simply taking the mean of all measured  Table 1). The lines and envelopes show best fit and 68% confidence interval of the fit. The fitted rating curves and the goodness-of-fit statistics are given in the supporting information.

10.1029/2020JF005554
Journal of Geophysical Research: Earth Surface values, whereas R 2 ¼ 1 indicates that model reproduces all measured values perfectly (Nash & Sutcliffe, 1970). In our case, R 2 was 0.64, 0.53, and 0.45 for SSC, D 50 , and %OC, respectively, indicating reasonably high model efficiency. Below, we provide the final assessment of the model utility by comparing the flux and mean sediment composition values calculated here with estimates derived using simpler approaches.

The Need for a Hydrodynamic Sediment Transport Model
Our results indicate that, at least in the case of the Irrawaddy and the Salween, the sampled sediments frequently deviate from the expected Rousean behavior; that is, sampled sand concentration does not always increase with depth ( Figure 3). It is therefore reasonable to ask whether a Rouse-based hydrodynamic sediment transport model is required, and whether a simple averaging of all sediment samples, such as employed previously by Robinson et al. (2007) for the Irrawaddy, would yield flux and mean sediment composition estimates that are indistinguishable from the more complex hydrodynamic modeling approach employed in this study. To do this, we have compared the instantaneous sediment and POC fluxes, as well as mean grain size parameters calculated using the different approaches (including previously published rating curves and %OC values), shown in Table 2. Given that we collected sediment samples at roughly consistent depth percentiles (typically around 5-25-50-75-95% or 5-50-95% of total depth), as well as at several different lateral locations across the channel, we consider our sample set to be reasonably uniform in both dimensions of the channel cross section. Taking a simple average of the sampled SSC values and multiplying by the total ADCP-measured discharge has yielded sediment flux estimates that ranged from ∼40% lower during the dry season to ∼50% higher during the wet season, compared to Rouse model results for both rivers. Similarly, the mean grain size parameters (D 50 and D 84 ) were frequently overestimated or underestimated, depending on the particular cross section, reflecting the fact that simple-mean estimates fail to accurately account for sand transport in the near-bed region. Finally, using simple means of measured values significantly overestimated the POC flux by anywhere between 40% and 95% during the wet season for both rivers. Given the large size and discharge of the two rivers, this would result in a non-negligible error of riverine carbon export on a globally relevant scale. This comparison shows how crucial it is to accurately account for hydrodynamic sorting of sediments in large and morphologically and hydrodynamically complex rivers.
Although the chemical composition of the transported sediments is outside of the scope of this study, similar averaging errors can significantly affect the calculated fluxes of chemical elements, which are highly sensitive to particle grain size, such as silicon (mostly contained in the coarser sand grains) and aluminum and iron (mostly contained in clay-sized particles). These sorting bias effects were well exemplified and quantified on an element-by-element basis by Lupker et al. (2011) for the Amazon and the Ganges rivers, respectively. Given the importance of hydrodynamic sorting for the SSC and POC values in the Irrawaddy and Salween, we therefore expect similarly significant bias in elemental (and isotopic) fluxes, to be explored in follow-up studies.

Temporal Integration of Sediment Flux and Composition
The mean SSC and POC values calculated at the four different sampling dates and discharge (Q w ) conditions for each river allowed SSC-Q w rating curves to be constructed (Figure 8). Using previously published monthly Irrawaddy discharge data over a 31-year period   (Furuichi et al., 2009), we can calculate the monthly sediment and POC fluxes (Figure 9) and mean sediment concentration, grain size, and organic  Furuichi et al., 2009, andChapman et al., 2015, respectively). For discharge, the thick line represents the 31-year monthly averages for the Irrawaddy, whereas the Salween monthly discharge was calculated using the Irrawaddy/Salween discharge ratio determined in the wet and dry seasons in this study (see Text S3 for details). In (b) and (c), the thick line shows the best estimate with shaded area as the 68% confidence interval propagated through all calculations (see Text S3).

10.1029/2020JF005554
Journal of Geophysical Research: Earth Surface carbon content (Figure 10), which can then be summed to obtain long-term average annual values (Table 3). Unfortunately, other than our measurements presented here, the only Salween discharge data available cover a period between May and October in 2004, previously published by Chapman et al. (2015). The only annual discharge value available for the Salween is 210 km 3 /year given by Meybeck and Ragu (1997), which has since been used in a number of publications on rivers in Myanmar, as well as global compilations of water, sediment, and chemical fluxes (e.g., Chapman et al., 2015;Gaillardet et al., 1999;Robinson et al., 2007). For this reason, we used our ADCP-measured discharge values, along with the average monthly Irrawaddy discharge, to re-estimate the monthly discharge of the Salween in proportion to Irrawaddy discharge, yielding a revised annual Salween discharge of 149 km 3 /year (see Text S3 for details).
Applying the rating curves shown in Figure 8 to the monthly discharge time series, we are able to calculate the monthly suspended sediment and particulate organic carbon concentrations, median grain size (Figure 10), and the sediment and POC fluxes (Figure 9; all values given in Table S4). As expected, the sediment composition and flux vary by more than an order of magnitude in both rivers, with the coarsening of the transported sediment and the highest fluxes during the monsoon: monthly mean SSC ranged from 0.20 to 1.1 g/L in the Irrawaddy and from 0.22 to 1.6 g/L in the Salween, with annual flux-weighted means of 0.9 ± 0.2 and 1:1 þ0:5 −0:4 g/L, respectively (1σ uncertainty; Table 3). Overall, the Salween sediments are finer (D 50 from 11 to 25 μm, compared to the Irrawaddy's 10 to 42 μm, with flux-weighted annual means of 21 þ5 −4 and 28 þ6 −5 μm, respectively. It must be noted that these mean SSC and D 50 values are conservative estimates and mean %OC concentration is an upper-end estimate, because the monthly resolution used here, combined with nonlinear rating curves (Figure 8), likely results in the underestimation of sediment flux during high discharge events over shorter time scales.
Due to its lower discharge, the Salween sediment flux of 159 þ78 −51 Mt/year is about half of the Irrawaddy's 326 þ91 −70 Mt/year, with bedload comprising ∼2% of each. However, because organic carbon concentration in the Salween is about twice that of the Irrawaddy (0.59 ± 0.13% vs. 0.29 ± 0.08%), both rivers deliver a similar POC flux of ∼1 Mt C/year to the ocean. Of the annual suspended sediment flux, 37% and 24% is composed of sand (particles larger than 63 μm) in the Irrawaddy and the Salween, respectively.

Comparison to Previously Published Annual Flux Estimates
Compared to other major global rivers, prior to this study there existed very little modern data on the water and sediment discharge by the Irrawaddy and the Salween. The most significant data set was collected in the 19th century by Gordon (1880), presenting 10 years of discharge and suspended sediment measurements on the Irrawaddy at a location close to our sampling site at Pyay. More recently, Robinson et al. (2007) collected additional sediment depth samples and re-evaluated the original Gordon data set, determining annual estimates of water discharge of 422 ± 41 km 3 /year and sediment flux of 364 ± 64 Mt/year. Subsequently, Furuichi et al. (2009) used 31 years of discharge data published by the Department of Hydrology and Meteorology (DHM) in Myanmar (the same data set was used in this study) to calculate annual discharge of 379 ± 47 km 3 /year, where the uncertainty was given as 1 standard deviation of interannual variability and is therefore an overestimate of actual uncertainty on the long-term average, which we recalculate here  Table 1). The thick line shows the best estimate with a 1σ uncertainty indicated by the envelope. Details of calculations are given in section 4.5 and Text S3, and the calculated monthly values are given in the supporting information.
as 1 standard error of the mean, equal to 9 km 3 /year (Table 3) for the same 31-year period. Furuichi et al. (2009) further used a sediment rating curve for the Irrawaddy developed by DHM to estimate an annual sediment flux of 325 ± 57 Mt/year, in good agreement with our results. However, it must be noted that neither the sampling protocol nor the data used to establish the rating curve given in Furuichi et al. (2009) are publicly available. Lastly, the Irrawaddy basin-wide erosion rate determined here (Table 3) is in excellent agreement with cosmogenic 10 Be-based erosion rate of 0.27 mm/year (Garzanti et al., 2016).
Similarly, we revised the Salween sediment flux from 180 Mt/year previously estimated by Robinson et al. (2007) using their Irrawaddy sediment rating curve, down to 159 þ78 −51 Mt/year, using the first rating curve for the Salween, presented here. We note that discharge monitoring of the Salween is necessary to further improve this estimate.
Our determined annual POC fluxes are significantly lower than the values previously presented in Bird et al. (2008): 0.55-1.55 versus 2.2-4.3 Mt C/year for the Irrawaddy and 0.46-1.79 versus 2.4-3.4 Mt C/year for the Salween, a twofold-to-fivefold reduction in each case. It is partly explained by the reduction in water discharge estimates, but the main reason is the significantly lower %OC measured in this study (Table 3; also see the supporting information for individual sample values), compared to the values determined by Bird et al. (2008). One possibility is that this difference represents an actual decrease in %OC over the past decade. However, a change of this magnitude is difficult to defend, considering the large area of both river basins, and the fact that the difference is of similar order for both rivers. We suggest that this discrepancy is likely the result of sampling methodology differences between Bird et al. (2008) and the present study. Bird et al. (2008) used a 2 L horizontal Van Dorn sampler, collecting sediment samples at 1 m depth from the surface, mid-depth, and 1 m depth from the bottom, measuring OC of 1.1-1.6 wt% during high-discharge monsoon conditions, with similar values in both Irrawaddy (at Pyay) and Salween (at Hpa-An) and almost constant throughout the water column, suggesting negligible hydrodynamic sorting. This observation is in stark disagreement with both our results and those of Gordon (1880). Although it is difficult to determine the exact reason for this discrepancy, we speculate that sand may not have been adequately sampled by the smaller 2 L volume sampler used by Bird et al. (2008) (vs. our 8.5 L sampler, where we took extreme care to rinse out and collect all sand particles during sample transfer). This reinforces the idea that thorough depth sampling and sediment flux modeling that accounts for hydrodynamic sorting is crucial for accurate flux estimates in large rivers, especially for elements such as carbon, whose concentrations are strongly coupled to sediment grain size (Figure 4).

Global Significance of the Irrawaddy-Salween System
Globally, using the values presented in this study, the Irrawaddy and the Salween exhibit some of the highest sediment fluxes (fifth and seventh worldwide, respectively; Figure 11) and area-normalized sediment yields (third and fourth, respectively, among world's 30 major global rivers with annual discharge >100 km 3 year −1 as compiled by Milliman & Farnsworth, 2011, and lower only than the Fly and Brahmaputra rivers). Compared to the nearby Ganges-Brahmaputra system, which is the main conveyor of Himalayan erosion products to the ocean, the Irrawaddy-Salween system sediment yield is very similar, and sediment flux is about 46% that of Ganges-Brahmaputra. In comparison, the Mekong River, also originating in the eastern Himalayan Syntaxis, used to deliver ∼150 Mt year −1 (Milliman & Farnsworth, 2011), which has decreased to 87 ± 28 Mt year −1 (approximately two and approximately four times lower than the current fluxes of the Salween and the Irrawaddy, respectively) over the past several decades due to damming and changes in precipitation across the basin (Darby et al., 2016).
Although it is difficult to assess the global significance of the Irrawaddy-Salween system due to uncertainty of the global sediment flux in general, comparing to the estimate of Milliman and Farnsworth (2011), the two rivers are an important source of sediment to the ocean, delivering 2-3% of the estimated 19,000 Mt year −1 total sediment and 0.8-1.2% of the 200 Mt C year −1 total (biospheric and petrogenic) POC  export to the ocean. It must be noted, however, that current sediment flux estimates may be inaccurate for a number of large global rivers, where values are derived from sparse sample sets, often of surface sediments only, lacking the depth sampling and hydrodynamic data required to obtain robust values. The significance of our results is further underlined by the fact that the Irrawaddy and the Salween are some of the last large river basins still relatively unaffected by damming. Only a few small dams have been built on some minor tributaries of both rivers, with their mainstems flowing freely from source to outlet (worldwide, the only other megarivers with free-flowing mainstems are the Amazon and the Congo (Grill et al., 2019)). Currently, the main anthropogenic pressures on the Irrawaddy and the Salween river basins, such as deforestation, agriculture, and sand mining, are likely to be net erosive, temporarily enhancing the sediment flux (Syvitski et al., 2005). However, large dams are planned on both rivers, which, if built, will trap large amounts of sediment. If the majority of the current sand flux is captured by dams or sand mining, this could reduce the total sediment export of these rivers by up to a third. Although the Irrawaddy delta has grown over the past several decades (Chen et al., 2020), a reduced sand supply would have serious detrimental effects on this densely populated and highly productive agricultural area, potentially replicating the crisis currently unfolding in the neighboring Mekong River basin (Hackney et al., 2020). Our results presented here thus establish an important pre-dam baseline that will allow to identify and quantify any future changes to the sediment export by the Irrawaddy and the Salween rivers.

Conclusions
In this contribution, we have presented a new semi-empirical hydrodynamic Rouse modeling approach to synoptically predict the two dimensional distribution of suspended sediment concentration, its physicochemical composition (grain size and organic carbon content), and flux in large, turbulent rivers with geomorophologically complex channels. We have applied this model to obtain spatially and temporally integrated estimates of the sediment composition and export flux of the Irrawaddy and Salween rivers in Southeast Asia. In comparison to the model, flux estimates derived from using simple means of evenly spaced depth point samples can result in errors of up to 50%. This demonstrates that synoptic (i.e., spatially highly resolved) sediment transport modeling is crucial for the accurate quantification of sediment composition and flux in large river channels, where wide sediment grain-size distributions and variable hydraulic conditions result in complex sediment transport patterns. 10.1029/2020JF005554

Journal of Geophysical Research: Earth Surface
Using the approach outlined above, we have calculated a total sediment flux of 485 (68% confidence interval of 364-654) Mt/year and a particulate carbon flux of 1.9 (1.0-3.3) Mt C/year for the Irrawaddy-Salween system, accounting, respectively, for 2-3% and 0.8-1.2% of the total global riverine export to the ocean. These results represent a ∼20% and a 60-80% reduction of sediment and POC fluxes, respectively, compared to previously best estimates, which were partly based on 19th century data. While some of this difference may potentially be accounted for by actual changes in deforestation, land use, and other anthropogenic pressures in the river basins, we suggest that most of the difference is likely methodological, stemming from the use of a robust hydrodynamic sediment transport model in the current study. We expect that the methods and results described here, when combined with chemical and isotopic analyses of sediments at these and other sites in the Irrawaddy and the Salween basins, will enable a deeper understanding of the sediment provenance, erosion, and chemical weathering dynamics in the region, with the ultimate aim of fully constraining the regional organic and inorganic carbon cycle.
While the upstream sediment supply remains relatively constant, our calibrated Rouse-model fits presented here allow the use of ADCP data to predict the spatial distribution of SSC and POC across each river channel in the near future. In turn, our calibrated SSC rating curves allow the prediction of total sediment flux with varying discharge. However, given that a number of large dams are planned on major tributaries and mainstems of both rivers, sediment supply to their respective lower basins is expected to change, if and once these dams are constructed. In this case, active, depth sampling-based monitoring of sediment fluxes will be required to accurately quantify the changing sediment flux and composition. In this case, the results of our current study provide an important pre-dam baseline against which future changes can be evaluated.

Data Availability Statement
All measurement and final model results are tabulated in the main text and in the supporting information tables. All data presented in this study, including sediment composition, raw, and processed ADCP data and individual cross-section model results, are available on the NERC EIDC data repository (https://doi. org/10.5285/86f17d61-141f-4500-9aa5-26a82aef0b33).