A new large-scale suspended sediment model and its application over the United States

. Suspended sediment plays a vital role in the regional and global cycling of carbon and nutrients by carrying carbon 25 and nutrients from headwaters into lowland rivers and the oceans. Sediment transport through river systems is often fundamentally modified by human activities like such as reservoir management. However, a physically based representation of sediment transport is still missing in most existing earth Earth system System models Models (ESMs), which are essential tools for modeling and predicting earth system changes. Here, we introduce a multi-process river sediment module for ESMs, which includes: 1) hillslope soil erosion and sediment discharge into streams; 2) sediment transport processes through river 30 networks; 3) reservoir operation based on the inflows from upstream areas and water demand from downstream areas; and 4) sediment trapping by reservoirs. All model parameters are estimated a priori without calibration. We apply this new sediment modeling framework to the contiguous United States and validate it against historical observations of monthly streamflow and sediment discharges at 35 river gauges. The model reasonably well captures the long-term balance and seasonal variations of suspended sediment in large river systems. Furthermore, our model results show that suspended sediment discharge in managed 35 rivers is affected more by reservoirs' direct trapping of sediment particles than by their flow regulation. This new sediment module enables future modeling of the transportation and transformation of carbon and nutrients carried by the fine sediment along the river-ocean continuum to close the global carbon and nutrients cycles.

Abstract. Suspended sediment plays a vital role in the regional and global cycling of carbon and nutrients by carrying carbon 25 and nutrients from headwaters into lowland rivers and the oceans. Sediment transport through river systems is often fundamentally modified by human activities like such as reservoir management. However, a physically based representation of sediment transport is still missing in most existing earth Earth system System models Models (ESMs), which are essential tools for modeling and predicting earth system changes. Here, we introduce a multi-process river sediment module for ESMs, which includes: 1) hillslope soil erosion and sediment discharge into streams; 2) sediment transport processes through river 30 networks; 3) reservoir operation based on the inflows from upstream areas and water demand from downstream areas; and 4) sediment trapping by reservoirs. All model parameters are estimated a priori without calibration. We apply this new sediment modeling framework to the contiguous United States and validate it against historical observations of monthly streamflow and sediment discharges at 35 river gauges. The model reasonably well captures the long-term balance and seasonal variations of suspended sediment in large river systems. Furthermore, our model results show that suspended sediment discharge in managed 35 rivers is affected more by reservoirs' direct trapping of sediment particles than by their flow regulation. This new sediment module enables future modeling of the transportation and transformation of carbon and nutrients carried by the fine sediment along the river-ocean continuum to close the global carbon and nutrients cycles.

Introduction
Fluvial suspended sediment (referred to as suspended sediment) is fine-grained particles that may be diffused throughout the vertical column of rivers via turbulence and transported along rivers (Garcia (Ed.), 2008). It typically consists of fine sand, silt, and clay, which can absorb carbon and nutrients. Recently, increasing evidence showed that suspended sediment transport plays a vital role in regional and global carbon and nutrient cycling (Berhe et al., 2018;Lal, 2003Lal, , 2004Ludwig & Probst, 60 1996;Maavara et al., 2020;Van Oost et al., 2007) by carrying carbon and nutrients from soil storage pools into rivers and eventually oceans or lakes. For instance, suspended sediment is a major source of carbon and nutrients sustaining coastal wetlands and deltas (Kirwan & Megonigal, 2013;Nienhuis et al., 2020). While understanding the global carbon and nutrient cycling has been one of the primary goals of Earth System Models (ESMs), few ESMs have incorporated the representation of suspended sediment and the associated carbon and nutrient fluxes, e.g., particulate organic and inorganic carbon, nitrogen, 65 and phosphorous.
In the field of Earth system modeling, there are two primary reasons for the lack of suspended sediment in current ESMs.
First, most ESMs have been developed with a conventional assumption that the lateral carbon fluxes through the land-riverocean continuum are not significant, compared to the carbon dynamics within the larger carbon pool in the atmosphere, land, and ocean or to the vertical carbon fluxes between the atmosphere and the underlying land and ocean. This assumption has 70 now been refuted by data-based empirical studies (Berhe et al., 2018;Lal, 2003Lal, , 2004Ludwig & Probst, 1996;Maavara et al., 2020;Van Oost et al., 2007). Second, process-based representations of hydrological fluxes, such as overland flow generation and routing, channel routing, and reservoir operations that drive the suspended sediment fluxes, are not included in many ESMs. The recent development of a physically-based large-scale surface flow routing model, Model for Scale Adaptive River Transport (MOSART) Li, Leung, Getirana, et al., 2015;Li, Leung, Tesfa, et al., 2015;Voisin, Li, et al., 75 2013;Voisin, Liu, et al., 2013) is one of the few attempts to overcome this obstacle. MOSART explicitly incorporates the processes of overland flow routing across hillslopes, channel routing through tributaries and main channel networks, and water management such as reservoir regulation and surface water withdrawal, thus providing a solid physical basis to represent suspended sediment dynamics in both natural and managed rivers. Therefore, building on the MOSART framework, our first objective is to advance Earth system modeling by introducing a new suspended sediment module in ESMs. 80 In the field of sediment modeling, there have been several large-scale models of suspended sediment accounting for the spatiotemporal heterogeneity in sediment supply and transport processes (Czuba et al., 2017;Ferguson et al., 2015;Patil et al., 2012;Pelletier, 2012;Schmitt et al., 2016Schmitt et al., , 2018Tsuruta et al., 2018;White et al., 2014). However, few have been designed for regional or global applications with explicit accounting for the effects of reservoirs prevalent in most river systems (Nilsson et al., 2005). Among the exceptions are BQART (Syvitski et al., , 2005 and WBMsed (Cohen et al., 2013(Cohen et al., , 201485 Moragoda & Cohen, 2020). In BQART, a multivariable regression relationship has been established for each basin between the long-term suspended sediment load (hereafter, "load" refers to sediment discharge/flux) at the basin outlet and several parameters accounting for the average climate, lithology, and anthropogenic conditions at the basin scale

Modeling Framework 120
As shown in Figure 1, MOSART-sediment is developed on top of the water module of MOSART (denoted as MOSARTwater). Here we provide a detailed description of MOSART-sediment, but only a summary description of MOSART-water for completeness since more details of the latter can be found in our previous publications Li, Leung, Getirana, et al., 2015;Li, Leung, Tesfa, et al., 2015;Voisin, Liu, et al., 2013).

MOSART-water 125
MOSART-water is a spatially distributed model which takes surface and subsurface runoff time series generated from a land model as inputs and performs runoff routing across hillslopes (overland flow routing) and through channels (channel routing) . The fundamental spatial units of MOSART-water can be regular latitude/longitude grid cells or watersheds, and the former is used in this study. Within each spatial unit, surface runoff is first routed across hillslopes and then discharged along with subsurface runoff into a "tributary sub-network" channel before entering the main channel. The 130 spatial units are thus linked via routing through the main channel network, which is constructed in a scale-consistent way across different spatial resolutions based on the hierarchical river tracing method (H. Wu et al., 2011Wu et al., , 2012. The overland flow routing is modeled using the kinematic wave method, while the channel routing can be simulated using either the kinematic wave method or the diffusion wave method. All model parameters are physically based and available globally as described in our previous studies Li, Leung, Getirana, et al., 2015;Li, Leung, Tesfa, et al., 2015). In these 135 applications of MOSART-water at the regional to global scales, no parameter calibration was conducted.
MOSART-water includes a water management module applied at regional and continental scales (Li, Leung, Tesfa, et al., 2015;Voisin, Liu, et al., 2013). The water management module represents two typical water management activities: reservoir regulation/operation, and local water extraction/diversion. The focus of our reservoir regulation module is not on any individual reservoir but on hundreds to thousands of reservoirs and their impacts on regional to global processes for earth Earth system 140 modeling. It is almost impossible to obtain the actual operation rules of all the reservoirs. Nevertheless, some critical information of large reservoirs (e.g., storage capacity no less than 0.1km 3 ), such as the storage capacity and primary operating purposes, is available in some global datasets (Lehner et al., 2011). These large reservoirs can be classified into different categories based on their purposes, such as irrigation, flood-control, navigation, etc. For each reservoir category, we followed a generic algorithm to derive an operation rule for each reservoir based on the hydroclimate conditions in the upstream areas 145 and the water demand in the downstream areas (Li, Leung, Tesfa, et al., 2015;Voisin, Liu, et al., 2013). For example, for flood control purpose, reservoir release should follow a flat curve. ; fFor irrigation purpose, reservoir release should follow the seasonal variation of irrigation water demand.

MOSART-sediment
Riverine sediment load can be divided into two components: bed-material load and wash load (Biedenharn et al., 2006). 150 Wash load is the fraction of the moving sediment that cannot be found in the bed material with an appropriate fraction (Einstein et al., 1940;Lane, 1947) whereas bed-material load is the part that can be found substantially in the bed material. In other words, the bed material load is from the bed material and controlled by the local hydrodynamics and the bed material grain size. Wash load follows the fluid flow well without much exchange with the bed material and is thus controlled by the upstream supply. The bed-material load can be further divided into bedload and suspended part of bed-material load (suspended bed-155 material load) (Garcia (Ed.), 2008). Bedload refers to the coarse particles being transported in the vicinity of the river bed.
Suspended bed-material load refers to particles being transported in suspension within the vertical column of the river flow and they interact with river beds usually after a relatively long suspended path in the fluid flow. Wash load often refers to the very fine particles, i.e., grain size less than 0.062 mm (Lane, 1947), which are usually in suspension, well mixed with water, and rarely reside on the river bed. Therefore, total suspended load consists of suspended bed-material load and wash load. This 160 study focuses on total suspended load.
Because of the difference in their sources, bed-material load and wash load are modelled separately, and then the bedmaterial load is further divided into suspended bed-material load and bedload. In most rivers, particularly those large rivers, bed-material load consists of sand and gravel particles typically sourced from riverbeds. Wash load consists of silt and clay particles usually sourced from surface erosion at hillslopes. The transport processes of bed-material load and wash load are 165 distinct and thus can be treated separately. Bed-material load is controlled by local hydraulic conditions and sediment particle size. Wash load is not constrained by local hydraulic conditions since channel flow usually has sufficient turbulent kinematic energy to suspend and transport wash load particles (silt and clay) (Garcia (Ed.), 2008;Moodie et al., 2020).. Three primary components thus underpin MOSART-sediment: 1) a recently developed soil erosion module (Tan et al., 2018(Tan et al., , 2020(Tan et al., , 2021 for wash load supplied from hillslopes driven by overland flow, 2) a new riverine sediment module that 170 treats wash load and bed-material load processes separately, and 3) a new reservoir module for reservoir effects on sediment transport. Next, we briefly describe the soil erosion module and provide details for the riverine and reservoir modules.

Hillslope Sediment Processes
The soil erosion module in the E3SM land model (ELM-Erosion) deals with two major processes: 1) detachment of soil particles (clay, silt, and sand) from topsoil (or soil erosion) by raindrops or overland flow, and 2) delivery of the detached soil 175 particles along with overland flow from its location of detachment to rivers, during which the sand particles may deposit back to soil surface due to the limited transport capacity of overland flow. Here we only provide a brief description of these two processes and refer to more details in Tan et al. (2018).
The ELM-Erosion sediment yield model was developed based on the simple Morgan-Morgan-Finney (MMF) model (Morgan, 2001;Morgan & Duzant, 2008), which is better suited for capturing the heterogeneous sediment yield at spatial and 180 temporal scales relevant to ESMs (several to hundreds of square kilometers in space and tens of minutes to hours in time) than other well-known models (Tan et al., 2018). This choice was made by applying and comparing the eight most extensively used sediment yield models, such as the Water Erosion Prediction Project (WEPP) model, at 454 small catchments in diverse environments across the United States, Canada, Puerto Rico, U.S. Virgin Islands and Guam. In our test that each sediment yield model was deployed at the same hydrological, vegetation and soil conditions,, the MMF model outperformed the other 185 models incan simulatingrepresent 53% of the spatial variability of sediment yield over these catchments, which is much better than the other seven models, such as the Water Erosion Prediction Project (WEPP) model. Our test also showed that the ELM-Erosion model couldthe MMF model reproduced the observed daily sediment yield at six cropland-dominated catchments, with a relative error as low as 10% (Tan et al., 2018). By implementing the MMF model in ELM as ELM-Erosion and driving it by ELM simulated throughfall, leaf drainage, overland flow, leaf and ground cover, canopy height, and soil cohesion, we 190 further tested its performance at the continental scale. FurthermoreAs shown in Tan et al. (2020), the model canis found to well capture the spatial variability of soil erosion and sediment yield among the U.S. states and large river basins (Tan et al., 2020). The ELM-Erosion model can account for the impacts of rainfall and runoff fluxes at the event scale and the impacts of vegetation dynamics such as changing leaf area index and canopy heights at the monthly or longer temporal scales. The ELM-Erosion has three parameters for calibration: a scaling factor parameter for rainfall-driven erosion, a scaling factor parameter 195 for runoff-driven erosion, and a scaling factor parameter for sediment transport capacity. Following Tan et al. (2020), in this study, the ELM-Erosion model is run at the 1/8-th degree spatial resolution and the 30-min time step. For each grid cell, the ELM-Erosion has three parameters for calibration: a scaling factor parameter for rainfall-driven erosion, a scaling factor parameter for runoff-driven erosion, and a scaling factor parameter for sediment transport capacity. The calibration is conducteded in two steps: 1) the two erosion related parameters are calibrated at the state level over the conterminous U.S. 200 using the National Resources Inventory (NRI) soil erosion benchmark data (the grid cells in the same state would have the same parameter values); and 2) with soil erosion constrained, the transport capability parameter is calibrated at the basin scale using the pre-dam sediment yield data of large river basins.
Following the original MMF model (Morgan, 2001), Note that in the ELM-Erosion model does not calculate sediment yield for different, the distribution of soil particle size, such as clay, silt and sand particles, as done in the revised MMF model 205 is not considered. (Morgan & Duzant, (2008). revised the original MMF model (Morgan, 2001) to account for the sediment yield for clay, silt and sand particles, respectively. ButIt is because we found thatin their comparison that this modification did not improve the model performance over the continental scale while introducing additional complexity and computational cost (Tan et al., 2018). Therefore, we chose to refine the MMF model directly for our sediment yield modeling. In this study, we assume that sediment yield to rivers caused by soil erosion mainly consists of silt and clay particles, or in other words in the 210 form of wash load (Garcia (Ed.), 2008;Patil et al., 2012). This assumption is reasonable because During a storm event, even through detachment could happen to either clay, silt, or sand soil particles,. tThe transport capacity of overland flow for the detached soil particles across hillslopes decreases exponentially with increasing soil particle size (Morgan & Duzant, 2008) and . aA majority of detached sand particles will thus deposit back to the topsoil instead of entering the rivers. Therefore, it is Field Code Changed reasonable to assume that wash load (silt and clay particles) is mainly from soil erosion at hillslopes (Garcia (Ed.), 2008;Patil 215 et al., 2012).

Riverine Sediment Processes
The riverine sediment module focuses on total suspended load and simulates suspended sediment process on top of the hydraulic conditions simulated by MOSART-water as aforementioned. The suspended load consists of the wash load and suspended bed-material load, which are treated separately. 220 Wash load is primarily carried by channel flow. For the sake of simplicity, we consider that the net sediment exchange between channel and floodplain is secondary as comparing to the riverine sediment discharge and thus negligible. The effects of reservoir processes on wash load are however nontrivial since reservoirs fundamentally slow down river flow in their immediate upstream channels. We detail the reservoir effects in Section 2.2.3.
Existing sediment transport formulas derived in sedimentology are mainly applicable to bed-material load depending on 225 the river bed conditions and do not account for wash load (Garcia (Ed.), 2008). Most rivers can be roughly classified as gravelbedded or sand-bedded rivers (Garcia (Ed.), 2008). Gravel-bedded rivers are usually located in mountainous headwater areas, where gravel and coarser material (e.g., cobbles and boulders pebbles and bounders) prevail in bed-material load and move mostly in bedload mode. Sand-bedded rivers are mostly in the lowland regions where the bed-material load is dominated by sand and finer particles, i.e., suspended bed-material load (Lamb & Venditti, 2016). 230 Since this study focuses on suspended sediment, we mainly consider sand-bedded rivers and choose the sediment transport formula from those more suitable for sand-bedded rivers. The classic Engelund-Hansen equation (Engelund & Hansen, 1967;Ma et al., 2020) is typically applicable over sand-bedded rivers and can be further extended to finer-grained rivers (Ma et al., 2017. It is one of the simplest formulas of bed-material load (Brownlie, 1983), and is relatively easy to estimate its parameters a priori, which is critically important for sediment simulation in an Earth system model setting at the regional or 235 larger scales. Despite its simplicity, the Engelund-Hansen equation is one of the most accurate sediment transport formulas based on extensive validation against both laboratory and field observations (Brownlie, 1983;Crosato et al., 2011;Crosato & Mosselman, 2009;Darby et al., 2002;Geleynse et al., 2011;Kleinhans et al., 2008;Lee et al., 1997;Ma et al., 2020;Mosselman, 1998;Nicholas, 2013;Schuurman et al., 2013;Simon & Darby, 1997;Van Der Wegen & Roelvink, 2008;Wiele et al., 1996). We therefore adopt the classic Engelund-Hansen equation as the governing equation for 240 riverine bed-material sediment transport. Its implementation within MOSART-sediment is outlined below.
where is the total resistance coefficient [-]., The Engelund-Hansen relation considers Cf as the sum of skin friction and form drag. One can use the hydraulic resistance relation of Engelund-Hansen to compute for water depth and then Cf, given water 250 discharge, channel slope, and D50. However, Tthe computational procedure will require an iterative technique such as Newton-Raphson. More details can be found in Garcia (2008). However,Here we used the Manning's friction to compute Cf in order to avoid the extra computational expense and keep consistencyt with the routing computation in MOSART-water. As such, Cf and can be calculated from channel Manning's roughness coefficient, , [s/m 1/3 ], and channel hydraulic radius (in this study approximated as channel water depth ℎ), 255 where is the density of water [10 3 kg/m 3 ]; is the bed shear stress [Pa] and = ℎ in steady and uniform flows; is the local channel bed slope [-]. * is the dimensionless shear stress (the Shields number), given as * = 50 = 2 50 (4) 260 where is the channel velocity [m/s], and is estimated using Manning's equation (Manning, 1891).
where ℎ is the slope of the energy grade line or energy gradient [-]. In the kinematic wave routing method, ℎ is taken as the local channel slope . In the diffusion wave routing method, ℎ is estimated based on the difference in the simulated water surface elevations between the current and downstream channels. 265 Note that the Engelund-Hansen equation estimates the bed-material load that is the sum of bedload and suspended bedmaterial load. We adopt the transport-mode-parameter approach (Greimann et al., 2008) to separate the suspended bedmaterial load from the bed-material load computed by the Engelund-Hansen equation. The transport-mode-parameter approach specifies how much of a sediment size class is transported as bedload, suspended bed-material load, or mixed load following the "allocation coefficient" concept (Holly & Rahuel, 1990). Using this approach, the portion of suspended sediment 270 in the bed-material load is given asIt is given as: where f (0 ≤ f ≤ 1) is the transport-mode-parameter, i.e., the portion of suspended sediment in the bed-material load. is the suspension parameter and can be calculated as  Several formulas are available for sediment settling velocity, such as those named after Stokes (Zhang & Xie, 1993), Zhang (Zhang & Xie, 1993), Van Rijn (Rijn, 1989), and Cheng (Cheng, 1997), respectively. Cheng (1997) evaluated these formulas and suggested that Cheng's formula has the highest degree of prediction accuracy. It applies to a wide range of 280 Reynolds numbers from the Stokes flow to the turbulent regime (Cheng, 1997). Wu and Wang (W. Wu & Wang, 2006) proposed a new formula and suggested that this new formula has a highest accuracy among several formula, followed by Cheng's formula. However, Cheng's formula is mathematically much simpler than that of Wu and Wang, and more suitable for large-scale applications. Thus, Cheng's formula is adopted in this study.
is the kinematic viscosity of water at 20 ºC (=1.036×10 -6 m 2 s -1 ), and * is a dimensionless particle diameter, which is defined as: * = ( 2 ) 1/3 50 (9) 290 The Engelund-Hansen equation assumes that sediment transport is always on capacity (equilibrium transport) Engelund & Hansen, 1967;Naito et al., 2019) This simplification is not unreasonable for large-scale river modeling since the travel distance from one reservoir to another is typically more than 10km. Such a spatial scale enables the sediment-laden flow to sufficiently interact with channel bed to reach transport capacity .
For comparison with the Engelund-Hansen equation, we also use another widely used formula to directly estimate the 295 suspended sediment discharge rate, named after Wu (W. Wu et al., 2000) and well verified in multiple studies (Fang et al., 2017;Rousseau et al., 2016;Su et al., 2013;Sun et al., 2015;W. Wu et al., 2005). Wu's equation is given as: * = 0.0000262[( * * − 1) ] 1.74 when * > * where * is the dimensionless sediment flux for suspended load; * is the dimensionless critical shear stress above which sediment is mobilized (the default value is 0.0386). 300

Reservoir Sediment Processes
The reservoir sediment module incorporates both direct and indirect reservoir effects on sediment. The direct effect is the trapping of sediment within the reservoirs, preventing it from moving downstream. It applies to both wash load and bedmaterial load. The indirect effect is mainly through the reservoir regulation function simulated by MOSART-water, which modifies the hydraulic conditions, i.e., flow velocity, depth, and surface hydraulic gradient, in the downstream rivers of 305 reservoirs. Therefore, it affects only bed-material load. Note that in this study we do not explicitly account for routing of water or suspended sediment within reservoirs for simplicity.
Reservoir trapping of sediment is mainly through fundamentally slowing down of the river flow upstream of dams. The slowing down of river flow in the reservoir essentially sets the transport capacity of sediment to near zero, leading to the settling of suspended sediment in the water column. Thus, the deposition rate (trapping rate) of suspended sediment is largely 310 dependent on the settling time of suspended sediment in the reservoir. To this end, Brune (1953) developed an empirical formula for reservoir trapping efficiency and validated it for the U.S. reservoirs. It was later adopted globally (Syvitski et al., , 2005Vörösmarty et al., 2003). Here we also adopt Brune's formula to estimate individual reservoir trapping efficiency, : where Δ is the increase of local water residence time due to the reservoir [years], estimated as the effective reservoir storage capacity divided by the mean annual inflow from the reservoir's upstream. Other empirical formulas for reservoir trapping efficiency, e.g., Lewis et al. (2013), were also proposed with more complexity but have not been tested at the regional or global scales. It is worth notingNote that the sediment trapping in the reservoir resulted from the settling of suspended sediment in a nearly-still water body, and it is thus applied to all the suspended sediment regardless of particle size. Brune 320 (1953) and Lewis et al. (2013) assumed that wash load (silt and clay particles) and the suspended bed-material load (fine sand particles) were trapped by reservoirs in a similar way and did not differentiate them. Here we adopt the same assumption for simplicity.
gives the fraction of suspended load that is trapped in a reservoir, and the rest is released to downstream along with water.
The reservoir regulation effect on sediment is mainly by regulating river discharge, e.g., reducing high flow during a wet 325 season and increasing low flow during a dry season. This flow regulation will modify riverine hydraulic conditions such as flow velocity and channel water depth, thus changing bed-material load transport capacity. Besides reservoir regulation, another water management activity represented in MOSART-water is surface water withdrawal which diverts channel water to somewhere else via artificial channels or pipes. For simplicity, we assume that surface water withdrawal reduces only water storage in channels but not sediment. 330

Study Area
MOSART-sediment is designed to be applicable at the regional to global scales. In this study, we apply it over the contiguous United States (CONUS) for demonstration and validation, as shown in Figure 2. This study area is smaller than our previous study (Li, Leung, Tesfa, et al., 2015), where we applied ELM-MOSART driven by the atmospheric forcing data 335 from the North American Land Data Assimilation System Data Phase 2 (NLDAS2) over central North America (Mitchell et al., 2004).

Inputs and Parameters
To support the MOSART-sediment application, we have derived several new parameters a priori over CONUS that are most critical to sediment transport, including sediment particle size, Manning's roughness coefficient, and channel geometry.
These new parameters are meaningful mostly within CONUS since they have been derived based on various national databases.
It is for this reason that our sediment modeling and analysis focus on CONUS. Note that, although the sediment modeling analysis focuses on CONUS, MOSART-water is still applied over the same larger domain as our previous study (Li, Leung, Tesfa, et al., 2015), which includes the Canadian portion of the Columbia River basin. Thus other MOSART-water parameters and inputs are adopted from this previous study. 345 The median bed-material sediment particle size, 50 , is the most important parameter for the Engelund-Hansen formula.
In an earlier study (Abeshu et al., 2021), we derived a national map of the median bed-material sediment particle size ( 50 ) over CONUS based on: 1) the observed sediment particle size data from USGS and the U.S. Army Corps of EngineeringEngineers, 2) the comprehensive channel and catchment attributes from the National Hydrography Dataset Plus (NHDplus) (McKay et al., 2012;Schwarz et al., 2018), and 3) a predictive model established using state-of-the-art machine 350 learning techniques (Chen & Guestrin, 2016). The 50 map from Abeshu et al. (2021) is in a vector format, i.e., each flowline in NHDplus is assigned a value of the median bed-material sediment particle size. In this study, the vector-format map is converted to the grid-based format at the 1/8 th -degree resolution to be consistent with the other parameters, as shown in Figure   2a.
The channel geometry data, such as bankfull channel width and depth, are shown in Figures 2b and 2c. They are derived 355 following a data-driven hydraulic geometry study (Bieger et al., 2015), which provided the empirical formulas for the whole contiguous U.SCONUS. Manning's roughness coefficient for river bed is derived based on an empirical formula estimating Manning's roughness coefficient as a linear function of channel water depth distribution within a river basin (Getirana et al., 2014;Luo et al., 2017), as shown in Figure 2d. Channel bed slope (Figure 2e) is essentially from the medium resolution NHDplus. In NHDplus, each flowline represents one river segment and is provided a channel slope value (along with many 360 other channel attributes). We convert these NHDplus channel slope values from a vector format to a grid-based format here.
First, the NHDplus flowlines are georeferenced to the MOSART lat/lon grids at the 1/8 th -degree resolution in two steps: 1) local matching and 2) neighborhood matching. For the local matching, we adjust the geographic locations of the downstream ends of the NHDplus flowlines to the MOSART grids center. The drainage areas associated with the NHDplus flowlines (in a vector-based river network) are then compared against those associated with the corresponding MOSART grids (in a grid-365 based river network). The neighborhood matching is performed only for those grids where the margin of error is significant (>20% of NHDplus drainage area). Neighborhood matching is to re-adjust the flowlines from the current MOSART grid cell to one of the 8-neighboring cells where the NHDplus-MOSART drainage area difference is minimal. Then we associate each MOSART-water grid with a single NHDplus flowline. Hence, in the case of multiple flowlines per grid, which is very common, we choose the flowline that has a drainage area closest to the MOSART grid. After the above procedure, each MOSART grid 370 is associated with a flow line. The channel slope of each MOSART grid is then taken as the channel slope value from the corresponding flow line. Reservoir storage capacity (Figure 2f) is already available from our previous study (Li, Leung, Tesfa, et al., 2015) but also critical for sediment modeling, so we display them here for completeness.
The observed monthly suspended sediment discharge and streamflow from the United States Geological Survey (USGS) river gauges are used to validate the model performance. We choose to validate the model at the monthly scale instead of finer 375 scales due to two considerations: 1) The biases in the streamflow simulations by ESMs and even other land surface models are already noticeable at the monthly scale (Getirana et al., 2014(Getirana et al., , 2021Li, Leung, Getirana, et al., 2015;Mizukami et al., 2021).
2) The sampling of riverine sediment by USGS is not continuous like the streamflow measurement (taken automatically at an hourly or even shorter time step), but rather infrequently due to the required field visits and subsequent sample analysis (Groten & Johnson, 2018). There is thus considerable uncertainty in the daily suspended sediment data from USGS. Figure 3a shows 380 the locations of 35 selected USGS river gauges along with the river networks (displayed in the form of accumulated upstream drainage area at 1/8 th -degree resolution). Table 1 also provides more details of these selected gauges. These 35 USGS gauges are selected based on three criteria: 1) a no-more-than-20% difference between the actual upstream drainage area (reported by USGS) and the area represented within a grid-based river network, 2) at least five years of both monthly streamflow and suspended sediment observations in our study period 1990-2012 (see Section 4.1 for the reasons of choosing this period), and 385 3) relative bias of MOSART-water streamflow simulation (using the full modeling capacity) no more than 50% in 1990-2012.
The last criterion is to control the impacts of runoff and streamflow simulation uncertainty to a certain level. There are 63 gauges meeting the first two criteria. When considering the third criterion, the number drops from 63 to 35, suggesting that there is still large uncertainty in the runoff and streamflow simulation in E3SM. In fact, the uncertainty in runoff and streamflow simulation remains a challenge in Earth System modeling (Lawrence et al., 2019) and fully addressing it is beyond 390 the scope of this study. Among the selected gauges, their upstream drainage areas vary between 1225km 2 to 1850000km 2 , covering a range of spatial scales. Three representative USGS gauges (shown as the red stars in Figure 3b) are chosen out of the 35 gauges for more detailed analysis, including USGS05586100 Illinois River at Valley City in Illinois (ILV), USGS06807000 Missouri River at Nebraska City in Nebraska (MON), and USGS07022000 Mississippi River at Thebes in Illinois (MST). The choice of these three gauges 395 is based on two considerations: 1) these gauges are subject to different levels of reservoir effects, and 2) the drainage areas of these gauges are relatively larger than the others, hence more representing large rivers. The upstream drainage areas for ILV, MON, and MST are 6.93×10 5 km 2 , 1.06×10 6 km 2 , and 1.85×10 6 km 2 , respectively. In this study, we use Δ , the increase of average river water residence time due to all the reservoirs in the upstream areas of a lat/lon grid, to approximate the total reservoir effects including flow regulation and trapping, i.e., the larger the Δ value, the stronger the reservoir effects. 400 Note that Δ is different from Δ in Eqn. (11). For any lat/lon grid, Δ is calculated as ratio of the sum of storage capacity of all the reservoirs in its upstream areas over the long-term average inflow to this grid. Δ closely corresponds to the reservoir regulation index, the ratio of total upstream reservoir storage capacity to the annual streamflow volume (annual mean flow multiplied by the time of a year), which has been used as an indicator of reservoir flow regulation effects . As shown in Figure 3b, MON and ILV are subject to strong and weak reservoir effects, respectively. MST is 405 downstream of MON and ILV and subject to moderate reservoir effects.

Numerical Experiments
The simulation period is 1979-2012 due to the availability of model inputs, particularly the monthly water demand data that were derived in our previous study (Li, Leung, Tesfa, et al., 2015). ELM-Erosion is applied to 1979-2012 driven by the NLDAS2 atmospheric forcing at the 1/8 th -degree resolution to generate both daily sediment yield and daily runoff time series 410 at each 1/8 th -degree grid as the inputs for MOSART-sediment. All the soil erosion parameters are adopted from Tan et al. ( 2021). Note that Tan et al. (2021) accounted for both climate and land-use change impacts on soil erosion. In this study, we do not consider land-use change, i.e., assuming static landuse. Using the daily inputs from ELM-Erosion, we run MOSARTsediment in 1979-2012 at a daily time step and a 1/8 th -degree spatial resolution. We consider 1979-1989 as the model spin-up period and carry out the model validation and analysis in 1990-2012. 415 To isolate the relative contributions of reservoir regulation and trapping processes to suspended sediment transport, three major numerical simulations are designed: 1) sim_nat, MOSART-sediment is run under natural river conditions only, i.e., no reservoir regulation or trapping; 2) sim_wm_only, MOSART-sediment is run with the water management option, but turning off reservoir trapping; 3) sim_wm_trapping, MOSART-sediment is run with the water management option to include both flow regulation and reservoir trapping effects. The individual effects of reservoir regulation and trapping can be estimated by 420 comparing sim_nat with sim_wm_only and comparing the sim_wm_only and sim_wm_trapping, respectively. Each simulation is driven by the same sediment yield time series for the period of 1979-2012. Note that the suspended bed-material load in these three major experiments is estimated using the Engelund-Hansen equation, and main channel routing is simulated using the diffusion wave routing method.
We also perform additional numerical experiments to understand the sensitivity of modeling results to the key model 425 parameters, compare the Engelund-Hansen equation with the Wu formula, and compare the diffusion wave routing method with the kinematic wave routing method. More details of these additional numerical experiments will be provided later. Note that in all these numerical experiments, no parameter calibration is conducted for the riverine component. All related parameters have been estimated a priori.

Model Validation
MOSART-sediment captures the significant spatial heterogeneity in total suspended load, wash load, and suspended bedmaterial load, respectively, as shown in Figure 4 where the maps are generated from the sim_wm_trapping simulation results.
The spatial patterns of total suspended load are determined by the patterns of wash load, suspended bed-material load, and their relative dominance at different locations. It appears that suspended bed-material load dominates in the western U.S., e.g., 435 the Columbia and Colorado River basins. Wash load dominates in east of the Rocky Mountains, e.g., the Mississippi River basin, and coastal areas of northeastern U.S. We use the percentage of wash load in the total suspended load to illustrate the relative dominance of wash load at each lat/lon grid, as shown in Figure 4e. In this study, the model simulated spatial average of the wash load percentage is 62.2%, with higher values in the central and eastern U.S. and coastal western U.S. The dominance of wash load in the simulated suspended load is consistent with previous data-driven studies. For example, Sadeghi 440 and Singh (2017) analyzed the observed suspended load at 24 USGS gauges distributed throughout CONUS and suggested that wash load dominate in the majority of these gauges, i.e., the percentage of wash load in the total suspended load is 79.085 ± 11.343% over these gauges.
The spatial heterogeneity of simulated wash load is controlled mainly by that of sediment yield and the spatial distribution of reservoirs. More sediment yield can be from more humid or mountainous areas since it is triggered primarily by surface 445 runoff and raindrops. Reservoir trapping significantly reduces the amount of wash load carried to the rivers downstream of reservoirs. The supply of sand sediment is mostly from the in-channel process, as a function of grain size, channel geometry, and hydraulic conditions such as the shear stress and channel velocity. Roughly, the sand discharge is large when the channel water depth is deep, the channel bed slope is steep, or Manning's roughness coefficient is small. Reservoir regulation reduces the high flow during the wet season and then reduces suspended bed-material load by reducing channel water depth and 450 velocity in the rivers downstream of reservoirs. However it enhances the low flow during the dry season and increases suspended bed-material load in the rivers downstream of reservoirs, which is consistent with the findings reported from literature (Naito et al., 2019;Nittrouer & Viparelli, 2014). Reservoir trapping appears to have little effect on the simulated suspended bed-material load. Reservoirs prevent the upstream suspended bed-material load from going downstream. But the reduced suspended bed-material load is supplemented by the local channel erosion in downstream rivers. Recall the assumption 455 underpinning the Engelund-Hansen equation that the local supply of suspended bed-material load is unlimited.
MOSART-sediment reasonably well captures the long-term average suspended sediment discharge values across CONUS, as shown in Figure 5 for comparison between the model simulated and observed annual mean values in 1990-2012 at multiple USGS gauges. Here the comparison is also shown between the sim_nat and sim_wm_trapping simulations. The model performs very well in simulating the long-term streamflow over the selected 35 gauges. The difference in streamflow between 460 sim_nat and sim_wm_trapping is negligible, which is reasonable because 1) reservoir trapping does not affect streamflow at all and 2) reservoir regulation changes mostly intra-annual variability of streamflow and sometimes inter-annual variability, but the long-term mean streamflow is largely unaffected. For suspended sediment, the simulated annual mean values are overall close to those observed, as indicated by the close alignment of the dots with the 1:1 line, particularly for sim_wm_trapping. Figure 5 to help visualize the overall patterns. The R 2 values between the simulations and 465 observations are 1.0 for the streamflow from both sim_nat and sim_wm_trapping, and 0.6 for the suspended sediment from both sim_nat and sim_wm_trapping respectively. WBMsed was validated over CONUS at 11 USGS gauges in 199711 USGS gauges in -2007 with an overall satisfactory performance (see Fig. 7 in Cohen et al., 2013). WBMsed nevertheless overestimated the long-term average suspended load over the four largest gauges on the Mississippi River (here "largest" means the gauges with the largest suspended load observations), whilst our model seems to have less overestimation over the largest gauges, e.g., Mississippi 470

Fitted lines are included in
River at Thebes (USGS ID 07022000). Comparing Fig. 7 in Cohen et al. (2013) with Fig. 5 in this study, one can infer that our model performance is comparable to WBMsed.
We further use the Kling-Gupta coefficient (KGE) (Gupta et al., 2009) and the normalized-root-mean-square (NRMSE) to evaluate the model simulated long-term average across multiple gauges. KGE is a comprehensive metric that integrates the effectiveness of previous extensively used metrics such as root-mean-square-error and the Nash-Sutcliffe coefficient (NSE). 475 It thus captures both the magnitude and timing (or phase) difference between two temporal or spatial series. KGE ranges from negative infinity to 1.0, with 1.0 indicating perfect model performance (Gupta et al., 2009). It has been suggested that a KGE value higher than -0.41 implies that a model performs better than the long-term mean benchmark, i.e., equivalent to NSE=0 (Knoben et al., 2019). NRMSE is calculated as the ratio of RMSE over the long-term mean of a temporal or spatial series. For streamflow, KGE increases from 0.89 to 0.96 and NRMSE decreases from 0.24 to 0.12 by including the reservoir effects, i.e., 480 from sim_nat to sim_wm_trapping. For suspended sediment, KGE increases from -0.40 to 0.69, and NRMSE decreases from 2.35 to 0.74 by including the reservoir effects. Overall, Figure 5 suggests a very good model performance in capturing the long-term suspended sediment balance across different spatial scales. It also suggests that adding the reservoir effects improves the simulated long-term suspended sediment loads more significantly than the simulated long-term streamflow over the selected gauges. 485 We further evaluate the model performance at the monthly scale using the KGE values between the simulated and observed monthly time series. We roughly classify the KGE values at the 35 gauges into three categories: Category I, KGE≤-0.41, poor model performance; Category II, -0.41<KGE≤0, decent model performance; Category III, KGE>0, good model performance.  Table 1 lists the specific KGE and NRMSE values between the monthly simulated and observed streamflow and suspended sediment time series at each gauge based on the sim_wm_trapping simulation. Note that the availability of monthly observations varies between different gauges. Overall, MOSART-sediment has a decent or good performance at over 25 of the 35 gauges, indicating that the model is able to capture reasonably well the intra-annual variability of suspended sediment in large river systems. 495 Compared to streamflow, model biases in the suspended sediment are noticeably larger. Figure 5 suggests that MOSARTsediment overall overestimates annual mean suspended sediment load in many gauges. Figure 6 also suggests poor model performance at some gauges. There are several possible reasons for these biases: 1) We have only included ~1800 large reservoirs (storage capacity no less than 0.1km 3 ) with critical information (storage capacity, location, major purposes, etc.) publicly available from the Global Reservoir and Dam Database (GRanD) (Lehner et al., 2011). Relatively small reservoirs in 500 CONUS are not included in the model inputs due to a lack of necessary information, which will be addressed in future work.
2) There are noticeable biases already in the runoff and streamflow simulation at the monthly or shorter time scales, as shown in Figure 7, which propagate to the suspended sediment simulations. 3) There are uncertainties in the model structure and model parameters such as sediment particle size, channel slope, etc., which we will discuss next.

Model Structure Uncertainty
MOSART-sediment simulates hillslope, riverine, and reservoir processes. Here we focus on the riverine and reservoir processes. Uncertainties in the hillslope process (e.g., soil erosion) were discussed in previous studies (Tan et al., 2018(Tan et al., , 2021. For the simulated riverine processes, we mainly analyze the uncertainties in the sediment governing equations and routing methods. The accuracy of suspended sediment simulation could be affected by the choice of sediment transport equations and 510 simulation of riverine hydraulic conditions. Many sediment transport equations, like such as Engelund-Hansen and Wu's equations, are well-established in the laboratory or single-reach scale. The implementations of such sediment transport equations at a regional scale for morphology simulation purposes are also well established (Naito et al., 2019;Nittrouer et al., 2012). Their applications in simulating suspended sediment at the continental scale are nevertheless rarely reported. Note our   Figure 7b includes all of the 35 gauges. Over some gauges, the Engelund-Hansen equation (plus the transport-mode-parameter approach) performs better than Wu's formula, particularly as indicated by the significantly decreased KGE and increased NRMSE values over some gauges. The long-term KGE and NRMSE values 525 between the simulated and observed suspended load over the selected gauges are 0.29 and 1.76 respectively when using Wu's formula and diffusion wave method, both are worse than those using the Engelund-Hansen equation and diffusion wave method. Therefore, the Engelund-Hansen equation (plus the transport-mode-parameter approach) overall outperforms Wu's formula.
As for the routing methods, it appears that the difference between the sediment simulations using the diffusion wave 530 routing method and those using the kinematic wave routing method is overall less than the difference between the Engelund-Hansen and Wu's formulas. The difference between the monthly streamflow simulations using different routing methods is even more negligible (Figure not shown). At the 35 selected stations located mostly upstream of river outlets, the backwater effects are likely not dominant. In flat locations such as close to the river mouths, however, the backwater effects are expected to be significant, and as such using the diffusion wave routing method should lead to better performance in not only the 535 streamflow but also the sediment simulations.
For the reservoir process, we analyze the uncertainties in the reservoir regulation and trapping processes. Figure 8 illustrates the total impacts of the two reservoir processes by contrasting the sim_nat and sim_wm_trapping simulations results at the three representative gauges subject to different levels of reservoir effects. At ILV, whose upstream drainage area is subject to relatively small reservoir effects, there is a negligible difference between the simulated streamflow from sim_nat 540 and sim_wm_trapping, indicating relatively minor water management activities in the upstream areas due to a relatively smaller number of reservoirs and low water demand. At MON, its upstream drainage area is subject to strong reservoir effects.
The simulated streamflow is significantly reduced from sim_nat to sim_wm_trapping particularly over the high-flow periods in the summer, mainly due to the intensive surface water extraction in the upstream drainage area and the strong reservoir regulation effect. At MST, whose upstream drainage area is subject to moderate reservoir effects, the simulated streamflow is 545 also noticeably reduced from sim_nat to sim_wm_trapping during the high-flow periods. At all these gauges, the simulated suspended load is significantly reduced from sim_nat to sim_wm_trapping by reservoir trapping. In the upstream area of ILV, there are only a few reservoirs simulated. These reservoirs, however, can still trap a large portion of suspended sediment, suggesting that the effects of reservoirs on suspended sediment discharge are likely more significant than those on streamflow. Figure 9 displays the change of KGE for the monthly streamflow and suspended load at all selected USGS gauges from 550 sim_nat to sim_wm_trapping. After adding the reservoir effects (sim_wm_trapping), the model performance is significantly improved (KGE increases more than 0.05) at 5 and 23 of the 35 USGS gauges for streamflow and suspended load simulations, respectively, and worsened at 12 and 2 USGS gauges, respectively. For the remaining gauges, the change of model performance is negligible. For streamflow, the model performance is already quite good even for sim_nat in terms of the monthly streamflow simulation (see Fig. 5 and Table 1). The uncertainties in the runoff generation, routing structures and 555 parameters appear to overshadow the benefits of adding water management, hence leading to worsened model performance at 12 gauges. The improvement of sediment simulation by adding both reservoir regulation (as part of water management) and trapping is obviously more dominant. Figures 8 and 9 only show the combined effects of reservoir regulation and trapping, as interpreting their individual effects is not straightforward.
The net effects of reservoir regulation and trapping are further examined in Figure 10 by comparing the mean monthly 560 streamflow and suspended sediment discharges among the sim_nat (blue lines), sim_wm_only (green lines), sim_wm_trapping (black lines) simulations. By comparing the blue and green lines, one can see the effects of flow regulation caused by water management. Overall, the net effects of reservoir operations on streamflow are to reduce the high-flow in the wet season (May to August) and increase the low-flow in the dry season, and these effects are relatively minor on the flow discharge at ILV, significant at MON, and moderate at MST. Reservoir regulation also reduces the suspended sediment 565 discharge, but the effects are only noticeable during the high-flow period and negligible over the low-flow period when the sediment load is already low. By comparing the green and black lines in Figure 10, one can see the net effects of reservoir trapping. Reservoir trapping largely reduces the magnitude of suspended sediment discharge across all three gauges, including ILV where the reservoir regulation effect seems minor. After adding reservoir trapping, the simulated sediment discharge is much closer to the observations at all three gauges. 570 The effect of reservoir regulation on streamflow can be propagated from one reservoir to another downstream in a cascading manner. This is because reservoirs do not store water in the long term but only hold it for a short period (i.e., a few months) and then release it downstream. The effects of reservoir trapping may propagate from upstream to downstream for wash load but not for suspended bed-material load. For wash load, reservoirs intercept and keep most of it without releasing it to the downstream reservoirs. Hence, the trapping effect propagates downstream. For suspended bed-material load, even if 575 reservoirs trap most of it from upstream inflow, the downstream river channel beds and in-channel bars will be eroded to supply suspended bed-material load to maintain a dynamic equilibrium state between the in-channel erosion and deposition.
This way, the transport rate of bed-material load in the river downstream of a reservoir is only limited by the local hydraulic conditions, i.e., according to the Engelund-Hansen equation.
Although including reservoir trapping improves the simulation of suspended sediment discharge, there are still some 580 discrepancies between the simulations and observations, such as the underestimation of sediment discharge in January-April (at ILV and MST) and overestimation in July-October (at all three gauges). These discrepancies between the simulated and observed seasonal variations in the suspended sediment discharge are likely not due to the streamflow biases. Besides possible uncertainties in the sediment model parameters, another probable cause of the model biases in sediment discharge is that the Brune formula adopted in this study estimates the long-term average reservoir trapping coefficient without explicitly 585 accounting for the seasonal variation in the trapping efficiency (Brune, 1953).

Model Parameter Uncertainty
Here we present model parameter uncertainty analysis focusing on three important and representative parameters: median bed-material sediment particle size, 50 , channel slope, , and Manning's Roughness coefficient for channels, . is involved in multiple riverine sediment processes (e.g., Eqn. (1), (4), (8), and (9)). and affect the sediment processes both 590 directly and indirectly (e.g., Eqn. (3)-(5)). Other parameters are not included here for simplicity. Figure 11 shows the changes in the simulated long-term average suspended load in 1990-2012 (based on sim_wm_trapping) after increasing or decreasing , , and by 10% respectively. The spatial patterns of these changes are largely controlled by the spatial patterns of the simulated wash load and suspended bed-material load (see Figure 4). , , and are mostly affecting the riverine hydraulic conditions and not much wash load transport. The changes in the total 595 suspended load are thus only significant over places dominated by the suspended bed-material load. For places dominated by wash load, the changes in the total suspended load are overall small or even negligible. In fact, the median percentage changes in the suspended load are mainly in the range of -2%~2% for all three parameters (Figure 13). Comparing among the three parameters, the simulation of the suspended bed-material load is most sensitive to the changes in , followed by , and least sensitive to the changes in . 600

Impacts of reservoirs on suspended sediment discharge to the coasts
Despite the uncertainties, MOSART-sediment is a useful tool to shed light on some interesting questions. For example, how will reservoirs affect fresh water and suspended sediment discharges to the coasts? Figure 13 shows the model-simulated impacts of water management on freshwater discharge (Figure 13a) and the effects of reservoirs on suspended sediment discharge (Figure 13b) to the coasts. Note that in Figure 13b, we focus on wash load only instead of total suspended load, i.e., 605 excluding the suspended bed-material load. There are three reasons: 1) Wash load consists of very fine particles and is a major carrier of carbon and nutrients from land into rivers and then to the coasts. 2) Wash load is dominant in most rivers. 3) Wash load is subject to the impacts of both upstream drainage areas (e.g., land use change, water management, precipitation) and local channel conditions (e.g., channel slope and, 50 ), whilst suspended bed-material load is mostly dominated by the local channel conditions. 610 According to our simulations, water management moderately reduces the freshwater discharge to the coasts, particularly for the smaller rivers on the southwest coasts ( Figure 13a). The major cause of reduction in freshwater discharge is surface water extraction (which is mainly controlled by water demand intensity) since reservoir regulation normally does not reduce streamflow on a decadal scale (e.g., in the period of 1990-2012). The reduction of freshwater discharge to the southwest coasts is thus consistent with the high irrigation water demand over the western U.S., e.g., in the Colorado Rriver basin. 615 The delivery of suspended sediment from the continental U.S. to the coasts is reduced more significantly than streamflow ( Figure 13b). According to our simulation results, the reduction of suspended load by reservoirs is high (e.g., over 95%) in some river basins such as Colorado and Brazos River but medium or low in other rivers such as the Mississippi and Columbia River. The former rivers have a large proportion of wash load relative to the total suspended load (indicated by the colors of the circles). As such, the reservoir trapping effect plays a vital role. In many small rivers close to the coasts but not subject to 620 reservoir effects, most of their suspended sediment discharge is delivered to the coasts without being trapped during the transport. We do not find any notable relationship between the percentage reductions in streamflow vs. those in suspended sediment discharge at the river mouths. For the former, surface water extraction is the major cause; whilst for the latter, reservoir trapping is the major cause. These two major causes do not directly interact with each other, suggesting that the reduction of suspended sediment discharge may not be closely linked to the reduction of streamflow. These inferences, 625 however, are made based on the simulation results, and further investigations using more comprehensive observational datasets and modeling approaches are needed in the future. Moreover, we caution that the wash load estimates (with sim_wm_trapping) do not account for the effects of small reservoirs, and are likely higher than the real conditions. Nevertheless, our estimates of wash load from large rivers to the coasts are comparable with previously reported numbers.
Within the conterminous United States, the greatest amount of sediment is transported by rivers within the Gulf of Mexico 630 drainage system. Meade and Moody (2010) reported that in 1987-2003, the Mississippi River system on average discharged 172 Million-tons/year total suspended load to the coast, which is moderately less than our wash load estimation of 228.23 Million-tons/year from sim_wm_trapping (used in all the comparisons below unless stated otherwise). Holeman (1968) estimated that the Columbia River discharged about 9 Million tons of total suspended load per year to the Pacific Ocean, and 15 years later, this magnitude was reported to be reduced to 8 Million tons/year (Meade & Milliman, 1983). Wise et al. (2007)  be 16 Million-tons/year, a much-reduced value. Our wash load estimate is 1.55 Million-tons/year, even further reduced, which is possible given that there has been a decreasing trend of suspended sediment discharge in the U.S. river systems over the past few decades due to improved soil erosion controls (Meade & Moody, 2010).

Summary and Conclusion
In this study, a new large-scale suspended sediment module has been developed on top of MOSART within E3SM. Both 645 reservoir trapping and regulation mechanisms have been incorporated. Validation of the model simulations against the observations from 35 USGS gauges located across the contiguous United States shows reasonable skill in simulating the longterm average suspended sediment discharge (KGE=0.69, NRMSE=0.74) and monthly suspended load (KGE>0 in 18 of the 35 gauges). We show that adding the reservoir effects improves the model simulations, particularly at the annual and monthly scales. As such, our model performance is comparable to WBMsed without any calibration. We carry out the model uncertainty 650 analysis in terms of the model structure and parameter uncertainties. We further show that reservoir trapping exerts dominant effects (over reservoir regulation) on suspended sediment discharge through large river systems to the coasts.
The analyses help identify several opportunities to further advance large-scale suspended sediment modeling, including but not limited to: (1) improving the runoff and streamflow simulation within the Earth system modeling framework; (2) effectively accounting for the intra-or inter-annual variations in the reservoir trapping efficiency over multiple reservoirs at 655 the regional or global scales; (3) accounting for more realistic reservoir stage-volume relationships and their impacts on reservoir trapping; (43) including the effects of small or navigational reservoirs on streamflow and sediment simulations; (54) improving the accuracy of channel bed slope estimation; (6) accounting for sediment exchanges between channel and floodplain more realistically. Moreover, our separate treatment of wash load and bed-material load will make it easier to isolate various climate and human impacts on riverine sediment. For instance, water management such as local channel water 660 withdrawal and reservoir operation will primarily affect bed-material load by modifying hydraulic conditions. Land management such as land-use change and farming activities will affect wash load by changing topsoil erosional conditions. Lastly, for the sediment particle size, we have developed adopted a spatial map of median bed-material sediment particle size to support the parameterization of the sediment governing equations but there is still room to enhance the accuracy and availability of the sediment particle data through improved understanding and more data collection. Each of these opportunities 665 above requires a chain of efforts to achieve and is thus beyond the scope of this work.
While future opportunities exist to further improve MOSART-sediment, this study presents a new process-based sediment modeling framework from a coupled earthEarth-human system perspective. It has the potential to enable investigations that were not possible on the suspended sediment delivery to the coasts under various disturbances induced by variability and changes in climate and human activities. For example, the models can be used to explore how extreme events such as floods 670 and droughts and their long term changes may influence sediment delivery to the coast and hence the biogeochemistry of the coastal ecosystem, and the role of soil erosion and sediment transport in atmospheric CO2 through changes in the terrestrial, riverine, and ocean biogeochemistry. Finally and importantly, the new sediment module lays the foundation for modeling transportation and transformation of particulate C, N, and P through the river-ocean continuum to help close the global C, N, and P cycles within Earth system models. 675

Data Availability
The data used in this study will be available upon request to the corresponding author. 680

Author Contribution
HL developed the code, carried out the simulations and wrote the manuscript draft. HL, ZT, HM, ZZ, SZ, SC all contributed to the methodology. HL, ZZ and LRL performed the analysis. GA, TZ and DX helped with the model inputs and observational data processing. HL and GA prepared the figures. All authors contributed to the writing and discussion.

Competing Interest 685
The authors declare that they have no conflict of interest.  List of Figures   Fig. 1 Modeling framework of MOSART-sediment. 1030       Fig. 11 Spatial patterns of sensitivity of median bed-material size (d50) (a, b), channel slope (rslp) (c, d), and Manning's roughness for channels (nr) (e, f). At each grid, the percentage change of suspended load is calculated as 100%*(S1-S0)/S0, where S1 and S0 are the average suspended load 1990-2012 simulated with the parameter changing by 10%, and the average suspended load simulated with the original parameters, respectively.  Figure 13 represents the lat/lon grid at the river mouth of a river system. The larger the dot size, the larger the freshwater or suspended sediment discharge simulated with Sim_wm_trapping. The color of each dot represents the percentage of freshwater or suspended sediment discharge that is reduced by water management or reservoirs, i.e., (sim_nat -Sim_wm_trapping)/sim_nat.    Figure 13 represents the lat/lon grid at the river mouth of a river system. The larger the dot size, the larger the freshwater or suspended 1195 sediment discharge simulated with Sim_wm_trapping. The color of each dot represents the percentage of freshwater or suspended sediment discharge that is reduced by water management or reservoirs, i.e., (sim_nat -Sim_wm_trapping)/sim_nat.